From a27688ddc038590b96f5fa03e2e14a996a4a7e5e Mon Sep 17 00:00:00 2001 From: Shahram Najm Date: Thu, 13 Aug 2015 15:41:00 +0100 Subject: [PATCH] Issue error on rotated grids rather than incorrect results --- src/grib_nearest_class_regular.c | 581 ++++++++++++++++--------------- 1 file changed, 297 insertions(+), 284 deletions(-) diff --git a/src/grib_nearest_class_regular.c b/src/grib_nearest_class_regular.c index 00a089eb3..c5e3fbb92 100644 --- a/src/grib_nearest_class_regular.c +++ b/src/grib_nearest_class_regular.c @@ -92,319 +92,332 @@ static void init_class(grib_nearest_class* c) static int init(grib_nearest* nearest,grib_handle* h,grib_arguments* args) { - grib_nearest_regular* self = (grib_nearest_regular*) nearest; - self->Ni = grib_arguments_get_name(h,args,self->cargs++); - self->Nj = grib_arguments_get_name(h,args,self->cargs++); - self->i=(int*)grib_context_malloc(h->context,2*sizeof(int)); - self->j=(int*)grib_context_malloc(h->context,2*sizeof(int)); - return 0; + grib_nearest_regular* self = (grib_nearest_regular*) nearest; + self->Ni = grib_arguments_get_name(h,args,self->cargs++); + self->Nj = grib_arguments_get_name(h,args,self->cargs++); + self->i=(int*)grib_context_malloc(h->context,2*sizeof(int)); + self->j=(int*)grib_context_malloc(h->context,2*sizeof(int)); + return 0; } #if 0 static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon,unsigned long flags, - double* outlats,double* outlons, - double *values,double *distances,int* indexes, size_t *len) { - grib_nearest_regular* self = (grib_nearest_regular*) nearest; - int ret=0,kk=0,ii=0,jj=0; - size_t nvalues=0; + double inlat, double inlon,unsigned long flags, + double* outlats,double* outlons, + double *values,double *distances,int* indexes, size_t *len) { + grib_nearest_regular* self = (grib_nearest_regular*) nearest; + int ret=0,kk=0,ii=0,jj=0; + size_t nvalues=0; - long iradius; - double radius; + long iradius; + double radius; + + if( (ret = grib_get_long(h,self->radius,&iradius))!= GRIB_SUCCESS) + return ret; + radius=((double)iradius)/1000.0; + + if (!nearest->h || (flags & GRIB_NEAREST_SAME_DATA)==0 || nearest->h!=h) { + grib_iterator* iter=NULL; + double lat=0,lon=0; + + if( (ret = grib_get_size(h,self->values_key,&nvalues))!= GRIB_SUCCESS) + return ret; + nearest->values_count = nvalues; + if (nearest->values) grib_context_free(nearest->context,nearest->values); + nearest->values = grib_context_malloc(h->context,nvalues*sizeof(double)); + if (!nearest->values) return GRIB_OUT_OF_MEMORY; + + ret=grib_get_double_array_internal( h,self->values_key, + nearest->values,&(nearest->values_count)); + if (ret!=GRIB_SUCCESS) grib_context_log(nearest->context,GRIB_LOG_ERROR, + "nearest: unable to get values array"); + + if (!nearest->h || (flags & GRIB_NEAREST_SAME_GRID)==0) { + double dummy=0; + double olat=1.e10, olon=1.e10; + int ilat=0,ilon=0; + long n=0; + + if( (ret = grib_get_long(h,self->Ni,&n))!= GRIB_SUCCESS) + return ret; + self->lons_count=n; + + if( (ret = grib_get_long(h,self->Nj,&n))!= GRIB_SUCCESS) + return ret; + self->lats_count=n; + + if (self->lats) grib_context_free(nearest->context,self->lats); + self->lats=grib_context_malloc( nearest->context, + self->lats_count* sizeof(double)); + if (!self->lats) return GRIB_OUT_OF_MEMORY; + + if (self->lons) grib_context_free(nearest->context,self->lons); + self->lons=grib_context_malloc( nearest->context, + self->lons_count*sizeof(double)); + if (!self->lons) return GRIB_OUT_OF_MEMORY; + + iter=grib_iterator_new(h,0,&ret); + if (ret) { + grib_context_log(nearest->context,GRIB_LOG_ERROR,"unable to create iterator"); + return ret; + } + while(grib_iterator_next(iter,&lat,&lon,&dummy)) { + if (olat != lat) { + self->lats[ilat++]=lat; + olat=lat; + } + if (ilonlons_count && olon != lon) { + self->lons[ilon++]=lon; + olon=lon; + } + } + grib_iterator_delete(iter); + + } + nearest->h=h; + + } + + if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT)==0 + || (flags & GRIB_NEAREST_SAME_GRID)==0) { + + grib_binary_search(self->lats,self->lats_count-1,inlat, + &(self->j[0]),&(self->j[1])); + grib_binary_search(self->lons,self->lons_count-1,inlon, + &(self->i[0]),&(self->i[1])); + if (!self->distances) + self->distances=(double*)grib_context_malloc( nearest->context,4*sizeof(double)); + if (!self->k) + self->k=(int*)grib_context_malloc( nearest->context,4*sizeof(int)); + kk=0; + for (ii=0;ii<2;ii++) { + for (jj=0;jj<2;jj++) { + self->k[kk]=self->i[ii]+self->lons_count*self->j[jj]-1; + self->distances[kk]=grib_nearest_distance(radius,inlon,inlat, + self->lons[self->i[ii]],self->lats[self->j[jj]]); + kk++; + } + } + } + + kk=0; + for (ii=0;ii<2;ii++) { + for (jj=0;jj<2;jj++) { + distances[kk]=self->distances[kk]; + outlats[kk]=self->lats[self->j[jj]]; + outlons[kk]=self->lons[self->i[ii]]; + values[kk]=nearest->values[self->k[kk]]; + indexes[kk]=self->k[kk]; + kk++; + } + } + + return GRIB_SUCCESS; +} +#else +static int is_rotated_grid(grib_handle* h) +{ + long is_rotated = 0; + int err = grib_get_long(h, "is_rotated_grid", &is_rotated); + if (!err && is_rotated ) return 1; + return 0; +} + +static int find(grib_nearest* nearest, grib_handle* h, + double inlat, double inlon,unsigned long flags, + double* outlats,double* outlons, + double *values,double *distances,int* indexes, size_t *len) +{ + grib_nearest_regular* self = (grib_nearest_regular*) nearest; + int ret=0,kk=0,ii=0,jj=0; + size_t nvalues=0; + long iradius; + double radius; - if( (ret = grib_get_long(h,self->radius,&iradius))!= GRIB_SUCCESS) - return ret; - radius=((double)iradius)/1000.0; - - if (!nearest->h || (flags & GRIB_NEAREST_SAME_DATA)==0 || nearest->h!=h) { grib_iterator* iter=NULL; double lat=0,lon=0; - if( (ret = grib_get_size(h,self->values_key,&nvalues))!= GRIB_SUCCESS) - return ret; - nearest->values_count = nvalues; - if (nearest->values) grib_context_free(nearest->context,nearest->values); - nearest->values = grib_context_malloc(h->context,nvalues*sizeof(double)); - if (!nearest->values) return GRIB_OUT_OF_MEMORY; + while (inlon<0) inlon+=360; + while (inlon>360) inlon-=360; - ret=grib_get_double_array_internal( h,self->values_key, - nearest->values,&(nearest->values_count)); - if (ret!=GRIB_SUCCESS) grib_context_log(nearest->context,GRIB_LOG_ERROR, - "nearest: unable to get values array"); + if( (ret = grib_get_size(h,self->values_key,&nvalues))!= GRIB_SUCCESS) + return ret; + nearest->values_count = nvalues; + + if (grib_is_missing(h,self->radius,&ret)) { + grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->radius); + return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; + } + + if( (ret = grib_get_long(h,self->radius,&iradius))!= GRIB_SUCCESS) + return ret; + radius=((double)iradius)/1000.0; if (!nearest->h || (flags & GRIB_NEAREST_SAME_GRID)==0) { - double dummy=0; - double olat=1.e10, olon=1.e10; - int ilat=0,ilon=0; - long n=0; + double dummy=0; + double olat=1.e10, olon=1.e10; + int ilat=0,ilon=0; + long n=0; - if( (ret = grib_get_long(h,self->Ni,&n))!= GRIB_SUCCESS) - return ret; - self->lons_count=n; - - if( (ret = grib_get_long(h,self->Nj,&n))!= GRIB_SUCCESS) - return ret; - self->lats_count=n; - - if (self->lats) grib_context_free(nearest->context,self->lats); - self->lats=grib_context_malloc( nearest->context, - self->lats_count* sizeof(double)); - if (!self->lats) return GRIB_OUT_OF_MEMORY; - - if (self->lons) grib_context_free(nearest->context,self->lons); - self->lons=grib_context_malloc( nearest->context, - self->lons_count*sizeof(double)); - if (!self->lons) return GRIB_OUT_OF_MEMORY; - - iter=grib_iterator_new(h,0,&ret); - if (ret) { - grib_context_log(nearest->context,GRIB_LOG_ERROR,"unable to create iterator"); - return ret; - } - while(grib_iterator_next(iter,&lat,&lon,&dummy)) { - if (olat != lat) { - self->lats[ilat++]=lat; - olat=lat; + if (grib_is_missing(h,self->Ni,&ret)) { + grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->Ni); + return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; } - if (ilonlons_count && olon != lon) { - self->lons[ilon++]=lon; - olon=lon; + + if (grib_is_missing(h,self->Nj,&ret)) { + grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->Nj); + return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; } - } - grib_iterator_delete(iter); + + /* Support for rotated grids not yet implemented */ + if (is_rotated_grid(h)) { + grib_context_log(h->context,GRIB_LOG_ERROR, + "Nearest neighbour functionality is not supported for rotated grids."); + return GRIB_NOT_IMPLEMENTED; + } + + if ((ret = grib_get_long(h,self->Ni,&n))!= GRIB_SUCCESS) + return ret; + self->lons_count=n; + + if ((ret = grib_get_long(h,self->Nj,&n))!= GRIB_SUCCESS) + return ret; + self->lats_count=n; + + if (self->lats) grib_context_free(nearest->context,self->lats); + self->lats=(double*)grib_context_malloc( nearest->context, + self->lats_count* sizeof(double)); + if (!self->lats) return GRIB_OUT_OF_MEMORY; + + if (self->lons) grib_context_free(nearest->context,self->lons); + self->lons=(double*)grib_context_malloc( nearest->context, + self->lons_count*sizeof(double)); + if (!self->lons) return GRIB_OUT_OF_MEMORY; + + iter=grib_iterator_new(h,0,&ret); + while(grib_iterator_next(iter,&lat,&lon,&dummy)) { + if (olat != lat) { + Assert( ilat < self->lats_count ); + self->lats[ilat++]=lat; + olat=lat; + } + if (ilonlons_count && olon != lon) { + self->lons[ilon++]=lon ; + olon=lon; + } + } + grib_iterator_delete(iter); } nearest->h=h; - } + if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT)==0 + || (flags & GRIB_NEAREST_SAME_GRID)==0) { + int nearest_lons_found=0; - if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT)==0 - || (flags & GRIB_NEAREST_SAME_GRID)==0) { + if (self->lats[self->lats_count-1] > self->lats[0]) { + if (inlatlats[0] || inlat>self->lats[self->lats_count-1]) + return GRIB_OUT_OF_AREA; + } else { + if (inlat > self->lats[0] || inlat < self->lats[self->lats_count-1]) + return GRIB_OUT_OF_AREA; + } + + if (self->lons[self->lons_count-1] > self->lons[0]) { + if (inlonlons[0] || inlon>self->lons[self->lons_count-1]) { + /* try to scale*/ + if (inlon>0) inlon-=360; + else inlon+=360; + + if (inlonlons[0] || inlon>self->lons[self->lons_count-1]) { + if ( self->lons[0]+360-self->lons[self->lons_count-1]<= + self->lons[1]-self->lons[0]) { + /*it's a global field in longitude*/ + self->i[0]=0; + self->i[1]=self->lons_count-1; + nearest_lons_found=1; + } else + return GRIB_OUT_OF_AREA; + } + } + } else { + if (inlon>self->lons[0] || inlonlons[self->lons_count-1]) { + /* try to scale*/ + if (inlon>0) inlon-=360; + else inlon+=360; + if (self->lons[0]-self->lons[self->lons_count-1]-360 <= + self->lons[0]-self->lons[1]) { + /*it's a global field in longitude*/ + self->i[0]=0; + self->i[1]=self->lons_count-1; + nearest_lons_found=1; + } else if (inlon>self->lons[0] || inlonlons[self->lons_count-1]) + return GRIB_OUT_OF_AREA; + } + } + + grib_binary_search(self->lats,self->lats_count-1,inlat, + &(self->j[0]),&(self->j[1])); + + if (!nearest_lons_found) + grib_binary_search(self->lons,self->lons_count-1,inlon, + &(self->i[0]),&(self->i[1])); + + if (!self->distances) + self->distances=(double*)grib_context_malloc( nearest->context,4*sizeof(double)); + if (!self->k) + self->k=(int*)grib_context_malloc( nearest->context,4*sizeof(int)); + kk=0; + for (jj=0;jj<2;jj++) { + for (ii=0;ii<2;ii++) { + self->k[kk]=self->i[ii]+self->lons_count*self->j[jj]; + self->distances[kk]=grib_nearest_distance(radius,inlon,inlat, + self->lons[self->i[ii]],self->lats[self->j[jj]]); + kk++; + } + } + } - grib_binary_search(self->lats,self->lats_count-1,inlat, - &(self->j[0]),&(self->j[1])); - grib_binary_search(self->lons,self->lons_count-1,inlon, - &(self->i[0]),&(self->i[1])); - if (!self->distances) - self->distances=(double*)grib_context_malloc( nearest->context,4*sizeof(double)); - if (!self->k) - self->k=(int*)grib_context_malloc( nearest->context,4*sizeof(int)); kk=0; - for (ii=0;ii<2;ii++) { - for (jj=0;jj<2;jj++) { - self->k[kk]=self->i[ii]+self->lons_count*self->j[jj]-1; - self->distances[kk]=grib_nearest_distance(radius,inlon,inlat, - self->lons[self->i[ii]],self->lats[self->j[jj]]); - kk++; - } - } - } - kk=0; - for (ii=0;ii<2;ii++) { + /* + * Brute force algorithm: + * First unpack all the values into an array. Then when we need the 4 points + * we just index into this array so no need to call grib_get_double_element_internal + * + * if (nearest->values) grib_context_free(nearest->context,nearest->values); + * nearest->values = grib_context_malloc(h->context,nvalues*sizeof(double)); + * if (!nearest->values) return GRIB_OUT_OF_MEMORY; + * ret = grib_get_double_array(h, self->values_key, nearest->values ,&nvalues); + * if (ret) return ret; + */ + for (jj=0;jj<2;jj++) { - distances[kk]=self->distances[kk]; - outlats[kk]=self->lats[self->j[jj]]; - outlons[kk]=self->lons[self->i[ii]]; - values[kk]=nearest->values[self->k[kk]]; - indexes[kk]=self->k[kk]; - kk++; - } - } + for (ii=0;ii<2;ii++) { + distances[kk]=self->distances[kk]; + outlats[kk]=self->lats[self->j[jj]]; + outlons[kk]=self->lons[self->i[ii]]; + 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]]; */ - return GRIB_SUCCESS; -} -#else -static int find(grib_nearest* nearest, grib_handle* h, - double inlat, double inlon,unsigned long flags, - double* outlats,double* outlons, - double *values,double *distances,int* indexes, size_t *len) -{ - grib_nearest_regular* self = (grib_nearest_regular*) nearest; - int ret=0,kk=0,ii=0,jj=0; - size_t nvalues=0; - long iradius; - double radius; - - grib_iterator* iter=NULL; - double lat=0,lon=0; - - while (inlon<0) inlon+=360; - while (inlon>360) inlon-=360; - - if( (ret = grib_get_size(h,self->values_key,&nvalues))!= GRIB_SUCCESS) - return ret; - nearest->values_count = nvalues; - - if (grib_is_missing(h,self->radius,&ret)) { - grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->radius); - return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; - } - - if( (ret = grib_get_long(h,self->radius,&iradius))!= GRIB_SUCCESS) - return ret; - radius=((double)iradius)/1000.0; - - if (!nearest->h || (flags & GRIB_NEAREST_SAME_GRID)==0) { - double dummy=0; - double olat=1.e10, olon=1.e10; - int ilat=0,ilon=0; - long n=0; - - if (grib_is_missing(h,self->Ni,&ret)) { - grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->Ni); - return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; + indexes[kk]=self->k[kk]; + kk++; + } } - if (grib_is_missing(h,self->Nj,&ret)) { - grib_context_log(h->context, GRIB_LOG_DEBUG,"Key '%s' is missing", self->Nj); - return ret ? ret : GRIB_GEOCALCULUS_PROBLEM; - } - - if ((ret = grib_get_long(h,self->Ni,&n))!= GRIB_SUCCESS) - return ret; - self->lons_count=n; - - if ((ret = grib_get_long(h,self->Nj,&n))!= GRIB_SUCCESS) - return ret; - self->lats_count=n; - - if (self->lats) grib_context_free(nearest->context,self->lats); - self->lats=(double*)grib_context_malloc( nearest->context, - self->lats_count* sizeof(double)); - if (!self->lats) return GRIB_OUT_OF_MEMORY; - - if (self->lons) grib_context_free(nearest->context,self->lons); - self->lons=(double*)grib_context_malloc( nearest->context, - self->lons_count*sizeof(double)); - if (!self->lons) return GRIB_OUT_OF_MEMORY; - - iter=grib_iterator_new(h,0,&ret); - while(grib_iterator_next(iter,&lat,&lon,&dummy)) { - if (olat != lat) { - Assert( ilat < self->lats_count ); - self->lats[ilat++]=lat; - olat=lat; - } - if (ilonlons_count && olon != lon) { - self->lons[ilon++]=lon ; - olon=lon; - } - } - grib_iterator_delete(iter); - - } - nearest->h=h; - - if (!self->distances || (flags & GRIB_NEAREST_SAME_POINT)==0 - || (flags & GRIB_NEAREST_SAME_GRID)==0) { - int nearest_lons_found=0; - - if (self->lats[self->lats_count-1] > self->lats[0]) { - if (inlatlats[0] || inlat>self->lats[self->lats_count-1]) - return GRIB_OUT_OF_AREA; - } else { - if (inlat > self->lats[0] || inlat < self->lats[self->lats_count-1]) - return GRIB_OUT_OF_AREA; - } - - if (self->lons[self->lons_count-1] > self->lons[0]) { - if (inlonlons[0] || inlon>self->lons[self->lons_count-1]) { - /* try to scale*/ - if (inlon>0) inlon-=360; - else inlon+=360; - - if (inlonlons[0] || inlon>self->lons[self->lons_count-1]) { - if ( self->lons[0]+360-self->lons[self->lons_count-1]<= - self->lons[1]-self->lons[0]) { - /*it's a global field in longitude*/ - self->i[0]=0; - self->i[1]=self->lons_count-1; - nearest_lons_found=1; - } else - return GRIB_OUT_OF_AREA; - } - } - } else { - if (inlon>self->lons[0] || inlonlons[self->lons_count-1]) { - /* try to scale*/ - if (inlon>0) inlon-=360; - else inlon+=360; - if (self->lons[0]-self->lons[self->lons_count-1]-360 <= - self->lons[0]-self->lons[1]) { - /*it's a global field in longitude*/ - self->i[0]=0; - self->i[1]=self->lons_count-1; - nearest_lons_found=1; - } else if (inlon>self->lons[0] || inlonlons[self->lons_count-1]) - return GRIB_OUT_OF_AREA; - } - } - - grib_binary_search(self->lats,self->lats_count-1,inlat, - &(self->j[0]),&(self->j[1])); - - if (!nearest_lons_found) - grib_binary_search(self->lons,self->lons_count-1,inlon, - &(self->i[0]),&(self->i[1])); - - if (!self->distances) - self->distances=(double*)grib_context_malloc( nearest->context,4*sizeof(double)); - if (!self->k) - self->k=(int*)grib_context_malloc( nearest->context,4*sizeof(int)); - kk=0; - for (jj=0;jj<2;jj++) { - for (ii=0;ii<2;ii++) { - self->k[kk]=self->i[ii]+self->lons_count*self->j[jj]; - self->distances[kk]=grib_nearest_distance(radius,inlon,inlat, - self->lons[self->i[ii]],self->lats[self->j[jj]]); - kk++; - } - } - } - - kk=0; - -/* - * Brute force algorithm: - * First unpack all the values into an array. Then when we need the 4 points - * we just index into this array so no need to call grib_get_double_element_internal - * - * if (nearest->values) grib_context_free(nearest->context,nearest->values); - * nearest->values = grib_context_malloc(h->context,nvalues*sizeof(double)); - * if (!nearest->values) return GRIB_OUT_OF_MEMORY; - * ret = grib_get_double_array(h, self->values_key, nearest->values ,&nvalues); - * if (ret) return ret; - */ - - for (jj=0;jj<2;jj++) { - for (ii=0;ii<2;ii++) { - distances[kk]=self->distances[kk]; - outlats[kk]=self->lats[self->j[jj]]; - outlons[kk]=self->lons[self->i[ii]]; - 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]]; */ - - indexes[kk]=self->k[kk]; - kk++; - } - } - - return GRIB_SUCCESS; + return GRIB_SUCCESS; } #endif - static int destroy(grib_nearest* nearest) { - grib_nearest_regular* self = (grib_nearest_regular*) nearest; - if (self->lats) grib_context_free(nearest->context,self->lats); - if (self->lons) grib_context_free(nearest->context,self->lons); - if (self->i) grib_context_free(nearest->context,self->i); - if (self->j) grib_context_free(nearest->context,self->j); - if (self->k) grib_context_free(nearest->context,self->k); - if (self->distances) grib_context_free(nearest->context,self->distances); - return GRIB_SUCCESS; + grib_nearest_regular* self = (grib_nearest_regular*) nearest; + if (self->lats) grib_context_free(nearest->context,self->lats); + if (self->lons) grib_context_free(nearest->context,self->lons); + if (self->i) grib_context_free(nearest->context,self->i); + if (self->j) grib_context_free(nearest->context,self->j); + if (self->k) grib_context_free(nearest->context,self->k); + if (self->distances) grib_context_free(nearest->context,self->distances); + return GRIB_SUCCESS; } -