ECC-478: Scaling values using grib_set ignores values=9999. when bitmapPresent=0

This commit is contained in:
Shahram Najm 2017-07-31 15:40:31 +01:00
parent 52557b0290
commit 061e84ff92
5 changed files with 138 additions and 94 deletions

View File

@ -137,60 +137,63 @@ static void init_class(grib_accessor_class* c)
static void init(grib_accessor* a,const long l, grib_arguments* args)
{
int n=0;
grib_accessor_offset_values* self= (grib_accessor_offset_values*)a;
self->values=grib_arguments_get_name(grib_handle_of_accessor(a),args,n++);
self->missingValue=grib_arguments_get_name(grib_handle_of_accessor(a),args,n++);
a->flags |= GRIB_ACCESSOR_FLAG_FUNCTION;
a->length=0;
int n=0;
grib_accessor_offset_values* self= (grib_accessor_offset_values*)a;
self->values=grib_arguments_get_name(grib_handle_of_accessor(a),args,n++);
self->missingValue=grib_arguments_get_name(grib_handle_of_accessor(a),args,n++);
a->flags |= GRIB_ACCESSOR_FLAG_FUNCTION;
a->length=0;
}
static int unpack_double(grib_accessor* a, double* val, size_t *len)
{
int ret=0;
*val=0;
*len =1;
return ret;
int ret=0;
*val=0;
*len =1;
return ret;
}
static int pack_double(grib_accessor* a, const double* val, size_t *len)
{
double* values=NULL;
size_t size=0;
double missingValue=0;
int ret=0,i=0;
grib_accessor_offset_values* self= (grib_accessor_offset_values*)a;
grib_context* c=a->context;
grib_handle* h=grib_handle_of_accessor(a);
double* values=NULL;
size_t size=0;
double missingValue=0;
long missingValuesPresent = 0;
int ret=0,i=0;
grib_accessor_offset_values* self= (grib_accessor_offset_values*)a;
grib_context* c=a->context;
grib_handle* h=grib_handle_of_accessor(a);
if (*val==0) return GRIB_SUCCESS;
if (*val==0) return GRIB_SUCCESS;
if((ret = grib_get_double_internal(h,self->missingValue,&missingValue))
!= GRIB_SUCCESS) {
if((ret = grib_get_double_internal(h,self->missingValue,&missingValue)) != GRIB_SUCCESS) {
return ret;
}
}
if((ret=grib_get_long_internal(h,"missingValuesPresent",&missingValuesPresent)) != GRIB_SUCCESS) {
return ret;
}
if ( (ret=grib_get_size(h,self->values,&size)) != GRIB_SUCCESS) return ret;
if ( (ret=grib_get_size(h,self->values,&size)) != GRIB_SUCCESS) return ret;
values=(double*)grib_context_malloc(c,size*sizeof(double));
if (!values) return GRIB_OUT_OF_MEMORY;
values=(double*)grib_context_malloc(c,size*sizeof(double));
if (!values) return GRIB_OUT_OF_MEMORY;
if((ret = grib_get_double_array_internal(h,self->values,values,&size))
!= GRIB_SUCCESS) {
if((ret = grib_get_double_array_internal(h,self->values,values,&size)) != GRIB_SUCCESS) {
grib_context_free(c,values);
return ret;
}
}
for (i=0;i<size;i++)
if (values[i]!=missingValue) values[i]+=*val;
for (i=0;i<size;i++) {
if (missingValuesPresent) {
if (values[i]!=missingValue) values[i]+=*val;
} else {
values[i]+=*val;
}
}
if((ret = grib_set_double_array_internal(h, self->values,values,size)) != GRIB_SUCCESS) return ret;
if((ret = grib_set_double_array_internal(h, self->values,values,size))
!= GRIB_SUCCESS) return ret;
grib_context_free(c,values);
grib_context_free(c,values);
return GRIB_SUCCESS;
return GRIB_SUCCESS;
}

View File

@ -137,59 +137,63 @@ static void init_class(grib_accessor_class* c)
static void init(grib_accessor* a,const long l, grib_arguments* args)
{
int n=0;
grib_accessor_scale_values* self= (grib_accessor_scale_values*)a;
self->values=grib_arguments_get_name(grib_handle_of_accessor(a),args,n++);
self->missingValue=grib_arguments_get_name(grib_handle_of_accessor(a),args,n++);
a->flags |= GRIB_ACCESSOR_FLAG_FUNCTION;
a->length=0;
int n=0;
grib_accessor_scale_values* self= (grib_accessor_scale_values*)a;
self->values=grib_arguments_get_name(grib_handle_of_accessor(a),args,n++);
self->missingValue=grib_arguments_get_name(grib_handle_of_accessor(a),args,n++);
a->flags |= GRIB_ACCESSOR_FLAG_FUNCTION;
a->length=0;
}
static int unpack_double(grib_accessor* a, double* val, size_t *len)
static int unpack_double(grib_accessor* a, double* val, size_t *len)
{
int ret=0;
*val=1;
*len =1;
return ret;
int ret=0;
*val=1;
*len =1;
return ret;
}
static int pack_double(grib_accessor* a, const double* val, size_t *len)
{
double* values=NULL;
double missingValue=0;
size_t size=0;
int ret=0,i=0;
grib_accessor_scale_values* self= (grib_accessor_scale_values*)a;
grib_context* c=a->context;
grib_handle* h=grib_handle_of_accessor(a);
double* values=NULL;
double missingValue=0;
long missingValuesPresent = 0;
size_t size=0;
int ret=0,i=0;
grib_accessor_scale_values* self= (grib_accessor_scale_values*)a;
grib_context* c=a->context;
grib_handle* h=grib_handle_of_accessor(a);
if (*val==1) return GRIB_SUCCESS;
if (*val==1) return GRIB_SUCCESS;
if((ret = grib_get_double_internal(h,self->missingValue,&missingValue))
!= GRIB_SUCCESS) {
if((ret = grib_get_double_internal(h,self->missingValue,&missingValue)) != GRIB_SUCCESS) {
return ret;
}
}
if((ret=grib_get_long_internal(h,"missingValuesPresent",&missingValuesPresent)) != GRIB_SUCCESS) {
return ret;
}
if ( (ret=grib_get_size(h,self->values,&size)) != GRIB_SUCCESS) return ret;
if ( (ret=grib_get_size(h,self->values,&size)) != GRIB_SUCCESS) return ret;
values=(double*)grib_context_malloc(c,size*sizeof(double));
if (!values) return GRIB_OUT_OF_MEMORY;
values=(double*)grib_context_malloc(c,size*sizeof(double));
if (!values) return GRIB_OUT_OF_MEMORY;
if((ret = grib_get_double_array_internal(h,self->values,values,&size))
!= GRIB_SUCCESS) {
if((ret = grib_get_double_array_internal(h,self->values,values,&size)) != GRIB_SUCCESS) {
grib_context_free(c,values);
return ret;
}
}
for (i=0;i<size;i++)
if (values[i]!=missingValue) values[i]*=*val;
for (i=0;i<size;i++) {
if (missingValuesPresent) {
if (values[i]!=missingValue) values[i]*=*val;
} else {
values[i]*=*val;
}
}
if((ret = grib_set_double_array_internal(h, self->values,values,size))
!= GRIB_SUCCESS) return ret;
if((ret = grib_set_double_array_internal(h, self->values,values,size)) != GRIB_SUCCESS) return ret;
grib_context_free(c,values);
grib_context_free(c,values);
return GRIB_SUCCESS;
return GRIB_SUCCESS;
}

View File

@ -170,6 +170,7 @@ static int unpack_double(grib_accessor* a, double* val, size_t *len)
size_t size=0, real_size=0;
double max,min,avg,sd,value,skew,kurt, m2=0,m3=0,m4=0;
double missing=0;
long missingValuesPresent = 0;
size_t number_of_missing=0;
grib_context* c=a->context;
grib_handle* h=grib_handle_of_accessor(a);
@ -181,32 +182,43 @@ static int unpack_double(grib_accessor* a, double* val, size_t *len)
grib_context_log(a->context,GRIB_LOG_DEBUG,
"grib_accessor_statistics: computing statistics for %d values",size);
if((ret=grib_get_double(grib_handle_of_accessor(a),self->missing_value,&missing))
!= GRIB_SUCCESS) return ret;
if((ret=grib_get_double(h,self->missing_value,&missing)) != GRIB_SUCCESS) return ret;
if((ret=grib_get_long_internal(h,"missingValuesPresent",&missingValuesPresent)) != GRIB_SUCCESS) return ret;
values=(double*)grib_context_malloc_clear(c,size*sizeof(double));
if (!values) return GRIB_OUT_OF_MEMORY;
if((ret = grib_get_double_array_internal(h,self->values,values,&size))
!= GRIB_SUCCESS) {
if((ret = grib_get_double_array_internal(h,self->values,values,&size)) != GRIB_SUCCESS) {
grib_context_free(c,values);
return ret;
}
number_of_missing=0;
i=0;
while (i<size && values[i] == missing ) {i++;number_of_missing++;}
max=values[i];
min=values[i];
avg=values[i];
for (i=number_of_missing+1;i<size;i++) {
value=values[i];
if (value > max && value != missing) max=value;
if (value < min && value != missing) min=value;
if (value != missing) avg+=value;
else number_of_missing++;
if (missingValuesPresent) {
i=0;
while (i<size && values[i] == missing) {i++;number_of_missing++;}
max=values[i];
min=values[i];
avg=values[i];
for (i=number_of_missing+1;i<size;i++) {
value=values[i];
if (value > max && value != missing) max=value;
if (value < min && value != missing) min=value;
if (value != missing) avg+=value;
else number_of_missing++;
}
} else {
max=values[0];
min=values[0];
avg=values[0];
for (i=1; i<size; i++) {
value = values[i];
if (value > max) max=value;
if (value < min) min=value;
avg += value;
}
}
/*printf("stats.......... number_of_missing=%ld\n", number_of_missing);*/
/* Don't divide by zero if all values are missing! */
if (size != number_of_missing) {
avg/=(size-number_of_missing);
@ -214,7 +226,11 @@ static int unpack_double(grib_accessor* a, double* val, size_t *len)
sd=0; skew=0; kurt=0;
for (i=0;i<size;i++) {
if (values[i] != missing) {
int valueNotMissing = 1;
if (missingValuesPresent) {
valueNotMissing = (values[i] != missing);
}
if (valueNotMissing) {
double v = values[i] - avg;
double tmp = v*v;
m2 += tmp;

View File

@ -35,5 +35,26 @@ ${tools_dir}/grib_filter statistics.filter ${data_dir}/$file > statistics.out
diff statistics.out ${data_dir}/statistics.out.good
done
rm -f statistics.out statistics.filter
rm -f statistics.out statistics.filter || true
# GRIB with no missing values but some entries = 9999
# See ECC-478
# ---------------------------------------------------
input=${data_dir}/lfpw.grib1
temp1=temp1.statistics.grib
temp2=temp2.statistics.grib
stats=`${tools_dir}/grib_get -w count=50 -F%.2f -n statistics $input`
[ "$stats" = "10098 0 1064.19 3066.07 2.57004 4.60965 0" ]
# Scaling values in presence of real 9999 values
${tools_dir}/grib_set -s scaleValuesBy=0.5 $input $temp1
${tools_dir}/grib_set -s missingValue=1.0E34,scaleValuesBy=0.5 $input $temp2
${tools_dir}/grib_compare $temp1 $temp2
# Offsetting values in presence of real 9999 values
${tools_dir}/grib_set -s offsetValuesBy=0.5 $input $temp1
${tools_dir}/grib_set -s missingValue=1.0E34,offsetValuesBy=0.5 $input $temp2
${tools_dir}/grib_compare $temp1 $temp2
rm -f $temp1 $temp2

View File

@ -48,9 +48,9 @@ int main(int argc, const char *argv[])
while( (h = grib_handle_new_from_file(NULL,in,&e)) != NULL )
{
long bitmap;
long missingValuesPresent;
double delta;
GRIB_CHECK(grib_get_long(h,"bitMapIndicator",&bitmap),argv[i]);
GRIB_CHECK(grib_get_long(h,"missingValuesPresent",&missingValuesPresent),argv[i]);
GRIB_CHECK(grib_get_size(h,"values",&size),argv[i]);
@ -60,7 +60,7 @@ int main(int argc, const char *argv[])
values = (double*)malloc(size*sizeof(double));
last_size = size;
if (!values) {
fprintf(stderr, "Failed to allocate memory for values (size=%ld)\n", size);
fprintf(stderr, "Failed to allocate memory for values (%ld bytes)\n", size*sizeof(double));
exit(1);
}
}
@ -70,7 +70,7 @@ int main(int argc, const char *argv[])
for(j = 0; j < count; j++)
intervals[j] = 0;
if(bitmap != 255)
if(missingValuesPresent)
{
double missing;
GRIB_CHECK(grib_get_double(h,"missingValue",&missing),argv[i]);