mirror of https://github.com/ecmwf/eccodes.git
ECC-904, ECC-905 and ECC-906
This commit is contained in:
parent
54dc9b028c
commit
0e67d66921
|
@ -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;
|
||||
|
@ -96,4 +103,3 @@ meta isOctahedral octahedral_gaussian(N, Ni, PLPresent, pl) = 0 : no_copy,dump;
|
|||
|
||||
meta gaussianGridName gaussian_grid_name(N, Ni, isOctahedral);
|
||||
alias gridName=gaussianGridName;
|
||||
|
||||
|
|
|
@ -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);
|
||||
|
|
|
@ -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;
|
||||
|
@ -437,7 +565,8 @@ static int unpack_long(grib_accessor* a, long* val, size_t *len)
|
|||
if (get_number_of_data_values(h, &numDataValues) == GRIB_SUCCESS) {
|
||||
if (*val != numDataValues) {
|
||||
if (h->context->debug)
|
||||
printf("ECCODES DEBUG number_of_points_gaussian: LEGACY MODE activated. Count(=%ld) changed to size(values)\n",*val);
|
||||
printf("ECCODES DEBUG number_of_points_gaussian: LEGACY MODE activated. "
|
||||
"Count(=%ld) changed to num values(=%ld)\n",*val,(long)numDataValues);
|
||||
*val = numDataValues;
|
||||
}
|
||||
}
|
||||
|
|
|
@ -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 */
|
||||
|
|
|
@ -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 )
|
||||
|
|
|
@ -14,6 +14,7 @@
|
|||
|
||||
|
||||
#include "grib_api_internal.h"
|
||||
#include <float.h>
|
||||
#include <math.h>
|
||||
|
||||
/*
|
||||
|
@ -126,58 +127,33 @@ static size_t count_subarea_points(grib_handle* h, get_reduced_row_proc get_redu
|
|||
return result;
|
||||
}
|
||||
|
||||
/* ECC-747 */
|
||||
static int iterate_reduced_gaussian_subarea_algorithm2(grib_iterator* iter, grib_handle* h,
|
||||
/* Search for 'x' in the array 'xx' (the index of last element being 'n') and return index in 'j' */
|
||||
static void binary_search(double xx[], const unsigned long n, double x, long *j)
|
||||
{
|
||||
/*This routine works only on descending ordered arrays*/
|
||||
#define EPSILON 1e-3
|
||||
|
||||
unsigned long ju,jm,jl;
|
||||
jl=0;
|
||||
ju=n;
|
||||
while (ju-jl > 1) {
|
||||
jm=(ju+jl) >> 1;
|
||||
if (fabs(x-xx[jm]) < EPSILON) {
|
||||
/* found something close enough. We're done */
|
||||
*j=jm;
|
||||
return;
|
||||
}
|
||||
if (x < xx[jm]) jl=jm;
|
||||
else ju=jm;
|
||||
}
|
||||
*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, i=0;
|
||||
double d=0;
|
||||
double olon_first, olon_last;
|
||||
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter;
|
||||
get_reduced_row_proc get_reduced_row = &grib_get_reduced_row;
|
||||
|
||||
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++) {
|
||||
const double delta = 360.0/pl[j];
|
||||
row_count=0;
|
||||
grib_get_reduced_row_p(pl[j],lon_first,lon_last, &row_count,&olon_first,&olon_last);
|
||||
for(i=0; i<row_count; ++i) {
|
||||
double lon2 = olon_first + i * delta;
|
||||
if(iter->e >= iter->nv){
|
||||
/* 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]=lon2;
|
||||
self->las[iter->e]=lats[j+l];
|
||||
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;
|
||||
|
@ -187,13 +163,11 @@ static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h,
|
|||
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 */
|
||||
}
|
||||
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: sub-area num points=%ld\n", (long)np);
|
||||
printf("ECCODES DEBUG grib_iterator_class_gaussian_reduced: Legacy sub-area num points=%ld\n", (long)np);
|
||||
}
|
||||
|
||||
/*find starting latitude */
|
||||
|
@ -210,12 +184,9 @@ static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h,
|
|||
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);
|
||||
}
|
||||
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;
|
||||
}
|
||||
|
||||
|
@ -231,14 +202,80 @@ static int iterate_reduced_gaussian_subarea(grib_iterator* iter, grib_handle* h,
|
|||
}
|
||||
return err;
|
||||
}
|
||||
|
||||
/* ECC-747 */
|
||||
static int iterate_reduced_gaussian_subarea_algorithm2(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, size_t numlats)
|
||||
{
|
||||
int err = 0;
|
||||
long l = 0;
|
||||
size_t j = 0;
|
||||
long row_count=0, i=0;
|
||||
double olon_first, olon_last;
|
||||
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)iter;
|
||||
get_reduced_row_proc get_reduced_row = &grib_get_reduced_row;
|
||||
|
||||
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 */
|
||||
binary_search(lats, numlats-1, lat_first, &l);
|
||||
Assert(l < numlats);
|
||||
|
||||
#if 0
|
||||
for(il=0; il<numlats; ++il) {
|
||||
const double diff = fabs(lat_first-lats[il]);
|
||||
if (diff < min_d) {
|
||||
min_d = diff;
|
||||
l = il; /* index of the latitude */
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
iter->e=0;
|
||||
for (j=0;j<plsize;j++) {
|
||||
const double delta = 360.0/pl[j];
|
||||
row_count=0;
|
||||
grib_get_reduced_row_p(pl[j],lon_first,lon_last, &row_count,&olon_first,&olon_last);
|
||||
for(i=0; i<row_count; ++i) {
|
||||
double lon2 = olon_first + i * delta;
|
||||
if(iter->e >= iter->nv){
|
||||
/* 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]=lon2;
|
||||
DebugAssert(j+l < numlats);
|
||||
self->las[iter->e]=lats[j+l];
|
||||
iter->e++;
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
static int iterate_reduced_gaussian_subarea_wrapper(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)
|
||||
double* lats, long* pl, size_t plsize, size_t numlats)
|
||||
{
|
||||
return iterate_reduced_gaussian_subarea_algorithm2(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize);
|
||||
return iterate_reduced_gaussian_subarea_algorithm2(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize, numlats);
|
||||
|
||||
#if 0
|
||||
/* Try legacy approach, if that fails try the next algorithm */
|
||||
|
@ -258,6 +295,7 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
|
|||
double angular_precision = 1.0/1000000.0;
|
||||
double* lats;
|
||||
size_t plsize=0;
|
||||
size_t numlats=0;
|
||||
long* pl;
|
||||
long max_pl=0;
|
||||
long nj=0,order=0,i;
|
||||
|
@ -291,7 +329,8 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
|
|||
if (editionNumber == 1) angular_precision = 1.0/1000;
|
||||
}
|
||||
|
||||
lats=(double*)grib_context_malloc(h->context,sizeof(double)*order*2);
|
||||
numlats = order*2;
|
||||
lats=(double*)grib_context_malloc(h->context,sizeof(double)*numlats);
|
||||
if (!lats) return GRIB_OUT_OF_MEMORY;
|
||||
if((ret = grib_get_gaussian_latitudes(order, lats)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
@ -323,7 +362,7 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
|
|||
is_global = is_gaussian_global(lat_first, lat_last, lon_first, lon_last, max_pl, lats, angular_precision);
|
||||
if ( !is_global ) {
|
||||
/*sub area*/
|
||||
ret = iterate_reduced_gaussian_subarea_wrapper(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize);
|
||||
ret = iterate_reduced_gaussian_subarea_wrapper(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize, numlats);
|
||||
} else {
|
||||
/*global*/
|
||||
iter->e=0;
|
||||
|
@ -340,7 +379,7 @@ static int init(grib_iterator* iter, grib_handle* h, grib_arguments* args)
|
|||
/*grib_context_log(h->context,GRIB_LOG_ERROR, "Failed to initialise reduced Gaussian iterator (global)");*/
|
||||
/*return GRIB_WRONG_GRID;*/
|
||||
/*Try now as NON-global*/
|
||||
ret = iterate_reduced_gaussian_subarea_wrapper(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize);
|
||||
ret = iterate_reduced_gaussian_subarea_wrapper(iter, h, lat_first, lon_first, lat_last, lon_last, lats, pl, plsize, numlats);
|
||||
if (ret !=GRIB_SUCCESS) grib_context_log(h->context,GRIB_LOG_ERROR, "Failed to initialise reduced Gaussian iterator (global)");
|
||||
goto finalise;
|
||||
}
|
||||
|
|
|
@ -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;
|
||||
}
|
||||
|
@ -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]];
|
||||
}
|
||||
|
|
Loading…
Reference in New Issue