diff --git a/src/grib_api_prototypes.h b/src/grib_api_prototypes.h
index 165622bd8..362345ec6 100644
--- a/src/grib_api_prototypes.h
+++ b/src/grib_api_prototypes.h
@@ -1048,6 +1048,8 @@ int grib_get_gaussian_latitudes(long trunc, double* lats);
 int is_gaussian_global(double lat1, double lat2, double lon1, double lon2, long num_points_equator, const double* latitudes, double angular_precision);
 void rotate(const double inlat, const double inlon, const double angleOfRot, const double southPoleLat, const double southPoleLon, double* outlat, double* outlon);
 void unrotate(const double inlat, const double inlon, const double angleOfRot, const double southPoleLat, const double southPoleLon, double* outlat, double* outlon);
+double geographic_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2);
+double geographic_distance_ellipsoid(double major, double minor, double lon1, double lat1, double lon2, double lat2);
 
 /* grib_handle.c */
 grib_section* grib_section_create(grib_handle* h, grib_accessor* owner);
@@ -1376,8 +1378,6 @@ 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_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_geography.c b/src/grib_geography.c
index a7cb2e0e7..3c60f9f28 100644
--- a/src/grib_geography.c
+++ b/src/grib_geography.c
@@ -4140,3 +4140,73 @@ void unrotate(const double inlat, const double inlon,
     *outlat = ret_lat;
     *outlon = ret_lon;
 }
+
+#define RADIAN(x) ((x)*acos(0.0) / 90.0)
+
+/* radius is in km, angles in degrees */
+/* Spherical Law of Cosines */
+double geographic_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2)
+{
+    double rlat1 = RADIAN(lat1);
+    double rlat2 = RADIAN(lat2);
+    double rlon1 = lon1;
+    double rlon2 = lon2;
+    double a;
+
+    if (lat1 == lat2 && lon1 == lon2) {
+        return 0.0; /* the two points are identical */
+    }
+
+    if (rlon1 >= 360) rlon1 -= 360.0;
+    rlon1 = RADIAN(rlon1);
+    if (rlon2 >= 360) rlon2 -= 360.0;
+    rlon2 = RADIAN(rlon2);
+
+    a = sin(rlat1) * sin(rlat2) + cos(rlat1) * cos(rlat2) * cos(rlon2 - rlon1);
+    DebugAssert(a >= -1 && a <= 1);
+
+    return radius * acos(a);
+}
+
+/* major and minor axes in km, angles in degrees */
+double geographic_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;
+}
diff --git a/src/grib_nearest.c b/src/grib_nearest.c
index 91014d12f..99788473e 100644
--- a/src/grib_nearest.c
+++ b/src/grib_nearest.c
@@ -114,76 +114,6 @@ void grib_binary_search(double xx[], const unsigned long n, double x,
     }
 }
 
-#define RADIAN(x) ((x)*acos(0.0) / 90.0)
-
-/* radius is in km, angles in degrees */
-/* Spherical Law of Cosines */
-double grib_nearest_distance_spherical(double radius, double lon1, double lat1, double lon2, double lat2)
-{
-    double rlat1 = RADIAN(lat1);
-    double rlat2 = RADIAN(lat2);
-    double rlon1 = lon1;
-    double rlon2 = lon2;
-    double a;
-
-    if (lat1 == lat2 && lon1 == lon2) {
-        return 0.0; /* the two points are identical */
-    }
-
-    if (rlon1 >= 360) rlon1 -= 360.0;
-    rlon1 = RADIAN(rlon1);
-    if (rlon2 >= 360) rlon2 -= 360.0;
-    rlon2 = RADIAN(rlon2);
-
-    a = sin(rlat1) * sin(rlat2) + cos(rlat1) * cos(rlat2) * cos(rlon2 - rlon1);
-    DebugAssert(a >= -1 && a <= 1);
-
-    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,
@@ -447,7 +377,7 @@ int grib_nearest_find_generic(
                 /* Ignore latitudes too far from our point */
             }
             else {
-                double dist = grib_nearest_distance_spherical(radius, inlon, inlat, lon, lat);
+                double dist = geographic_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 6f9346cb0..3be537bea 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_spherical(radius, inlon, inlat,
+                self->distances[kk] = geographic_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 6b52e188b..6b5c0b119 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_spherical(radius, inlon, inlat,
+                self->distances[kk] = geographic_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_spherical(radius, inlon, inlat,
+                self->distances[kk] = geographic_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 327bacafd..72a13d426 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_spherical(radius,inlon,inlat,
+                self->distances[kk]=geographic_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_spherical(radius, inlon, inlat,
+                self->distances[kk] = geographic_distance_spherical(radius, inlon, inlat,
                                                             self->lons[self->i[ii]], self->lats[self->j[jj]]);
                 kk++;
             }