ECC-1258: find_nearest returns wrong point if searched location is at a grid point

src/grib_geography.c

@@ -4155,14 +4155,15 @@ double geographic_distance_spherical(double radius, double lon1, double lat1, do
     if (lat1 == lat2 && lon1 == lon2) {
         return 0.0; /* the two points are identical */
     }
-
     if (rlon1 >= 360) rlon1 -= 360.0;
     rlon1 = RADIAN(rlon1);
     if (rlon2 >= 360) rlon2 -= 360.0;
     rlon2 = RADIAN(rlon2);
 
     a = sin(rlat1) * sin(rlat2) + cos(rlat1) * cos(rlat2) * cos(rlon2 - rlon1);
-    DebugAssert(a >= -1 && a <= 1);
+    /* ECC-1258: sometimes 'a' can be very slightly outside the range [-1,1] */
+    if (a > 1.0) a = 1.0;
+    if (a < -1.0) a = -1.0;
 
     return radius * acos(a);
 }

tests/CMakeLists.txt

@@ -85,6 +85,7 @@ if( HAVE_BUILD_TOOLS )
         grib_ecc-1212
         grib_ecc-1230
         grib_ecc-1255
+        grib_ecc-1258
         bufr_ecc-1028
         bufr_ecc-1195
         bufr_json_samples

tests/grib_ecc-1255.sh

@@ -10,14 +10,15 @@
 
 . ./include.sh
 set -u
-REDIRECT=/dev/null
+
 label="grib_ecc-1255-test"
 tempGrib=temp.$label.grib
 tempFilt1=temp.1.$label.filt
 tempFilt2=temp.2.$label.filt
-
 sample_grib2=$ECCODES_SAMPLES_PATH/GRIB2.tmpl
 
+# Non-EFAS GRIB
+# -----------------
 cat > $tempFilt1 <<EOF
  print "is_efas=[is_efas]";
 EOF
@@ -25,6 +26,8 @@ result=`${tools_dir}/grib_filter $tempFilt1 $sample_grib2`
 [ "$result" = "is_efas=0" ]
 
 
+# EFAS GRIB
+# -----------
 cat > $tempFilt2 <<EOF
  set setLocalDefinition = 1;
  set localDefinitionNumber = 41; # sets is_efas to 1
@@ -32,11 +35,9 @@ cat > $tempFilt2 <<EOF
 EOF
 ${tools_dir}/grib_filter -o $tempGrib $tempFilt2 $sample_grib2
 
-#${tools_dir}/grib_ls -p is_efas $tempGrib
-
 result=`${tools_dir}/grib_filter $tempFilt1 $tempGrib`
 [ "$result" = "is_efas=1" ]
 
 
+# Clean up
 rm -f $tempGrib $tempFilt1 $tempFilt2
-

tests/grib_ecc-1258.sh

@@ -0,0 +1,38 @@
+#!/bin/sh
+# (C) Copyright 2005- ECMWF.
+#
+# This software is licensed under the terms of the Apache Licence Version 2.0
+# which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+# 
+# In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
+# virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
+#
+
+. ./include.sh
+set -u
+
+label="grib_ecc-1258-test"
+tempFilt=temp.$label.filt
+tempGrib=temp.$label.grib
+tempOut=temp.$label.txt
+sample_grib1=$ECCODES_SAMPLES_PATH/GRIB1.tmpl
+
+cat > $tempFilt <<EOF
+ set Ni = 16;
+ set Nj = 21;
+ set latitudeOfFirstGridPoint= -29600;
+ set longitudeOfFirstGridPoint = 134000;
+ set latitudeOfLastGridPoint= -37600;
+ set longitudeOfLastGridPoint= 140000;
+ set iDirectionIncrement = 400;
+ set jDirectionIncrement= 400;
+ write;
+EOF
+
+${tools_dir}/grib_filter -o $tempGrib $tempFilt $sample_grib1
+${tools_dir}/grib_ls  -l -32.4,137.6  $tempGrib > $tempOut
+grep -q "Grid Point chosen #3 index=121" $tempOut
+grep -q "index=121 .*distance=0.0" $tempOut
+
+
+rm -f $tempFilt $tempGrib $tempOut