core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

ECC-906: grib_get_data not working correctly with old-style sub-areas of reduced grids produced by PRODGEN

This commit is contained in:
Shahram Najm 2019-04-15 12:05:31 +01:00
parent 18595a0d10 2dcb46ed5a
commit 298fb4e4ac
7 changed files with 235 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
definitions/grib1/grid_definition_gaussian.def
View File

@ -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,
N,
# With legacy mode support
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;

11
definitions/grib2/template.3.gaussian.def
View File

@ -89,9 +89,10 @@ meta gaussianGridName gaussian_grid_name(N, Ni, isOctahedral);
alias gridName=gaussianGridName;
# Useful for sub-areas
# meta numberOfExpectedPoints number_of_points_gaussian(Ni,Nj,PLPresent,pl,
# N,
# latitudeOfFirstGridPointInDegrees,longitudeOfFirstGridPointInDegrees,
# latitudeOfLastGridPointInDegrees,longitudeOfLastGridPointInDegrees) : dump;
# For sub-areas
# Uses 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);

128
src/grib_accessor_class_number_of_points_gaussian.c
View File

@ -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,132 @@ 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 algorithm */
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;

1
src/grib_api_prototypes.h
View File

@ -809,6 +809,7 @@ long accessor_raw_get_offset(grib_accessor *a);
/* grib_gaussian_reduced.c */
void grib_get_reduced_row_wrapper(grib_handle *h, long pl, double lon_first, double lon_last, long *npoints, long *ilon_first, long *ilon_last);
void grib_get_reduced_row(long pl, double lon_first, double lon_last, long *npoints, long *ilon_first, long *ilon_last);
void grib_get_reduced_row_legacy(long pl, double lon_first, double lon_last, long *npoints, long *ilon_first, long *ilon_last);
void grib_get_reduced_row_p(long pl, double lon_first, double lon_last, long *npoints, double *olon_first, double *olon_last);
/* grib_accessor_class_abstract_vector.c */

8
src/grib_gaussian_reduced.c
View File

@ -310,9 +310,9 @@ void grib_get_reduced_row_wrapper(grib_handle* h, long pl, double lon_first, dou
*/
}
#if 0
/* This was the legacy way of counting the points. Now deprecated */
static void grib_get_reduced_row1(long pl, double lon_first, double lon_last, long* npoints, long* ilon_first, long* ilon_last )
/* 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;
long irange;
@ -393,14 +393,12 @@ static void grib_get_reduced_row1(long pl, double lon_first, double lon_last, lo
pl,*npoints,range,*ilon_first,*ilon_last,irange);
#endif
}
}
if (*ilon_first<0) *ilon_first+=pl;
return;
}
#endif
/* New method based on eckit Fractions and matching MIR count */
void grib_get_reduced_row(long pl, double lon_first, double lon_last, long* npoints, long* ilon_first, long* ilon_last )

120
src/grib_iterator_class_gaussian_reduced.c
View File

@ -149,6 +149,60 @@ static void binary_search(double xx[], const unsigned long n, double x, long *j)
*j=jl;
}
/* Use legacy way to compute the iterator latitude/longitude values */
static int iterate_reduced_gaussian_subarea_legacy(grib_iterator* iter, grib_handle* h,
double lat_first, double lon_first,
double lat_last, double lon_last,
double* lats, long* pl, size_t plsize)
{
int err = 0;
int l = 0;
size_t j = 0;
long row_count=0;
double d=0;
long ilon_first, ilon_last, i;
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter;
get_reduced_row_proc get_reduced_row = &grib_get_reduced_row;
get_reduced_row = &grib_get_reduced_row_legacy; /* legacy algorithm */
if (h->context->debug) {
const size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last);
printf("ECCODES DEBUG grib_iterator_class_gaussian_reduced: Legacy sub-area num points=%ld\n", (long)np);
}
/*find starting latitude */
d = fabs(lats[0] - lats[1]);
while (fabs(lat_first-lats[l]) > d ) {l++;}
iter->e=0;
for (j=0;j<plsize;j++) {
long k = 0;
row_count=0;
get_reduced_row(pl[j],lon_first,lon_last, &row_count,&ilon_first,&ilon_last);
/*printf("iterate_reduced_gaussian_subarea %ld %g %g count=%ld, (i1=%ld, i2=%ld)\n",pl[j],lon_first,lon_last, row_count, ilon_first,ilon_last);*/
if (ilon_first>ilon_last) ilon_first-=pl[j];
for (i=ilon_first;i<=ilon_last;i++) {
if(iter->e >= iter->nv){
size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last);
grib_context_log(h->context,GRIB_LOG_ERROR,
"Reduced Gaussian iterator (sub-area legacy). Num points=%ld, size(values)=%ld", np, iter->nv);
return GRIB_WRONG_GRID;
}
self->los[iter->e]=((i)*360.0)/pl[j];
self->las[iter->e]=lats[j+l];
iter->e++;
k++;
if (k >= row_count) {
/* Ensure we exit the loop and only process 'row_count' points */
break;
}
}
}
return err;
}
/* ECC-747 */
static int iterate_reduced_gaussian_subarea_algorithm2(grib_iterator* iter, grib_handle* h,
double lat_first, double lon_first,
@ -202,69 +256,19 @@ static int iterate_reduced_gaussian_subarea_algorithm2(grib_iterator* iter, grib
iter->e++;
}
}
return err;
}
#if 0
/* ECC-445: Try to compute the iterator latitude/longitude values. If algorithm2 is set, try a different point count */
static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h,
double lat_first, double lon_first,
double lat_last, double lon_last,
double* lats, long* pl, size_t plsize, int algorithm2)
{
int err = 0;
int l = 0;
size_t j = 0;
long row_count=0;
double d=0;
long ilon_first, ilon_last, i;
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter;
get_reduced_row_proc get_reduced_row = &grib_get_reduced_row;
if (algorithm2) {
get_reduced_row = &grib_get_reduced_row2; /* switch to 2nd algorithm */
}
if (h->context->debug) {
const size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last);
printf("ECCODES DEBUG grib_iterator_class_gaussian_reduced: sub-area num points=%ld\n", (long)np);
}
/*find starting latitude */
d = fabs(lats[0] - lats[1]);
while (fabs(lat_first-lats[l]) > d ) {l++;}
iter->e=0;
for (j=0;j<plsize;j++) {
long k = 0;
row_count=0;
get_reduced_row(pl[j],lon_first,lon_last, &row_count,&ilon_first,&ilon_last);
/*printf("iterate_reduced_gaussian_subarea %ld %g %g count=%ld, (i1=%ld, i2=%ld)\n",pl[j],lon_first,lon_last, row_count, ilon_first,ilon_last);*/
if (ilon_first>ilon_last) ilon_first-=pl[j];
for (i=ilon_first;i<=ilon_last;i++) {
if(iter->e >= iter->nv){
if (algorithm2) {
/* Only print error message on the second pass */
size_t np = count_subarea_points(h, get_reduced_row, pl, plsize, lon_first, lon_last);
grib_context_log(h->context,GRIB_LOG_ERROR,
"Reduced Gaussian iterator (sub-area). Num points=%ld, size(values)=%ld", np, iter->nv);
}
return GRIB_WRONG_GRID;
}
self->los[iter->e]=((i)*360.0)/pl[j];
self->las[iter->e]=lats[j+l];
iter->e++;
k++;
if (k >= row_count) {
/* Ensure we exit the loop and only process 'row_count' points */
break;
}
if (iter->e != iter->nv) {
/* Fewer counted points in the sub-area than the number of data values */
const size_t legacy_count = count_subarea_points(h, grib_get_reduced_row_legacy, pl, plsize, lon_first, lon_last);
if (iter->nv == legacy_count) {
/* Legacy (produced by PRODGEN/LIBEMOS) */
return iterate_reduced_gaussian_subarea_legacy(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize);
} else {
/* TODO: A gap exists! Not all values can be mapped. Inconsistent grid or error in calculating num. points! */
}
}
return err;
}
#endif
static int iterate_reduced_gaussian_subarea_wrapper(grib_iterator* iter, grib_handle* h,
double lat_first, double lon_first,

30
src/grib_nearest_class_reduced.c
View File

@ -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,6 +201,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
}
while(lon>360) lon-=360;
if(!self->global) { /* ECC-756 */
if (!is_legacy_grib) /*TODO*/
if (lon>180 && lon<360) lon-=360;
}
self->lons[ilon++]=lon;
@ -210,6 +224,10 @@ static int find(grib_nearest* nearest, grib_handle* h,
if (self->global) {
while (inlon<0) inlon+=360;
while (inlon>360) inlon-=360;
} else {
/* TODO: Experimental */
if (!is_legacy_grib)
if (inlon>180 && inlon<360) inlon-=360;
}
ilat=self->lats_count;
@ -247,11 +265,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 +301,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 +319,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 +352,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]];
}