mirror of https://github.com/ecmwf/eccodes.git
GRIB-704: Support Octahedral Gaussian grids
This commit is contained in:
Shahram Najm
parent
c2f8d5b42f
commit
b33696cf6d
|
|
@ -153,124 +153,132 @@ static void init_class(grib_accessor_class* c)
|
|||
|
||||
static void init(grib_accessor* a,const long l, grib_arguments* c)
|
||||
{
|
||||
int n=0;
|
||||
grib_accessor_number_of_points_gaussian* self = (grib_accessor_number_of_points_gaussian*)a;
|
||||
self->ni = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->nj = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->plpresent = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->pl = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->order = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->lat_first = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->lon_first = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->lat_last = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->lon_last = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
a->flags |= GRIB_ACCESSOR_FLAG_READ_ONLY;
|
||||
a->flags |= GRIB_ACCESSOR_FLAG_FUNCTION;
|
||||
a->length=0;
|
||||
int n=0;
|
||||
grib_accessor_number_of_points_gaussian* self = (grib_accessor_number_of_points_gaussian*)a;
|
||||
self->ni = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->nj = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->plpresent = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->pl = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->order = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->lat_first = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->lon_first = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->lat_last = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
self->lon_last = grib_arguments_get_name(a->parent->h,c,n++);
|
||||
a->flags |= GRIB_ACCESSOR_FLAG_READ_ONLY;
|
||||
a->flags |= GRIB_ACCESSOR_FLAG_FUNCTION;
|
||||
a->length=0;
|
||||
}
|
||||
|
||||
static int unpack_long(grib_accessor* a, long* val, size_t *len){
|
||||
int ret=GRIB_SUCCESS;
|
||||
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;
|
||||
float d;
|
||||
long* pl=NULL;
|
||||
long* plsave=NULL;
|
||||
int j=0,i=0;
|
||||
long row_count;
|
||||
long ilon_first=0,ilon_last=0;
|
||||
double lon_first_row=0,lon_last_row=0;
|
||||
static int unpack_long(grib_accessor* a, long* val, size_t *len)
|
||||
{
|
||||
int ret=GRIB_SUCCESS;
|
||||
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;
|
||||
|
||||
grib_accessor_number_of_points_gaussian* self = (grib_accessor_number_of_points_gaussian*)a;
|
||||
grib_context* c=a->parent->h->context;
|
||||
|
||||
grib_accessor_number_of_points_gaussian* self = (grib_accessor_number_of_points_gaussian*)a;
|
||||
grib_context* c=a->parent->h->context;
|
||||
if((ret = grib_get_long_internal(a->parent->h, self->ni,&ni)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
if((ret = grib_get_long_internal(a->parent->h, self->ni,&ni)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_long_internal(a->parent->h, self->nj,&nj)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
if((ret = grib_get_long_internal(a->parent->h, self->nj,&nj)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_long_internal(a->parent->h, self->plpresent,&plpresent)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
if((ret = grib_get_long_internal(a->parent->h, self->plpresent,&plpresent)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if (nj == 0) return GRIB_GEOCALCULUS_PROBLEM;
|
||||
|
||||
if (nj == 0) return GRIB_GEOCALCULUS_PROBLEM;
|
||||
if (plpresent) {
|
||||
long max_pl=0;
|
||||
float d = 0;
|
||||
int j=0;
|
||||
double lon_first_row=0,lon_last_row=0;
|
||||
|
||||
if (plpresent) {
|
||||
/*reduced*/
|
||||
if((ret = grib_get_long_internal(a->parent->h, self->order,&order)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(a->parent->h, self->lat_first,&lat_first)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(a->parent->h, self->lon_first,&lon_first)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(a->parent->h, self->lat_last,&lat_last)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(a->parent->h, self->lon_last,&lon_last)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
/*reduced*/
|
||||
if((ret = grib_get_long_internal(a->parent->h, self->order,&order)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(a->parent->h, self->lat_first,&lat_first)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(a->parent->h, self->lon_first,&lon_first)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(a->parent->h, self->lat_last,&lat_last)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_double_internal(a->parent->h, self->lon_last,&lon_last)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
lats=(double*)grib_context_malloc(a->parent->h->context,sizeof(double)*order*2);
|
||||
if((ret = grib_get_gaussian_latitudes(order, lats)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
lats=(double*)grib_context_malloc(a->parent->h->context,sizeof(double)*order*2);
|
||||
if((ret = grib_get_gaussian_latitudes(order, lats)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
if((ret = grib_get_size(a->parent->h,self->pl,&plsize)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
if((ret = grib_get_size(a->parent->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(a->parent->h,self->pl,pl, &plsize);
|
||||
pl=(long*)grib_context_malloc_clear(c,sizeof(long)*plsize);
|
||||
plsave=pl;
|
||||
grib_get_long_array_internal(a->parent->h,self->pl,pl, &plsize);
|
||||
|
||||
if (lon_last<0) lon_last+=360;
|
||||
if (lon_first<0) lon_first+=360;
|
||||
if (lon_last<0) lon_last+=360;
|
||||
if (lon_first<0) lon_first+=360;
|
||||
|
||||
d=fabs(lats[0]-lats[1]);
|
||||
if ( (fabs(lat_first-lats[0]) >= d ) ||
|
||||
(fabs(lat_last+lats[0]) >= d ) ||
|
||||
lon_first != 0 ||
|
||||
fabs(lon_last - (360.0-90.0/order)) > 90.0/order
|
||||
) {
|
||||
/*sub area*/
|
||||
/* 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]);
|
||||
if ( (fabs(lat_first-lats[0]) >= d ) ||
|
||||
(fabs(lat_last+lats[0]) >= d ) ||
|
||||
lon_first != 0 ||
|
||||
fabs(lon_last - (360.0-360.0/max_pl)) > 360.0/max_pl
|
||||
) {
|
||||
/*sub area*/
|
||||
#if EFDEBUG
|
||||
printf("-------- subarea fabs(lat_first-lats[0])=%g d=%g\n",fabs(lat_first-lats[0]),d);
|
||||
printf("-------- subarea fabs(lat_last+lats[0])=%g d=%g\n",fabs(lat_last+lats[0]),d);
|
||||
printf("-------- subarea lon_last=%g order=%ld 360.0-90.0/order=%g\n",
|
||||
lon_last,order,360.0-90.0/order);
|
||||
printf("-------- subarea lon_first=%g fabs(lon_last -( 360.0-90.0/order))=%g 90.0/order=%g\n",
|
||||
lon_first,fabs(lon_last - (360.0-90.0/order)),90.0/order);
|
||||
printf("-------- subarea fabs(lat_first-lats[0])=%g d=%g\n",fabs(lat_first-lats[0]),d);
|
||||
printf("-------- subarea fabs(lat_last+lats[0])=%g d=%g\n",fabs(lat_last+lats[0]),d);
|
||||
printf("-------- subarea lon_last=%g order=%ld 360.0-90.0/order=%g\n",
|
||||
lon_last,order,360.0-90.0/order);
|
||||
printf("-------- subarea lon_first=%g fabs(lon_last -( 360.0-90.0/order))=%g 90.0/order=%g\n",
|
||||
lon_first,fabs(lon_last - (360.0-90.0/order)),90.0/order);
|
||||
#endif
|
||||
*val=0;
|
||||
for (j=0;j<nj;j++) {
|
||||
row_count=0;
|
||||
*val=0;
|
||||
for (j=0;j<nj;j++) {
|
||||
row_count=0;
|
||||
#if EFDEBUG
|
||||
printf("-- %d ",j);
|
||||
printf("-- %d ",j);
|
||||
#endif
|
||||
grib_get_reduced_row(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;
|
||||
grib_get_reduced_row(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;
|
||||
#if EFDEBUG
|
||||
printf(" ilon_first=%ld lon_first=%.10e ilon_last=%ld lon_last=%.10e count=%ld row_count=%ld\n",
|
||||
ilon_first,lon_first_row,ilon_last,lon_last_row,*val,row_count);
|
||||
printf(" ilon_first=%ld lon_first=%.10e ilon_last=%ld lon_last=%.10e count=%ld row_count=%ld\n",
|
||||
ilon_first,lon_first_row,ilon_last,lon_last_row,*val,row_count);
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
} else {
|
||||
int i = 0;
|
||||
*val=0;
|
||||
for (i=0;i<plsize;i++) *val+=pl[i];
|
||||
}
|
||||
} else {
|
||||
*val=0;
|
||||
for (i=0;i<plsize;i++) *val+=pl[i];
|
||||
/*regular*/
|
||||
*val=ni*nj;
|
||||
}
|
||||
} else {
|
||||
/*regular*/
|
||||
*val=ni*nj;
|
||||
}
|
||||
#if EFDEBUG
|
||||
printf("DEBUG: number_of_points_gaussian=%ld plpresent=%ld plsize=%d\n",*val,plpresent,plsize);
|
||||
for (i=0;i<plsize;i++) printf(" DEBUG: pl[%d]=%ld\n",i,pl[i]);
|
||||
printf("DEBUG: number_of_points_gaussian=%ld plpresent=%ld plsize=%d\n",*val,plpresent,plsize);
|
||||
for (i=0;i<plsize;i++) printf(" DEBUG: pl[%d]=%ld\n",i,pl[i]);
|
||||
#endif
|
||||
if (lats) grib_context_free(c,lats);
|
||||
if (plsave) grib_context_free(c,plsave);
|
||||
return ret;
|
||||
if (lats) grib_context_free(c,lats);
|
||||
if (plsave) grib_context_free(c,plsave);
|
||||
return ret;
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -89,28 +89,27 @@ static void init_class(grib_iterator_class* c)
|
|||
|
||||
static int next(grib_iterator* i, double *lat, double *lon, double *val)
|
||||
{
|
||||
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)i;
|
||||
|
||||
grib_iterator_gaussian_reduced* self = (grib_iterator_gaussian_reduced*)i;
|
||||
if((long)i->e >= (long)(i->nv-1))
|
||||
return 0;
|
||||
i->e++;
|
||||
|
||||
if((long)i->e >= (long)(i->nv-1))
|
||||
return 0;
|
||||
i->e++;
|
||||
*lat = self->las[i->e];
|
||||
*lon = self->los[i->e];
|
||||
*val = i->data[i->e];
|
||||
|
||||
*lat = self->las[i->e];
|
||||
*lon = self->los[i->e];
|
||||
*val = i->data[i->e];
|
||||
|
||||
return 1;
|
||||
return 1;
|
||||
}
|
||||
|
||||
static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
|
||||
{
|
||||
int ret=GRIB_SUCCESS,j;
|
||||
int ret=GRIB_SUCCESS, j;
|
||||
double lat_first=0,lon_first=0,lat_last=0,lon_last=0,d=0;
|
||||
double* lats;
|
||||
size_t plsize=0;
|
||||
int l=0;
|
||||
long* pl;
|
||||
long max_pl=0;
|
||||
long nj=0,order=0,ilon_first,ilon_last,i;
|
||||
long row_count=0;
|
||||
grib_context* c=h->context;
|
||||
|
|
@ -145,6 +144,7 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
|
|||
if((ret = grib_get_size(h,spl,&plsize)) != GRIB_SUCCESS)
|
||||
return ret;
|
||||
|
||||
Assert(plsize);
|
||||
pl=(long*)grib_context_malloc(c,sizeof(long)*plsize);
|
||||
if (!pl) return GRIB_OUT_OF_MEMORY;
|
||||
|
||||
|
|
@ -158,12 +158,20 @@ static int init(grib_iterator* iter,grib_handle* h,grib_arguments* args)
|
|||
while (lon_last<0) lon_last+=360;
|
||||
while (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]);
|
||||
if ( (fabs(lat_first-lats[0]) >= d ) ||
|
||||
(fabs(lat_last+lats[0]) >= d ) ||
|
||||
lon_first != 0 ||
|
||||
fabs(lon_last - (360.0-90.0/order)) > 90.0/order
|
||||
fabs(lon_last - (360.0-360.0/max_pl)) > 360.0/max_pl
|
||||
) {
|
||||
int l=0;
|
||||
/*sub area*/
|
||||
/*find starting latitude */
|
||||
while (fabs(lat_first-lats[l]) > d ) {l++;}
|
||||
|
|
@ -215,4 +223,3 @@ static int destroy(grib_iterator* i)
|
|||
grib_context_free(c,self->los);
|
||||
return 1;
|
||||
}
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue