From 714415054d31ffc2523376aedcce7c59e1dbf825 Mon Sep 17 00:00:00 2001
From: Shahram Najm <masn@ecmwf.int>
Date: Thu, 28 Jul 2022 12:34:26 +0100
Subject: [PATCH] ECC-762 and ECC-1432

---
 data/ls.log                      | 12 ++++-----
 src/grib_nearest.c               |  5 +---
 src/grib_nearest_class_reduced.c | 45 +++++++++++++++++++++++++++++++-
 tests/grib_nearest_test.sh       | 16 ++++++++++++
 4 files changed, 67 insertions(+), 11 deletions(-)

diff --git a/data/ls.log b/data/ls.log
index 4e290e613..20d9557f4 100644
--- a/data/ls.log
+++ b/data/ls.log
@@ -75,17 +75,17 @@ index=1290
 index=1171
 index=1170
 reduced_gaussian_sub_area.grib1
-#3
-index=19985
-index=19984
+#1
 index=19665
+index=19985
 index=19664
+index=19984
 reduced_gaussian_sub_area.grib2
-#3
-index=19985
-index=19984
+#1
 index=19665
+index=19985
 index=19664
+index=19984
 reduced_gaussian_surface.grib1
 #4
 index=1291
diff --git a/src/grib_nearest.c b/src/grib_nearest.c
index a30245232..5a69622d5 100644
--- a/src/grib_nearest.c
+++ b/src/grib_nearest.c
@@ -342,10 +342,7 @@ int grib_nearest_find_generic(
         double lat1 = 0, lat2 = 0;     /* inlat will be between these */
         const double LAT_DELTA = 10.0; /* in degrees */
 
-        if (grib_is_missing(h, Ni_keyname, &ret)) {
-            grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", Ni_keyname);
-            return ret ? ret : GRIB_GEOCALCULUS_PROBLEM;
-        }
+        /* Note: If this is being called for a REDUCED grid, its Ni will be missing */
 
         if (grib_is_missing(h, Nj_keyname, &ret)) {
             grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", Nj_keyname);
diff --git a/src/grib_nearest_class_reduced.c b/src/grib_nearest_class_reduced.c
index 598985892..0b75b1167 100644
--- a/src/grib_nearest_class_reduced.c
+++ b/src/grib_nearest_class_reduced.c
@@ -130,17 +130,60 @@ static int init(grib_nearest* nearest, grib_handle* h, grib_arguments* args)
 
 typedef void (*get_reduced_row_proc)(long pl, double lon_first, double lon_last, long* npoints, long* ilon_first, long* ilon_last);
 
+static int find_global(grib_nearest* nearest, grib_handle* h,
+                double inlat, double inlon, unsigned long flags,
+                double* outlats, double* outlons, double* values,
+                double* distances, int* indexes, size_t* len);
+
 static int is_legacy(grib_handle* h)
 {
     long is_legacy = 0;
     return (grib_get_long(h, "legacyGaussSubarea", &is_legacy) == GRIB_SUCCESS && is_legacy == 1);
 }
 
-/* Old implementation in src/deprecated/grib_nearest_class_reduced.old */
+
 static int find(grib_nearest* nearest, grib_handle* h,
                 double inlat, double inlon, unsigned long flags,
                 double* outlats, double* outlons, double* values,
                 double* distances, int* indexes, size_t* len)
+{
+    int err = 0;
+    grib_nearest_reduced* self = (grib_nearest_reduced*)nearest;
+
+    if (self->global) {
+        err = find_global(nearest, h,
+                inlat, inlon, flags,
+                outlats, outlons, values,
+                distances, indexes, len);
+    }
+    else
+    {
+        /* ECC-762, ECC-1432: Use brute force generic algorithm
+         * for reduced grid subareas. Review in the future
+        */
+        int lons_count = 0; /*dummy*/
+
+        err = grib_nearest_find_generic(
+            nearest, h, inlat, inlon, flags,
+            self->values_key,
+            "Ni",
+            self->Nj,
+            &(self->lats),
+            &(self->lats_count),
+            &(self->lons),
+            &(lons_count),
+            &(self->distances),
+            outlats, outlons,
+            values, distances, indexes, len);
+    }
+    return err;
+}
+
+/* Old implementation in src/deprecated/grib_nearest_class_reduced.old */
+static int find_global(grib_nearest* nearest, grib_handle* h,
+                double inlat, double inlon, unsigned long flags,
+                double* outlats, double* outlons, double* values,
+                double* distances, int* indexes, size_t* len)
 {
     grib_nearest_reduced* self = (grib_nearest_reduced*)nearest;
     int ret = 0, kk = 0, ii = 0, jj = 0;
diff --git a/tests/grib_nearest_test.sh b/tests/grib_nearest_test.sh
index 4799d06aa..1ca9e379a 100755
--- a/tests/grib_nearest_test.sh
+++ b/tests/grib_nearest_test.sh
@@ -119,6 +119,22 @@ ${tools_dir}/grib_ls -l 0,0,1 $tempGrib > $temp
 grep -q "Grid Point chosen #2 index=480 " $temp
 
 
+# Subarea: Reduced Gaussian grid
+# ------------------------------
+input_grb=${data_dir}/reduced_gaussian_sub_area.grib1
+${tools_dir}/grib_get -F%3.6g -l 89.656,220,1 $input_grb > $temp
+grep -q "274.038" $temp
+${tools_dir}/grib_ls -l 89.656,220,1 $input_grb > $temp
+grep -q "Grid Point chosen #1 index=0 " $temp
+
+${tools_dir}/grib_ls -l 0.225,399.6,1 $input_grb > $temp
+grep -q "Grid Point chosen #1 index=53563 " $temp
+
+${tools_dir}/grib_get -F%3.6g -l 11.91,366.3,1 $input_grb > $temp
+grep -q "299.632" $temp
+${tools_dir}/grib_ls -l 11.91,366.3,1 $input_grb > $temp
+grep -q "Grid Point chosen #1 index=43089 " $temp
+
 
 # Clean up
 rm -f $temp $tempRef $tempGrib