GRIB-406: grib_get_data: Values outside of subarea specified

This commit is contained in:
Shahram Najm 2013-07-17 15:57:48 +01:00
parent 4dfda27c0b
commit da3ae73708
1 changed files with 13 additions and 12 deletions

View File

@ -90,15 +90,14 @@ static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){
grib_iterator_gaussian* self = (grib_iterator_gaussian*)i;
double *lats;
double laf;
double lal;
long trunc;
double laf; /* latitude of first point in degrees */
double lal; /* latitude of last point in degrees */
long trunc; /* number of parallels between a pole and the equator */
long lai;
long jScansPositively=0;
int size;
double start;
unsigned long istart=0;
int ret = GRIB_SUCCESS;
const char* latofirst = grib_arguments_get_name(h,args,self->carg++);
@ -106,7 +105,6 @@ static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){
const char* numtrunc = grib_arguments_get_name(h,args,self->carg++);
const char* s_jScansPositively = grib_arguments_get_name(h,args,self->carg++);
if((ret = grib_get_double_internal(h,latofirst, &laf))) return ret;
if((ret = grib_get_double_internal(h,latoflast, &lal))) return ret;
if((ret = grib_get_long_internal(h,numtrunc,&trunc))) return ret;
@ -128,14 +126,13 @@ static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){
/*
for(loi=(trunc*2)-1;loi>=0;loi--)
if(fabs(lats[loi] - lal) < glatPrecision) break;
for(j=(trunc*2)-1;j>0;j--) {
if(fabs(lats[j] - laf) < glatPrecision) break;
}
*/
binary_search(lats,size-1,start,&istart);
Assert(istart < size);
if (jScansPositively) {
for(lai=0;lai<self->nam;lai++) {
@ -152,21 +149,25 @@ static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){
grib_context_free(h->context,lats);
return ret;
}
static void binary_search(double xx[], const unsigned long n, double x, unsigned long *j)
{
/*These routine works only on descending ordered arrays*/
/*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 (x <= xx[jm]) jl=jm;
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;
}