Geographic distance calculations (spherical and ellipsoid)

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,

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);*/

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

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

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