Added legacyGaussSubarea and use this in nearest

definitions/grib1/grid_definition_gaussian.def

@@ -52,10 +52,17 @@ meta global global_gaussian(N,Ni,iDirectionIncrement,
                             longitudeOfLastGridPoint,
                             PLPresent,pl) = 0 : dump;
 
-meta numberOfDataPoints number_of_points_gaussian(Ni,Nj,PLPresent,pl,
+# With legacy mode support
-  N,
+meta numberOfDataPoints number_of_points_gaussian(Ni,Nj,PLPresent,pl,N,
   latitudeOfFirstGridPointInDegrees,longitudeOfFirstGridPointInDegrees,
-  latitudeOfLastGridPointInDegrees,longitudeOfLastGridPointInDegrees) : dump;
+  latitudeOfLastGridPointInDegrees,longitudeOfLastGridPointInDegrees,one) : dump;
+
+# Use the new algorithm for counting. No support for legacy mode
+meta numberOfDataPointsExpected number_of_points_gaussian(Ni,Nj,PLPresent,pl,N,
+  latitudeOfFirstGridPointInDegrees,longitudeOfFirstGridPointInDegrees,
+  latitudeOfLastGridPointInDegrees,longitudeOfLastGridPointInDegrees,zero) : dump;
+
+meta legacyGaussSubarea evaluate(numberOfDataPoints != numberOfDataPointsExpected);
 
 alias numberOfPoints=numberOfDataPoints;
 # alias numberOfExpectedPoints=numberOfDataPoints;

definitions/grib2/template.3.gaussian.def

@@ -89,9 +89,10 @@ meta gaussianGridName gaussian_grid_name(N, Ni, isOctahedral);
 alias gridName=gaussianGridName;
 
 
-# Useful for sub-areas
+# For sub-areas
-# meta numberOfExpectedPoints number_of_points_gaussian(Ni,Nj,PLPresent,pl,
+# Uses new algorithm for counting. No support for legacy mode
-#    N,
+meta numberOfDataPointsExpected number_of_points_gaussian(Ni,Nj,PLPresent,pl,N,
-#    latitudeOfFirstGridPointInDegrees,longitudeOfFirstGridPointInDegrees,
+  latitudeOfFirstGridPointInDegrees,longitudeOfFirstGridPointInDegrees,
-#    latitudeOfLastGridPointInDegrees,longitudeOfLastGridPointInDegrees) : dump;
+  latitudeOfLastGridPointInDegrees,longitudeOfLastGridPointInDegrees,zero) : dump;
 
+meta legacyGaussSubarea evaluate(numberOfDataPoints != numberOfDataPointsExpected);

src/grib_accessor_class_number_of_points_gaussian.c

@@ -32,6 +32,7 @@
    MEMBERS = const char* lon_first
    MEMBERS = const char* lat_last
    MEMBERS = const char* lon_last
+   MEMBERS = const char* support_legacy
    END_CLASS_DEF
  */
 
@@ -63,6 +64,7 @@ typedef struct grib_accessor_number_of_points_gaussian {
 	const char* lon_first;
 	const char* lat_last;
 	const char* lon_last;
+	const char* support_legacy;
 } grib_accessor_number_of_points_gaussian;
 
 extern grib_accessor_class* grib_accessor_class_long;
@@ -167,6 +169,7 @@ static void init(grib_accessor* a,const long l, grib_arguments* c)
     self->lon_first = grib_arguments_get_name(h,c,n++);
     self->lat_last = grib_arguments_get_name(h,c,n++);
     self->lon_last = grib_arguments_get_name(h,c,n++);
+    self->support_legacy = grib_arguments_get_name(h,c,n++);
     a->flags  |= GRIB_ACCESSOR_FLAG_READ_ONLY;
     a->flags |= GRIB_ACCESSOR_FLAG_FUNCTION;
     a->length=0;
@@ -303,7 +306,133 @@ static int get_number_of_data_values(grib_handle* h, size_t* numDataValues)
     return err;
 }
 
+static int unpack_long_with_legacy_support(grib_accessor* a, long* val, size_t *len);
+static int unpack_long_new(grib_accessor* a, long* val, size_t *len);
+
 static int unpack_long(grib_accessor* a, long* val, size_t *len)
+{
+    int ret=GRIB_SUCCESS;
+    long support_legacy=1;
+    grib_accessor_number_of_points_gaussian* self = (grib_accessor_number_of_points_gaussian*)a;
+    grib_handle* h = grib_handle_of_accessor(a);
+
+    if((ret = grib_get_long_internal(h, self->support_legacy,&support_legacy)) != GRIB_SUCCESS)
+        return ret;
+    if (support_legacy==1) return unpack_long_with_legacy_support(a, val, len);
+    else                   return unpack_long_new(a, val, len);
+}
+
+// ------------
+// New way
+static int unpack_long_new(grib_accessor* a, long* val, size_t *len)
+{
+    int ret=GRIB_SUCCESS;
+    int is_global = 0;
+    long ni=0,nj=0,plpresent=0,order=0;
+    size_t plsize=0;
+    double* lats={0,};
+    double lat_first,lat_last,lon_first,lon_last;
+    long* pl=NULL;
+    long* plsave=NULL;
+    long row_count;
+    long ilon_first=0,ilon_last=0;
+    double angular_precision = 1.0/1000000.0;
+    long editionNumber = 0;
+    grib_handle* h = grib_handle_of_accessor(a);
+
+    grib_accessor_number_of_points_gaussian* self = (grib_accessor_number_of_points_gaussian*)a;
+    grib_context* c=a->context;
+
+    if((ret = grib_get_long_internal(h, self->ni,&ni)) != GRIB_SUCCESS)
+        return ret;
+
+    if((ret = grib_get_long_internal(h, self->nj,&nj)) != GRIB_SUCCESS)
+        return ret;
+
+    if((ret = grib_get_long_internal(h, self->plpresent,&plpresent)) != GRIB_SUCCESS)
+        return ret;
+
+    if (nj == 0) return GRIB_GEOCALCULUS_PROBLEM;
+
+    if (grib_get_long(h, "editionNumber", &editionNumber)==GRIB_SUCCESS) {
+        if (editionNumber == 1) angular_precision = 1.0/1000;
+    }
+
+    if (plpresent) {
+        long max_pl=0;
+        float d = 0;
+        int j=0;
+        double lon_first_row=0,lon_last_row=0;
+
+        /*reduced*/
+        if((ret = grib_get_long_internal(h, self->order,&order)) != GRIB_SUCCESS)
+            return ret;
+        if((ret = grib_get_double_internal(h, self->lat_first,&lat_first)) != GRIB_SUCCESS)
+            return ret;
+        if((ret = grib_get_double_internal(h, self->lon_first,&lon_first)) != GRIB_SUCCESS)
+            return ret;
+        if((ret = grib_get_double_internal(h, self->lat_last,&lat_last)) != GRIB_SUCCESS)
+            return ret;
+        if((ret = grib_get_double_internal(h, self->lon_last,&lon_last)) != GRIB_SUCCESS)
+            return ret;
+
+        lats=(double*)grib_context_malloc(a->context,sizeof(double)*order*2);
+        if((ret = grib_get_gaussian_latitudes(order, lats)) != GRIB_SUCCESS)
+            return ret;
+
+        if((ret = grib_get_size(h,self->pl,&plsize)) != GRIB_SUCCESS)
+            return ret;
+
+        pl=(long*)grib_context_malloc_clear(c,sizeof(long)*plsize);
+        plsave=pl;
+        grib_get_long_array_internal(h,self->pl,pl, &plsize);
+
+        if (lon_last<0) lon_last+=360;
+        if (lon_first<0) lon_first+=360;
+
+        /* Find the maximum element of "pl" array, do not assume it's 4*N! */
+        /* This could be an Octahedral Gaussian Grid */
+        max_pl = pl[0];
+        for (j=1; j<plsize; j++) {
+            if (pl[j] > max_pl) max_pl = pl[j];
+        }
+
+        d=fabs(lats[0]-lats[1]);
+        is_global = 0; /* ECC-445 */
+
+        correctWestEast(max_pl, angular_precision, &lon_first, &lon_last);
+
+        if ( !is_global ) {
+            /*sub area*/
+            (void)d;
+            *val=0;
+            for (j=0;j<nj;j++) {
+                row_count=0;
+                grib_get_reduced_row_wrapper(h, pl[j],lon_first,lon_last,&row_count,&ilon_first,&ilon_last);
+                lon_first_row=((ilon_first)*360.0)/pl[j];
+                lon_last_row=((ilon_last)*360.0)/pl[j];
+                *val+=row_count;
+                (void)lon_last_row;
+                (void)lon_first_row;
+            }
+        } else {
+            int i = 0;
+            *val=0;
+            for (i=0;i<plsize;i++) *val+=pl[i];
+        }
+    } else {
+        /*regular*/
+        *val=ni*nj;
+    }
+    if (lats) grib_context_free(c,lats);
+    if (plsave) grib_context_free(c,plsave);
+
+    return ret;
+}
+
+
+// With Legacy support
+static int unpack_long_with_legacy_support(grib_accessor* a, long* val, size_t *len)
 {
     int ret=GRIB_SUCCESS;
     int is_global = 0;

src/grib_gaussian_reduced.c

@@ -310,7 +310,8 @@ void grib_get_reduced_row_wrapper(grib_handle* h, long pl, double lon_first, dou
      */
 }
 
-/* This was the legacy way of counting the points. Now deprecated */
+/* This was the legacy way of counting the points within a subarea of a Gaussian grid.
+   In the days of Prodgen/libemos */
 void grib_get_reduced_row_legacy(long pl, double lon_first, double lon_last, long* npoints, long* ilon_first, long* ilon_last)
 {
     double range=0,dlon_first=0,dlon_last=0;

src/grib_nearest_class_reduced.c

@@ -122,6 +122,14 @@ static int init(grib_nearest* nearest,grib_handle* h,grib_arguments* args)
     return 0;
 }
 
+typedef void (*get_reduced_row_proc)(long pl, double lon_first, double lon_last, long* npoints, long* ilon_first, long* ilon_last);
+
+static int is_legacy(grib_handle* h)
+{
+    long is_legacy = 0;
+    return (grib_get_long(h, "legacyGaussSubarea", &is_legacy)==GRIB_SUCCESS && is_legacy==1);
+}
+
 #if 1
 static int find(grib_nearest* nearest, grib_handle* h,
         double inlat, double inlon,unsigned long flags,
@@ -137,6 +145,11 @@ static int find(grib_nearest* nearest, grib_handle* h,
     long iradius;
     double radius;
     int ilat=0,ilon=0;
+    const int is_legacy_grib = is_legacy(h);
+    get_reduced_row_proc get_reduced_row_func = &grib_get_reduced_row;
+    if (is_legacy_grib) {
+        get_reduced_row_func = &grib_get_reduced_row_legacy;
+    }
 
     if( (ret =  grib_get_size(h,self->values_key,&nvalues))!= GRIB_SUCCESS)
         return ret;
@@ -188,7 +201,8 @@ static int find(grib_nearest* nearest, grib_handle* h,
             }
             while(lon>360) lon-=360;
             if(!self->global) {  /* ECC-756 */
-                if (lon>180 && lon<360) lon-=360;
+                if (!is_legacy_grib) //TODO
+                    if (lon>180 && lon<360) lon-=360;
             }
             self->lons[ilon++]=lon;
         }
@@ -247,11 +261,11 @@ static int find(grib_nearest* nearest, grib_handle* h,
             nlon=0;
             for (jj=0;jj<self->j[0];jj++) {
                 row_count=0;ilon_first=0;ilon_last=0;
-                grib_get_reduced_row_wrapper(h, pl[jj],self->lon_first,self->lon_last,&row_count,&ilon_first,&ilon_last);
+                get_reduced_row_func(pl[jj],self->lon_first,self->lon_last,&row_count,&ilon_first,&ilon_last);
                 nlon+=row_count;
             }
             row_count=0;ilon_first=0;ilon_last=0;
-            grib_get_reduced_row_wrapper(h, pl[self->j[0]],self->lon_first,self->lon_last,&row_count,&ilon_first,&ilon_last);
+            get_reduced_row_func(pl[self->j[0]],self->lon_first,self->lon_last,&row_count,&ilon_first,&ilon_last);
             nplm1=row_count-1;
         }
         lons=self->lons+nlon;
@@ -283,7 +297,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
         if (!nearest_lons_found) {
             if (!self->global) {
                 row_count=0;ilon_first=0;ilon_last=0;
-                grib_get_reduced_row_wrapper(h, pl[self->j[0]],self->lon_first,self->lon_last,&row_count,&ilon_first,&ilon_last);
+                get_reduced_row_func(pl[self->j[0]],self->lon_first,self->lon_last,&row_count,&ilon_first,&ilon_last);
             } else {
                 row_count=pl[self->j[0]];
             }
@@ -301,11 +315,11 @@ static int find(grib_nearest* nearest, grib_handle* h,
         } else {
             for (jj=0;jj<self->j[1];jj++) {
                 row_count=0;ilon_first=0;ilon_last=0;
-                grib_get_reduced_row_wrapper(h, pl[jj],self->lon_first,self->lon_last,&row_count,&ilon_first,&ilon_last);
+                get_reduced_row_func(pl[jj],self->lon_first,self->lon_last,&row_count,&ilon_first,&ilon_last);
                 nlon+=row_count;
             }
             row_count=0;ilon_first=0;ilon_last=0;
-            grib_get_reduced_row_wrapper(h, pl[self->j[1]],self->lon_first,self->lon_last,&nplm1,&ilon_first,&ilon_last);
+            get_reduced_row_func(pl[self->j[1]],self->lon_first,self->lon_last,&nplm1,&ilon_first,&ilon_last);
             nplm1--;
         }
         lons=self->lons+nlon;
@@ -334,7 +348,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
         if (!nearest_lons_found) {
             if (!self->global) {
                 row_count=0;ilon_first=0;ilon_last=0;
-                grib_get_reduced_row_wrapper(h, pl[self->j[1]],self->lon_first,self->lon_last,&row_count,&ilon_first,&ilon_last);
+                get_reduced_row_func(pl[self->j[1]],self->lon_first,self->lon_last,&row_count,&ilon_first,&ilon_last);
             } else {
                 row_count=pl[self->j[1]];
             }