core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

ECC-341: implement area extraction in bufr_filter for uncompressed data

This commit is contained in:
Shahram Najm 2018-01-10 14:20:06 +00:00
parent b80884e7ca
commit 0b2b9d8195
2 changed files with 100 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

140
src/grib_accessor_class_bufr_extract_area_subsets.c
View File

@ -183,84 +183,102 @@ static int select_area(grib_accessor* a)
ret=grib_get_long(h,"compressedData",&compressed);
if (ret) return ret;
double *lat=0;
double *lon=0;
size_t n;
double lonWest,lonEast,latNorth,latSouth;
long numberOfSubsets,i,latRank,lonRank;
grib_iarray* subsets;
long *subsets_ar=0;
size_t nsubsets=0;
char latstr[20]={0,};
char lonstr[20]={0,};
ret=grib_get_long(h,self->numberOfSubsets,&numberOfSubsets);
if (ret) return ret;
subsets=grib_iarray_new(c,numberOfSubsets,10);
ret=grib_set_long(h,"unpack",1);
if (ret) return ret;
if (compressed) {
double *lat=0;
double *lon=0;
size_t n;
double lonWest,lonEast,latNorth,latSouth;
long numberOfSubsets,i,latRank,lonRank;
grib_iarray* subsets;
long *subsets_ar=0;
size_t nsubsets=0;
char latstr[20]={0,};
char lonstr[20]={0,};
ret=grib_get_long(h,self->numberOfSubsets,&numberOfSubsets);
if (ret) return ret;
subsets=grib_iarray_new(c,numberOfSubsets,10);
ret=grib_set_long(h,"unpack",1);
if (ret) return ret;
ret=grib_get_long(h,self->extractAreaLongitudeRank,&lonRank);
if (ret) return ret;
sprintf(lonstr,"#%ld#longitude",lonRank);
ret=grib_get_long(h,self->extractAreaLatitudeRank,&latRank);
if (ret) return ret;
sprintf(latstr,"#%ld#latitude",latRank);
}
n=numberOfSubsets;
lat=(double*)grib_context_malloc_clear(c,sizeof(double)*numberOfSubsets);
n=numberOfSubsets;
lat=(double*)grib_context_malloc_clear(c,sizeof(double)*numberOfSubsets);
if (compressed) {
ret=grib_get_double_array(h,latstr,lat,&n);
if (ret) return ret;
if (n!=numberOfSubsets) return GRIB_INTERNAL_ERROR;
if (n>numberOfSubsets) {
/* n could be less than the number of subsets if all latitudes are the same */
return GRIB_INTERNAL_ERROR;
}
} else {
for(i=0; i<numberOfSubsets; ++i) {
sprintf(latstr,"#%ld#latitude", i+1);
ret=grib_get_double(h,latstr,&(lat[i]));
if (ret) return ret;
}
}
lon=(double*)grib_context_malloc_clear(c,sizeof(double)*numberOfSubsets);
lon=(double*)grib_context_malloc_clear(c,sizeof(double)*numberOfSubsets);
if (compressed) {
ret=grib_get_double_array(h,lonstr,lon,&n);
if (ret) return ret;
if (n!=numberOfSubsets) return GRIB_INTERNAL_ERROR;
ret=grib_get_double(h,self->extractAreaWestLongitude,&lonWest);
if (ret) return ret;
ret=grib_get_double(h,self->extractAreaEastLongitude,&lonEast);
if (ret) return ret;
ret=grib_get_double(h,self->extractAreaNorthLatitude,&latNorth);
if (ret) return ret;
ret=grib_get_double(h,self->extractAreaSouthLatitude,&latSouth);
if (ret) return ret;
for (i=0;i<numberOfSubsets;i++) {
/* printf("++++++ lat: %g <= %g <= %g lon: %g <= %g <= %g \n",latSouth,lat[i],latNorth,lonWest,lon[i],lonEast); */
if (lat[i]>=latSouth && lat[i]<=latNorth && lon[i]>=lonWest && lon[i]<=lonEast) {
grib_iarray_push(subsets,i+1);
/* printf("++++++++ %ld\n",i+1); */
}
if (n>numberOfSubsets) {
/* n could be less than the number of subsets if all longitudes are the same */
return GRIB_INTERNAL_ERROR;
}
nsubsets=grib_iarray_used_size(subsets);
ret=grib_set_long(h,self->extractedAreaNumberOfSubsets,nsubsets);
if (ret) return ret;
if (nsubsets!=0) {
subsets_ar=grib_iarray_get_array(subsets);
ret=grib_set_long_array(h,self->extractSubsetList,subsets_ar,nsubsets);
if (ret) return ret;
ret=grib_set_long(h,self->doExtractSubsets,1);
if (ret) return ret;
}
grib_context_free(c,lat);
grib_context_free(c,lon);
grib_iarray_delete(subsets);
subsets=0;
} else {
grib_context_log(c, GRIB_LOG_ERROR, "Area extraction not implemented for uncompressed BUFR messages");
return GRIB_NOT_IMPLEMENTED;
for(i=0; i<numberOfSubsets; ++i) {
sprintf(lonstr,"#%ld#longitude", i+1);
ret=grib_get_double(h,lonstr,&(lon[i]));
if (ret) return ret;
}
}
ret=grib_get_double(h,self->extractAreaWestLongitude,&lonWest);
if (ret) return ret;
ret=grib_get_double(h,self->extractAreaEastLongitude,&lonEast);
if (ret) return ret;
ret=grib_get_double(h,self->extractAreaNorthLatitude,&latNorth);
if (ret) return ret;
ret=grib_get_double(h,self->extractAreaSouthLatitude,&latSouth);
if (ret) return ret;
for (i=0;i<n;i++) {
/* printf("++++++ lat: %g <= %g <= %g lon: %g <= %g <= %g \n",latSouth,lat[i],latNorth,lonWest,lon[i],lonEast); */
if (lat[i]>=latSouth && lat[i]<=latNorth && lon[i]>=lonWest && lon[i]<=lonEast) {
grib_iarray_push(subsets,i+1);
/* printf("++++++++ %ld\n",i+1); */
}
}
nsubsets=grib_iarray_used_size(subsets);
ret=grib_set_long(h,self->extractedAreaNumberOfSubsets,nsubsets);
if (ret) return ret;
if (nsubsets!=0) {
subsets_ar=grib_iarray_get_array(subsets);
ret=grib_set_long_array(h,self->extractSubsetList,subsets_ar,nsubsets);
if (ret) return ret;
ret=grib_set_long(h,self->doExtractSubsets,1);
if (ret) return ret;
}
grib_context_free(c,lat);
grib_context_free(c,lon);
grib_iarray_delete(subsets);
subsets=0;
return ret;
}

21
tests/bufr_filter_extract_area.sh
View File

@ -136,4 +136,25 @@ EOF
diff $outputRef $outputFilt
# Uncompressed message
# ---------------------
inputBufr="delayed_repl_01.bufr"
outputBufr=${label}.${inputBufr}.out
cat > $fRules <<EOF
transient originalNumberOfSubsets = numberOfSubsets;
set extractAreaNorthLatitude = -21.0;
set extractAreaSouthLatitude = -25.0;
set extractAreaWestLongitude = 136;
set extractAreaEastLongitude = 154;
set extractAreaLongitudeRank=1;
set doExtractArea=1;
write;
print "extracted [numberOfSubsets] of [originalNumberOfSubsets] subsets";
assert(3 == extractedAreaNumberOfSubsets);
EOF
${tools_dir}/codes_bufr_filter -o $outputBufr $fRules $inputBufr
ns=`${tools_dir}/bufr_get -p numberOfSubsets $outputBufr`
[ $ns -eq 3 ]
rm -f $outputRef $outputFilt $outputBufr $fLog $fRules