From f575a0468c2b43f0cab66c11f2f04491ae72025b Mon Sep 17 00:00:00 2001 From: Shahram Najm Date: Sat, 6 Nov 2021 12:35:12 +0000 Subject: [PATCH] GRIB: nearest function on oblate spheroid --- src/grib_api_prototypes.h | 1 + src/grib_nearest.c | 50 +++++++++++++++++-------- src/grib_nearest_class_latlon_reduced.c | 13 ++----- src/grib_nearest_class_reduced.c | 13 ++----- src/grib_nearest_class_regular.c | 32 +++------------- tests/grib_nearest_test.sh | 25 +++++++++++++ 6 files changed, 72 insertions(+), 62 deletions(-) diff --git a/src/grib_api_prototypes.h b/src/grib_api_prototypes.h index a10dabe40..4e4cfa376 100644 --- a/src/grib_api_prototypes.h +++ b/src/grib_api_prototypes.h @@ -1380,6 +1380,7 @@ grib_box* grib_box_factory(grib_handle* h, grib_arguments* args); int grib_nearest_find(grib_nearest* nearest, const grib_handle* h, double inlat, double inlon, unsigned long flags, double* outlats, double* outlons, double* values, double* distances, int* indexes, size_t* len); int grib_nearest_init(grib_nearest* i, grib_handle* h, grib_arguments* args); int grib_nearest_delete(grib_nearest* i); +int grib_nearest_get_radius(grib_handle* h, double* radiusInKm); void grib_binary_search(double xx[], const unsigned long n, double x, int* ju, int* jl); int grib_nearest_find_multiple(const grib_handle* h, int is_lsm, const double* inlats, const double* inlons, long npoints, double* outlats, double* outlons, double* values, double* distances, int* indexes); int grib_nearest_find_generic(grib_nearest* nearest, grib_handle* h, double inlat, double inlon, unsigned long flags, diff --git a/src/grib_nearest.c b/src/grib_nearest.c index a07ab5904..f389c60b2 100644 --- a/src/grib_nearest.c +++ b/src/grib_nearest.c @@ -97,6 +97,37 @@ int grib_nearest_delete(grib_nearest* i) return 0; } +/* Get the radius in kilometres for nearest neighbour distance calculations */ +/* For an ellipsoid, approximate the radius using the average of the semimajor and semiminor axes */ +int grib_nearest_get_radius(grib_handle* h, double* radiusInKm) +{ + int err = 0; + long lRadiusInMetres; + double result = 0; + const char* s_radius = "radius"; + const char* s_minor = "earthMinorAxisInMetres"; + const char* s_major = "earthMajorAxisInMetres"; + + if ((err = grib_get_long(h, s_radius, &lRadiusInMetres)) == GRIB_SUCCESS) { + if (grib_is_missing(h, s_radius, &err) || lRadiusInMetres == GRIB_MISSING_LONG) { + grib_context_log(h->context, GRIB_LOG_DEBUG, "Key 'radius' is missing"); + return GRIB_GEOCALCULUS_PROBLEM; + } + result = ((double)lRadiusInMetres) / 1000.0; + } + else { + double minor = 0, major = 0; + if ((err = grib_get_double_internal(h, s_minor, &minor)) != GRIB_SUCCESS) return err; + if ((err = grib_get_double_internal(h, s_major, &major)) != GRIB_SUCCESS) return err; + if (grib_is_missing(h, s_minor, &err)) return GRIB_GEOCALCULUS_PROBLEM; + if (grib_is_missing(h, s_major, &err)) return GRIB_GEOCALCULUS_PROBLEM; + result = (major + minor) / 2.0; + result = result / 1000.0; + } + *radiusInKm = result; + return GRIB_SUCCESS; +} + void grib_binary_search(double xx[], const unsigned long n, double x, int* ju, int* jl) { @@ -278,7 +309,7 @@ int grib_nearest_find_generic( { int ret = 0, i = 0; size_t nvalues = 0, nneighbours = 0; - double radiusInMetres, radiusInKm; + double radiusInKm; grib_iterator* iter = NULL; double lat = 0, lon = 0; @@ -291,21 +322,8 @@ int grib_nearest_find_generic( return ret; nearest->values_count = nvalues; - /* 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_double(h, radius_keyname, &radiusInMetres)) == GRIB_SUCCESS && - !grib_is_missing(h, radius_keyname, &ret)) { - radiusInKm = radiusInMetres / 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; - radiusInMetres = (major + minor) / 2; - radiusInKm = radiusInMetres / 1000.0; - } + if ((ret = grib_nearest_get_radius(h, &radiusInKm)) != GRIB_SUCCESS) + return ret; neighbours = (PointStore*)grib_context_malloc(nearest->context, nvalues * sizeof(PointStore)); for (i = 0; i < nvalues; ++i) { diff --git a/src/grib_nearest_class_latlon_reduced.c b/src/grib_nearest_class_latlon_reduced.c index 5c4243718..02e31f6e8 100644 --- a/src/grib_nearest_class_latlon_reduced.c +++ b/src/grib_nearest_class_latlon_reduced.c @@ -119,22 +119,15 @@ static int find(grib_nearest* nearest, grib_handle* h, size_t nvalues = 0; grib_iterator* iter = NULL; double lat = 0, lon = 0; - long iradius; - double radius; + double radiusInKm; int ilat = 0, ilon = 0; 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) + if ((ret = grib_nearest_get_radius(h, &radiusInKm)) != 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 @@ -352,7 +345,7 @@ static int find(grib_nearest* nearest, grib_handle* h, kk = 0; for (jj = 0; jj < 2; jj++) { for (ii = 0; ii < 2; ii++) { - self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat, + self->distances[kk] = geographic_distance_spherical(radiusInKm, inlon, inlat, self->lons[self->k[kk]], self->lats[self->j[jj]]); kk++; } diff --git a/src/grib_nearest_class_reduced.c b/src/grib_nearest_class_reduced.c index 2a2235d1e..7fb2232f8 100644 --- a/src/grib_nearest_class_reduced.c +++ b/src/grib_nearest_class_reduced.c @@ -148,8 +148,7 @@ static int find(grib_nearest* nearest, grib_handle* h, size_t nvalues = 0; grib_iterator* iter = NULL; double lat = 0, lon = 0; - long iradius; - double radius; + double radiusInKm; int ilat = 0, ilon = 0; get_reduced_row_proc get_reduced_row_func = &grib_get_reduced_row; @@ -164,14 +163,8 @@ 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; - } - - if ((ret = grib_get_long(h, self->radius, &iradius)) != GRIB_SUCCESS) + if ((ret = grib_nearest_get_radius(h, &radiusInKm)) != 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 @@ -421,7 +414,7 @@ static int find(grib_nearest* nearest, grib_handle* h, kk = 0; for (jj = 0; jj < 2; jj++) { for (ii = 0; ii < 2; ii++) { - self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat, + self->distances[kk] = geographic_distance_spherical(radiusInKm, inlon, inlat, self->lons[self->k[kk]], self->lats[self->j[jj]]); kk++; } diff --git a/src/grib_nearest_class_regular.c b/src/grib_nearest_class_regular.c index 6849294cc..deea985d4 100644 --- a/src/grib_nearest_class_regular.c +++ b/src/grib_nearest_class_regular.c @@ -108,13 +108,10 @@ static int find(grib_nearest* nearest, grib_handle* h, grib_nearest_regular* self = (grib_nearest_regular*) nearest; int ret=0,kk=0,ii=0,jj=0; size_t nvalues=0; + double radiusInKm; - long iradius; - double radius; - - if( (ret = grib_get_long(h,self->radius,&iradius))!= GRIB_SUCCESS) + if ((ret = grib_nearest_get_radius(h, &radiusInKm)) != 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; @@ -231,8 +228,7 @@ static int find(grib_nearest* nearest, grib_handle* h, 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; /* radius in km */ + double radiusInKm; grib_iterator* iter = NULL; double lat = 0, lon = 0; @@ -248,24 +244,8 @@ static int find(grib_nearest* nearest, grib_handle* h, return ret; nearest->values_count = nvalues; - /* 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_nearest_get_radius(h, &radiusInKm)) != GRIB_SUCCESS) + return ret; /* 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 @@ -424,7 +404,7 @@ static int find(grib_nearest* nearest, grib_handle* h, 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] = geographic_distance_spherical(radius, inlon, inlat, + self->distances[kk] = geographic_distance_spherical(radiusInKm, inlon, inlat, self->lons[self->i[ii]], self->lats[self->j[jj]]); kk++; } diff --git a/tests/grib_nearest_test.sh b/tests/grib_nearest_test.sh index 03a35aa49..33bc2c16f 100755 --- a/tests/grib_nearest_test.sh +++ b/tests/grib_nearest_test.sh @@ -54,10 +54,35 @@ diff $tempRef $temp # ECC-1295: regular lat/lon on ellipsoid # ---------------------------------------- sample2=$ECCODES_SAMPLES_PATH/GRIB2.tmpl +${tools_dir}/grib_set -s shapeOfTheEarth=2 $sample2 $temp +grib_check_key_equals $sample2 earthIsOblate 0 +grib_check_key_equals $temp earthIsOblate 1 +${tools_dir}/grib_ls -l 0,0 $temp + +# reduced lat/lon on ellipsoid +# ---------------------------------------- +sample2=$ECCODES_SAMPLES_PATH/reduced_ll_sfc_grib2.tmpl +${tools_dir}/grib_set -s shapeOfTheEarth=4 $sample2 $temp +grib_check_key_equals $sample2 earthIsOblate 0 +grib_check_key_equals $temp earthIsOblate 1 +${tools_dir}/grib_ls -l 0,0 $temp + +# regular gaussian on ellipsoid +# ---------------------------------------- +sample2=$ECCODES_SAMPLES_PATH/regular_gg_pl_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 +# reduced gaussian on ellipsoid +# ---------------------------------------- +sample2=$ECCODES_SAMPLES_PATH/reduced_gg_pl_48_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