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

GRIB: nearest function on oblate spheroid

This commit is contained in:
Shahram Najm 2021-11-06 12:35:12 +00:00
parent 040ab8d66d
commit f575a0468c
6 changed files with 72 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/grib_api_prototypes.h
View File

@ -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,

50
src/grib_nearest.c
View File

@ -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) {

13
src/grib_nearest_class_latlon_reduced.c
View File

@ -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++;
}

13
src/grib_nearest_class_reduced.c
View File

@ -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++;
}

32
src/grib_nearest_class_regular.c
View File

@ -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++;
}

25
tests/grib_nearest_test.sh
View File

@ -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