ECC-499: function grib_nearest_find to get the indexes but not the values

examples/C/CMakeLists.txt

@@ -83,6 +83,7 @@ list( APPEND tests_extra
    grib_set_bitmap
    grib_list
    grib_get_data
+   grib_nearest
    grib_nearest_multiple
    set_missing
    bufr_attributes

examples/C/grib_nearest.c

@@ -21,7 +21,7 @@
 #include "eccodes.h"
 
 static void usage(const char* prog) {
-    printf("Usage: %s grib_file grib_file ...\n",prog);
+    printf("Usage: %s [-n] grib_file grib_file ...\n",prog);
     exit(1);
 }
 
@@ -31,6 +31,7 @@ int main(int argc, char** argv)
     long step=0;
     size_t nfiles;
     int i=0;
+    int no_values = 0; /* Default: decode the values */
     codes_fieldset* set=NULL;
     codes_handle* h=NULL;
     char param[20]={0,};
@@ -41,26 +42,35 @@ int main(int argc, char** argv)
     double distances[4]={0,};
     int indexes[4]={0,};
     char* order_by="param,step";
+    double* pValues = values;
 
     size_t size=4;
     double lat=-40,lon=15;
     int mode=0;
     int count;
+    int filearg_offset = 1;
     char** filenames;
     codes_nearest* nearest=NULL;
 
     if (argc < 2) usage(argv[0]);
-
     nfiles=argc-1;
+
+    if (strcmp(argv[1],"-n")==0) {
+        no_values = 1; /* No values decoding. Just indexes, distances etc */
+        nfiles -= 1;
+        pValues = NULL;
+        filearg_offset++;
+    }
+
     filenames=(char**)malloc(sizeof(char*)*nfiles);
     for (i=0;i<nfiles;i++)
-        filenames[i]=(char*)strdup(argv[i+1]);
+        filenames[i]=(char*)strdup(argv[i+filearg_offset]);
 
     set=codes_fieldset_new_from_files(0,filenames,nfiles,0,0,0,order_by,&err);
     CODES_CHECK(err,0);
 
-    printf("\nordering by %s\n",order_by);
-    printf("\n%d fields in the fieldset\n",codes_fieldset_count(set));
+    printf("ordering by %s\n",order_by);
+    printf("%d fields in the fieldset\n",codes_fieldset_count(set));
     printf("n,step,param\n");
 
     mode=CODES_NEAREST_SAME_GRID | CODES_NEAREST_SAME_POINT;
@@ -74,10 +84,16 @@ int main(int argc, char** argv)
         printf("%d %ld %s  ",count,step,param);
         if (!nearest) nearest=codes_grib_nearest_new(h,&err);
         CODES_CHECK(err,0);
-        CODES_CHECK(codes_grib_nearest_find(nearest,h,lat,lon,mode,lats,lons,values,distances,indexes,&size),0);
+        CODES_CHECK(codes_grib_nearest_find(nearest,h,lat,lon,mode,lats,lons,pValues,distances,indexes,&size),0);
+        if (no_values) {
+            printf("\nIdx\tlat\tlon\tdist\n");
+            for (i=0;i<4;i++) printf("%d\t%.2f\t%.2f\t%g\n",
+                    (int)indexes[i],lats[i],lons[i],distances[i]);
+        } else {
             printf("\nIdx\tlat\tlon\tdist\tval\n");
             for (i=0;i<4;i++) printf("%d\t%.2f\t%.2f\t%g\t%g\n",
                     (int)indexes[i],lats[i],lons[i],distances[i],values[i]);
+        }
         printf("\n");
 
         codes_handle_delete(h);

examples/C/grib_nearest.sh

@@ -0,0 +1,53 @@
+#!/bin/sh
+# Copyright 2005-2019 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
+
+label="grib_nearest"
+temp=$label.temp
+tempRef=$label.ref
+input_grb=${data_dir}/reduced_gaussian_pressure_level.grib1
+
+# Nearest with decoding the data values
+# --------------------------------------
+${examples_dir}/c_grib_nearest $input_grb > $temp
+cat > $tempRef <<EOF
+ordering by param,step
+1 fields in the fieldset
+n,step,param
+1 0 t  
+Idx	lat	lon	dist	val
+4839	-40.46	18.00	259.679	285.527
+4838	-40.46	15.00	51.5268	284.074
+4719	-37.67	18.00	366.445	289.621
+4718	-37.67	15.00	258.597	289.027
+
+EOF
+diff $tempRef $temp
+
+# Nearest without decoding the data values
+# ----------------------------------------
+${examples_dir}/c_grib_nearest -n $input_grb > $temp
+cat > $tempRef <<EOF
+ordering by param,step
+1 fields in the fieldset
+n,step,param
+1 0 t  
+Idx	lat	lon	dist
+4839	-40.46	18.00	259.679
+4838	-40.46	15.00	51.5268
+4719	-37.67	18.00	366.445
+4718	-37.67	15.00	258.597
+
+EOF
+diff $tempRef $temp
+
+# Clean up
+rm -f $temp $tempRef

src/grib_nearest.c

@@ -16,6 +16,9 @@
 
 #include "grib_api_internal.h"
 
+/* Note: The 'values' argument can be NULL in which case the data section will not be decoded
+ * See ECC-499
+ */
 int grib_nearest_find(
         grib_nearest *nearest, const grib_handle* ch,
         double inlat, double inlon,
@@ -147,17 +150,22 @@ int grib_nearest_find_multiple(
     double qoutlats[4]={0,};
     double qoutlons[4]={0,};
     double qvalues[4]={0,};
+    double* rvalues = NULL;
     int qindexes[4]={0,};
     int ret=0;
     long i=0;
     size_t len=4;
     int flags=GRIB_NEAREST_SAME_GRID | GRIB_NEAREST_SAME_DATA;
 
+    if (values) rvalues = qvalues;
+
     nearest=grib_nearest_new(h,&ret);
     if (ret!=GRIB_SUCCESS) return ret;
 
     if (is_lsm) {
         int noland=1;
+        /* ECC-499: In land-sea mask mode, 'values' cannot be NULL because we need to query whether >= 0.5 */
+        Assert(values);
         for (i=0;i<npoints;i++) {
             ret=grib_nearest_find(nearest,h,inlats[i],inlons[i],flags,qoutlats,qoutlons,
                     qvalues,qdistances,qindexes,&len);
@@ -181,9 +189,10 @@ int grib_nearest_find_multiple(
 
         }
     } else {
+        /* ECC-499: 'values' can be NULL */
         for (i=0;i<npoints;i++) {
             ret=grib_nearest_find(nearest,h,inlats[i],inlons[i],flags,qoutlats,qoutlons,
-                    qvalues,qdistances,qindexes,&len);
+                    rvalues,qdistances,qindexes,&len);
             min=qdistances[0];
             for (ii=0;ii<4;ii++) {
                 if ((min >= qdistances[ii])) {
@@ -193,7 +202,9 @@ int grib_nearest_find_multiple(
             }
             *poutlats=qoutlats[idx];poutlats++;
             *poutlons=qoutlons[idx];poutlons++;
+            if (values) {
                 *pvalues=qvalues[idx];pvalues++;
+            }
             *pdistances=qdistances[idx];pdistances++;
             *pindexes=qindexes[idx];pindexes++;
         }

src/grib_nearest_class_latlon_reduced.c

@@ -328,7 +328,9 @@ static int find(grib_nearest* nearest, grib_handle* h,
             distances[kk]=self->distances[kk];
             outlats[kk]=self->lats[self->j[jj]];
             outlons[kk]=self->lons[self->k[kk]];
+            if (values) { /* ECC-499 */
                 grib_get_double_element_internal(h,self->values_key,self->k[kk],&(values[kk]));
+            }
             indexes[kk]=self->k[kk];
             kk++;
         }

src/grib_nearest_class_reduced.c

@@ -382,7 +382,9 @@ static int find(grib_nearest* nearest, grib_handle* h,
             distances[kk]=self->distances[kk];
             outlats[kk]=self->lats[self->j[jj]];
             outlons[kk]=self->lons[self->k[kk]];
+            if (values) { /* ECC-499 */
                 grib_get_double_element_internal(h,self->values_key,self->k[kk],&(values[kk]));
+            }
             indexes[kk]=self->k[kk];
             kk++;
         }

src/grib_nearest_class_regular.c

@@ -420,7 +420,9 @@ static int find(grib_nearest* nearest, grib_handle* h,
                 outlats[kk] = new_lat;
                 outlons[kk] = new_lon;
             }
+            if (values) { /* ECC-499 */
                 grib_get_double_element_internal(h,self->values_key,self->k[kk],&(values[kk]));
+            }
             /* Using the brute force approach described above */
             /* Assert(self->k[kk] < nvalues); */
             /* values[kk]=nearest->values[self->k[kk]]; */