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

ECC-1295: GRIB: nearest function fails on regular lat/lon grid on oblate spheroid

This commit is contained in:
Shahram Najm 2021-11-03 16:17:48 +00:00
parent 86c6ee5b09
commit ae1206473e
2 changed files with 26 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
src/grib_nearest_class_regular.c
View File

@ -232,7 +232,7 @@ static int find(grib_nearest* nearest, grib_handle* h,
int ret = 0, kk = 0, ii = 0, jj = 0;
size_t nvalues = 0;
long iradius;
double radius;
double radius; /* radius in km */
grib_iterator* iter = NULL;
double lat = 0, lon = 0;
@ -248,14 +248,24 @@ static int find(grib_nearest* nearest, grib_handle* h,
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;
/* We need the radius to calculate the nearest distance. For an oblate earth
approximate this using the average of the semimajor and semiminor axes */
if ((ret = grib_get_long(h, self->radius, &iradius)) == GRIB_SUCCESS) {
if (grib_is_missing(h, self->radius, &ret) || iradius == GRIB_MISSING_LONG) {
grib_context_log(h->context, GRIB_LOG_DEBUG, "Key '%s' is missing", self->radius);
return GRIB_GEOCALCULUS_PROBLEM;
}
radius = ((double)iradius) / 1000.0;
}
else {
double minor = 0, major = 0;
if ((ret = grib_get_double_internal(h, "earthMinorAxisInMetres", &minor)) != GRIB_SUCCESS) return ret;
if ((ret = grib_get_double_internal(h, "earthMajorAxisInMetres", &major)) != GRIB_SUCCESS) return ret;
if (grib_is_missing(h, "earthMinorAxisInMetres", &ret)) return GRIB_GEOCALCULUS_PROBLEM;
if (grib_is_missing(h, "earthMajorAxisInMetres", &ret)) return GRIB_GEOCALCULUS_PROBLEM;
radius = (major + minor) / 2.0;
radius = radius / 1000.0;
}
if ((ret = grib_get_long(h, self->radius, &iradius)) != GRIB_SUCCESS)
return ret;
radius = ((double)iradius) / 1000.0;
/* Compute lat/lon info, create iterator etc if it's the 1st time or different grid.
* This is for performance: if the grid has not changed, we only do this once

8
tests/grib_nearest_test.sh
View File

@ -51,5 +51,13 @@ Idx lat lon dist
EOF
diff $tempRef $temp
# ECC-1295: regular lat/lon on ellipsoid
# ----------------------------------------
sample2=$ECCODES_SAMPLES_PATH/GRIB2.tmpl
${tools_dir}/grib_set -s shapeOfTheEarth=5 $sample2 $temp
grib_check_key_equals $sample2 earthIsOblate 0
grib_check_key_equals $temp earthIsOblate 1
${tools_dir}/grib_ls -l 0,0 $temp
# Clean up
rm -f $temp $tempRef