From defe48fe6c2b96a8bdf08ac716dcc19f06bf7d3e Mon Sep 17 00:00:00 2001 From: Shahram Najm Date: Sat, 13 Jun 2020 19:34:13 +0100 Subject: [PATCH] Geographic distance calculations (spherical and ellipsoid) --- src/grib_api_prototypes.h | 3 +- src/grib_nearest.c | 48 +++++++++++++++++++++++-- src/grib_nearest_class_latlon_reduced.c | 2 +- src/grib_nearest_class_reduced.c | 4 +-- src/grib_nearest_class_regular.c | 4 +-- 5 files changed, 53 insertions(+), 8 deletions(-) diff --git a/src/grib_api_prototypes.h b/src/grib_api_prototypes.h index bb8e8d182..165622bd8 100644 --- a/src/grib_api_prototypes.h +++ b/src/grib_api_prototypes.h @@ -1376,7 +1376,8 @@ int grib_nearest_find(grib_nearest* nearest, const grib_handle* h, double inlat, int grib_nearest_init(grib_nearest* i, grib_handle* h, grib_arguments* args); int grib_nearest_delete(grib_nearest* i); void grib_binary_search(double xx[], const unsigned long n, double x, int* ju, int* jl); -double grib_nearest_distance(double radius, double lon1, double lat1, double lon2, double lat2); +double grib_nearest_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2); +double grib_nearest_distance_ellipsoid(double major, double minor, double lon1, double lat1, double lon2, double lat2); 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, const char* values_keyname, const char* radius_keyname, diff --git a/src/grib_nearest.c b/src/grib_nearest.c index 81743fadd..0908457ca 100644 --- a/src/grib_nearest.c +++ b/src/grib_nearest.c @@ -116,7 +116,8 @@ void grib_binary_search(double xx[], const unsigned long n, double x, #define RADIAN(x) ((x)*acos(0.0) / 90.0) -double grib_nearest_distance(double radius, double lon1, double lat1, double lon2, double lat2) +/* radius is in km, angles in degrees */ +double grib_nearest_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2) { double rlat1 = RADIAN(lat1); double rlat2 = RADIAN(lat2); @@ -143,6 +144,49 @@ double grib_nearest_distance(double radius, double lon1, double lat1, double lon return radius * acos(a); } +/* major and minor axes in km, angles in degrees */ +double grib_nearest_distance_ellipsoid(double major, double minor, double lon1, double lat1, double lon2, double lat2) +{ + /* Lambert's formula */ + double rlat1 = RADIAN(lat1); + double rlat2 = RADIAN(lat2); + double rlon1 = RADIAN(lon1); + double rlon2 = RADIAN(lon2); + double deltaLat = rlat2 - rlat1; + double deltaLon = rlon2 - rlon1; + double sinDlat = sin(deltaLat/2.0); + double sinDlon = sin(deltaLon/2.0); + double sin2Dlat = sinDlat*sinDlat; + double sin2Dlon = sinDlon*sinDlon; + + double a = sin2Dlat + cos(rlat1) * cos(rlat2) * sin2Dlon; + double c = 2 * atan2(sqrt(a), sqrt(1.0-a)); + double f = (major - minor)/major; /*flattening*/ + + double latr1 = atan( (1.0-f)*tan(rlat1) ); /*Reduced latitude1*/ + double latr2 = atan( (1.0-f)*tan(rlat2) ); /*Reduced latitude2*/ + double P = (latr1+latr2)/2; + double Q = (latr2-latr1)/2; + double sinP = sin(P); + double sin2P = sinP*sinP; + double cosQ = cos(Q); + double cos2Q = cosQ*cosQ; + double cosc2 = cos(c/2); + double cos2c2 = cosc2*cosc2; + + double sinQ = sin(Q); + double sin2Q = sinQ*sinQ; + double cosP = cos(P); + double cos2P = cosP*cosP; + double sinc2 = sin(c/2); + double sin2c2 = sinc2*sinc2; + + double X = (c - sin(c))* sin2P * cos2Q / cos2c2; + double Y = (c + sin(c))*sin2Q*cos2P/sin2c2; + double dist = major * (c - f*(X+Y)/2); + return dist; +} + int grib_nearest_find_multiple( const grib_handle* h, int is_lsm, const double* inlats, const double* inlons, long npoints, @@ -409,7 +453,7 @@ int grib_nearest_find_generic( /* Ignore latitudes too far from our point */ } else { - dist = grib_nearest_distance(radius, inlon, inlat, lon, lat); + dist = grib_nearest_distance_spherical(radius, inlon, inlat, lon, lat); if (dist < min_dist) min_dist = dist; /*printf("Candidate: lat=%.5f lon=%.5f dist=%f Idx=%ld Val=%f\n",lat,lon,dist,the_index,the_value);*/ diff --git a/src/grib_nearest_class_latlon_reduced.c b/src/grib_nearest_class_latlon_reduced.c index caa000c08..6f9346cb0 100644 --- a/src/grib_nearest_class_latlon_reduced.c +++ b/src/grib_nearest_class_latlon_reduced.c @@ -346,7 +346,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] = grib_nearest_distance(radius, inlon, inlat, + self->distances[kk] = grib_nearest_distance_spherical(radius, 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 7467ea0cc..6b52e188b 100644 --- a/src/grib_nearest_class_reduced.c +++ b/src/grib_nearest_class_reduced.c @@ -409,7 +409,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] = grib_nearest_distance(radius, inlon, inlat, + self->distances[kk] = grib_nearest_distance_spherical(radius, inlon, inlat, self->lons[self->k[kk]], self->lats[self->j[jj]]); kk++; } @@ -548,7 +548,7 @@ static int find(grib_nearest* nearest, grib_handle* h, kk = 0; for (ii = 0; ii < 2; ii++) { for (jj = 0; jj < 2; jj++) { - self->distances[kk] = grib_nearest_distance(radius, inlon, inlat, + self->distances[kk] = grib_nearest_distance_spherical(radius, 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 7e789cc8c..327bacafd 100644 --- a/src/grib_nearest_class_regular.c +++ b/src/grib_nearest_class_regular.c @@ -193,7 +193,7 @@ static int find(grib_nearest* nearest, grib_handle* h, for (ii=0;ii<2;ii++) { for (jj=0;jj<2;jj++) { self->k[kk]=self->i[ii]+self->lons_count*self->j[jj]-1; - self->distances[kk]=grib_nearest_distance(radius,inlon,inlat, + self->distances[kk]=grib_nearest_distance_spherical(radius,inlon,inlat, self->lons[self->i[ii]],self->lats[self->j[jj]]); kk++; } @@ -409,7 +409,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] = grib_nearest_distance(radius, inlon, inlat, + self->distances[kk] = grib_nearest_distance_spherical(radius, inlon, inlat, self->lons[self->i[ii]], self->lats[self->j[jj]]); kk++; }