Geo variogram

gstools/variogram/estimator.pyx

@@ -10,7 +10,7 @@ import numpy as np
 cimport cython
 from cython.parallel import prange, parallel
 from libcpp.vector cimport vector
-from libc.math cimport fabs, sqrt
+from libc.math cimport fabs, sqrt, pow, sin, cos, atan2, M_PI
 cimport numpy as np
 
 
@@ -36,6 +36,27 @@ cdef inline double _distance_2d(
 ) nogil:
     return sqrt((x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j]))
 
+cdef inline double _distance_2d_haversine(
+    const double[:] x,
+    const double[:] y,
+    const double[:] z,
+    const int i,
+    const int j
+) nogil:
+    # x is latitude in degree
+    # y is longitude in degree
+    cdef double deg_2_rad = M_PI / 180.0
+    cdef double diff_lat = (x[j] - x[i]) * deg_2_rad
+    cdef double diff_lon = (y[j] - y[i]) * deg_2_rad
+    cdef double arg = (
+        pow(sin(diff_lat/2.0), 2) +
+        cos(x[i]*deg_2_rad) *
+        cos(x[j]*deg_2_rad) *
+        pow(sin(diff_lon/2.0), 2)
+    )
+    # earth radius for WGS84 ellipsoid r ~ 6371km
+    return 12742000.0 * atan2(sqrt(arg), sqrt(1.0-arg))
+
 cdef inline double _distance_3d(
     const double[:] x,
     const double[:] y,
@@ -117,7 +138,8 @@ def unstructured(
     const double[:] x,
     const double[:] y=None,
     const double[:] z=None,
-    str estimator_type='m'
+    str estimator_type='m',
+    str distance_type='e'
 ):
     if x.shape[0] != f.shape[0]:
         raise ValueError('len(x) = {0} != len(f) = {1} '.
@@ -138,6 +160,8 @@ def unstructured(
             raise ValueError('len(y) = {0} != len(f) = {1} '.
                              format(y.shape[0], f.shape[0]))
         distance = _distance_2d
+        if distance_type == 's':
+            distance = _distance_2d_haversine
     # 1d
     else:
         distance = _distance_1d

gstools/variogram/variogram.py

@@ -32,6 +32,18 @@ def _set_estimator(estimator):
         )
     return cython_estimator
 
+def _set_distance(distance):
+    """Translate the verbose Python distance identifier to single char."""
+    if distance.lower() == "euclidean":
+        cython_distance = "e"
+    elif distance.lower() == "great-circle":
+        cython_distance = "s"
+    else:
+        raise ValueError(
+            "Unknown variogram distance function " + str(distance)
+        )
+    return cython_distance
+
 
 def vario_estimate_unstructured(
     pos,
@@ -40,6 +52,7 @@ def vario_estimate_unstructured(
     sampling_size=None,
     sampling_seed=None,
     estimator="matheron",
+    distance="euclidean",
 ):
     r"""
     Estimates the variogram on a unstructured grid.
@@ -92,6 +105,13 @@ def vario_estimate_unstructured(
             * "cressie": an estimator more robust to outliers
 
         Default: "matheron"
+    distance : :class:`str`, optional
+        the distance function, possible choices:
+
+            * "euclidean": the distance in Euclidean spaces
+            * "great-circle": the great-circle distance on spheres
+
+        Default: "euclidean"
 
     Returns
     -------
@@ -116,11 +136,18 @@ def vario_estimate_unstructured(
             z = z[sampled_idx]
 
     cython_estimator = _set_estimator(estimator)
+    cython_distance = _set_distance(distance)
 
     return (
         bin_centres,
         unstructured(
-            field, bin_edges, x, y, z, estimator_type=cython_estimator
+            field,
+            bin_edges,
+            x,
+            y,
+            z,
+            estimator_type=cython_estimator,
+            distance_type=cython_distance,
         ),
     )
 

tests/test_variogram_unstructured.py

@@ -183,6 +183,15 @@ def test_sampling_3d(self):
         self.assertAlmostEqual(gamma[len(gamma) // 2], var, places=2)
         self.assertAlmostEqual(gamma[-1], var, places=2)
 
+    def test_spherical_dist(self):
+        lat = np.array((38.725267, 0.313611, 31.205167))
+        lon = np.array((-9.150019, 32.581111, -7.862263))
+        d = np.array((2, 1200, 2570))
+        bins = np.array((0.0, 844999.0, 6021999.0, 10000000.0))
+        _, gamma = vario_estimate_unstructured((lat, lon), d, bins, distance='great-circle')
+        self.assertAlmostEqual(gamma[0], 3297312.)
+        self.assertAlmostEqual(gamma[1], 828026.)
+
     def test_assertions(self):
         x = np.arange(0, 10)
         x_e = np.arange(0, 11)