Implemented directional variograms for unstructured

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, atan2
 cimport numpy as np
 
 
@@ -47,6 +47,49 @@ cdef inline double _distance_3d(
                 (y[i] - y[j]) * (y[i] - y[j]) +
                 (z[i] - z[j]) * (z[i] - z[j]))
 
+cdef inline bint _angle_test_1d(
+    const double[:] x,
+    const double[:] y,
+    const double[:] z,
+    const double[:] angles,
+    const double angles_tol,
+    const int i,
+    const int j
+) nogil:
+    return True
+
+cdef inline bint _angle_test_2d(
+    const double[:] x,
+    const double[:] y,
+    const double[:] z,
+    const double[:] angles,
+    const double angles_tol,
+    const int i,
+    const int j
+) nogil:
+    cdef double dx = x[i] - x[j]
+    cdef double dy = y[i] - y[j]
+
+    cdef double phi = atan2(dy,dx)
+    return fabs(phi - angles[0]) <= angles_tol
+
+
+cdef inline bint _angle_test_3d(
+    const double[:] x,
+    const double[:] y,
+    const double[:] z,
+    const double[:] angles,
+    const double angles_tol,
+    const int i,
+    const int j
+) nogil:
+    cdef double dx = x[i] - x[j]
+    cdef double dy = y[i] - y[j]
+    cdef double dz = z[i] - z[j]
+
+    cdef double theta = atan2(dz,sqrt(dx + dy))
+    cdef double phi = atan2(dy,dx)
+    return fabs(theta - angles[0]) <= angles_tol and fabs(phi - angles[1]) <= angles_tol
 
 cdef inline double estimator_matheron(const double f_diff) nogil:
     return f_diff * f_diff
@@ -110,13 +153,25 @@ ctypedef double (*_dist_func)(
     const int
 ) nogil
 
+ctypedef bint (*_angle_test_func)(
+    const double[:],
+    const double[:],
+    const double[:],
+    const double[:],
+    const double,
+    const int,
+    const int
+) nogil
+
 
 def unstructured(
     const double[:] f,
     const double[:] bin_edges,
     const double[:] x,
     const double[:] y=None,
     const double[:] z=None,
+    const double[:] angles=None,
+    const double angles_tol=0.0872665,
     str estimator_type='m'
 ):
     if x.shape[0] != f.shape[0]:
@@ -126,21 +181,35 @@ def unstructured(
         raise ValueError('len(bin_edges) too small')
 
     cdef _dist_func distance
+    cdef _angle_test_func angle_test
+
     # 3d
     if z is not None:
         if z.shape[0] != f.shape[0]:
             raise ValueError('len(z) = {0} != len(f) = {1} '.
                              format(z.shape[0], f.shape[0]))
         distance = _distance_3d
+        angle_test = _angle_test_3d
     # 2d
     elif y is not None:
         if y.shape[0] != f.shape[0]:
             raise ValueError('len(y) = {0} != len(f) = {1} '.
                              format(y.shape[0], f.shape[0]))
         distance = _distance_2d
+        angle_test = _angle_test_2d
     # 1d
     else:
         distance = _distance_1d
+        angle_test = _angle_test_1d
+
+    if angles is not None:
+        if z is not None and angles.size < 2:
+            raise ValueError('3d requested but only one angle given')
+        if y is not None and angles.size < 1:
+            raise ValueError('2d with angle requested but no angle given')
+
+    if angles_tol <= 0:
+        raise ValueError('tolerance for angle search masks must be > 0')
 
     cdef _estimator_func estimator_func = choose_estimator_func(estimator_type)
     cdef _normalization_func normalization_func = (
@@ -160,8 +229,9 @@ def unstructured(
             for k in range(j+1, k_max):
                 dist = distance(x, y, z, k, j)
                 if dist >= bin_edges[i] and dist < bin_edges[i+1]:
-                    counts[i] += 1
-                    variogram[i] += estimator_func(f[k] - f[j])
+                    if angles is None or angle_test(x, y, z, angles, angles_tol, k, j):
+                        counts[i] += 1
+                        variogram[i] += estimator_func(f[k] - f[j])
 
     normalization_func(variogram, counts)
     return np.asarray(variogram)

gstools/variogram/variogram.py

@@ -37,6 +37,8 @@ def vario_estimate_unstructured(
     pos,
     field,
     bin_edges,
+    angles=None,
+    angles_tol=0.0872665,
     sampling_size=None,
     sampling_seed=None,
     estimator="matheron",
@@ -75,6 +77,17 @@ def vario_estimate_unstructured(
         the spatially distributed data
     bin_edges : :class:`numpy.ndarray`
         the bins on which the variogram will be calculated
+    angles : :class:`numpy.ndarray`
+        the angles of the main axis to calculate the variogram for in radians
+        angle definitions from ISO standard 80000-2:2009 
+        for 1d this parameter will have no effect at all
+        for 2d supply one angle which is azimuth φ (ccw from +x in xy plane)
+        for 3d supply two angles which are inclination θ (cw from +z)
+        and azimuth φ (ccw from +x in xy plane)
+    angles_tol : class:`float`
+        the tolerance around the variogram angle to count a point as being
+        within this direction from another point (the angular tolerance around
+        the directional vector given by angles)
     sampling_size : :class:`int` or :any:`None`, optional
         for large input data, this method can take a long
         time to compute the variogram, therefore this argument specifies
@@ -100,7 +113,18 @@ def vario_estimate_unstructured(
     field = np.array(field, ndmin=1, dtype=np.double)
     bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double)
     x, y, z, dim = pos2xyz(pos, calc_dim=True, dtype=np.double)
+
     bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0
+
+    if angles is not None:
+        angles = np.array(angles, ndmin=1, dtype=np.double)
+        if angles.size == 0:
+            angles = np.append(angles, [0,0,0])
+        elif angles.size == 1:
+            angles = np.append(angles, [0,0])
+        elif angles.size == 2:
+            angles = np.append(angles, [0])
+
 
     if sampling_size is not None and sampling_size < len(field):
         sampled_idx = np.random.RandomState(sampling_seed).choice(
@@ -118,7 +142,7 @@ def vario_estimate_unstructured(
     return (
         bin_centres,
         unstructured(
-            field, bin_edges, x, y, z, estimator_type=cython_estimator
+            field, bin_edges, x, y, z, angles, angles_tol, estimator_type=cython_estimator
         ),
     )
 

tests/test_variogram_unstructured.py

@@ -5,12 +5,38 @@
 
 import unittest
 import numpy as np
-from gstools import vario_estimate_unstructured
-
+from gstools import vario_estimate_unstructured, Exponential, SRF
+import gstools as gs
 
 class TestVariogramUnstructured(unittest.TestCase):
     def setUp(self):
-        pass
+
+        # this code generates the testdata for the 2D rotation test cases
+        x = np.random.RandomState(19970221).rand(20) * 10.0 - 5
+        y = np.zeros_like(x)
+        model = gs.Exponential(dim=2, var=2, len_scale=8)
+        srf = gs.SRF(model, mean=0, seed=19970221)
+        field = srf((x, y))
+
+        bins = np.arange(10)
+        bin_center, gamma = gs.vario_estimate_unstructured((x, ), field, bins)
+        idx = np.argsort(x)
+        self.test_data_rotation_1 = {'gamma': gamma, 'x': x[idx], 'field': field[idx], 'bins': bins, 'bin_center': bin_center}
+
+        # CODE ABOVE SHOULD GENERATE THIS DATA
+        # x = np.array([
+        # -4.86210059, -4.1984934 , -3.9246953 , -3.28490663, -2.16332379,
+        # -1.87553275, -1.74125124, -1.27224687, -1.20931578, -0.2413368 ,
+        #  0.03200921,  1.17099773,  1.53863105,  1.64478688,  2.75252136,
+        #  3.3556915 ,  3.89828775,  4.21485964,  4.5364357 ,  4.79236969]),
+        # field = np.array([
+        # -1.10318365, -0.53566629, -0.41789049, -1.06167529,  0.38449961,
+        # -0.36550477, -0.98905552, -0.19352766,  0.16264266,  0.26920833,
+        #  0.05379665,  0.71275006,  0.36651935,  0.17366865,  1.20022343,
+        #  0.79385446,  0.69456069,  1.0733393 ,  0.71191592,  0.71969766])
+        # gamma_exp = np.array([
+        # 0.14260989, 0.18301197, 0.25855841, 0.29990083, 0.67914526,
+        # 0.60136535, 0.92875492, 1.46910435, 1.10165104])
 
     def test_doubles(self):
         x = np.arange(1, 11, 1, dtype=np.double)
@@ -217,6 +243,121 @@ def test_assertions(self):
             ValueError, vario_estimate_unstructured, [x], field_e, bins
         )
 
+    def test_angles_2D_x2x(self):
+
+        x = self.test_data_rotation_1['x']
+        field = self.test_data_rotation_1['field']
+        gamma_exp = self.test_data_rotation_1['gamma']
+        bins = self.test_data_rotation_1['bins']
+        y = np.zeros_like(x)
+
+        # test case 1.)
+        #    all along x axis on x axis
+
+        bin_centres, gamma = vario_estimate_unstructured(
+            (x, y),
+            field,
+            bins,
+            angles=[0]
+        )
+
+        for i in range(gamma.size):
+            self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3)
+
+    def test_angles_2D_y2x(self):
+
+        x = self.test_data_rotation_1['x']
+        field = self.test_data_rotation_1['field']
+        gamma_exp = self.test_data_rotation_1['gamma']
+        bins = self.test_data_rotation_1['bins']
+        y = np.zeros_like(x)
 
+        # test case 2.)
+        #    all along y axis on y axis but calculation for x axis
+
+        bin_centres, gamma = vario_estimate_unstructured(
+            (y, x),
+            field,
+            bins,
+            angles=[0]
+        )
+
+        for i in range(gamma.size):
+            self.assertAlmostEqual(0, gamma[i], places=3)
+
+    def test_angles_2D_y2y(self):
+
+        x = self.test_data_rotation_1['x']
+        field = self.test_data_rotation_1['field']
+        gamma_exp = self.test_data_rotation_1['gamma']
+        bins = self.test_data_rotation_1['bins']
+        y = np.zeros_like(x)
+
+        # test case 3.)
+        #    all along y axis on y axis and calculation for y axis
+
+        bin_centres, gamma = vario_estimate_unstructured(
+            (y, x),
+            field,
+            bins,
+            angles=[np.pi/2.]
+        )
+
+        for i in range(gamma.size):
+            self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3)
+
+    def test_angles_2D_xy2x(self):
+
+        x = self.test_data_rotation_1['x']
+        field = self.test_data_rotation_1['field']
+        gamma_exp = self.test_data_rotation_1['gamma']
+        bins = self.test_data_rotation_1['bins']
+        y = np.zeros_like(x)
+
+
+        # test case 4.)
+        #    data along 45deg axis but calculation for x axis
+
+        ccos, csin = np.cos(np.pi/4.), np.sin(np.pi/4.)
+
+        xr = [xx * ccos - yy * csin for xx, yy in zip(x, y)]
+        yr = [xx * csin + yy * ccos for xx, yy in zip(x, y)]
+
+        bin_centres, gamma = vario_estimate_unstructured(
+            (xr, yr),
+            field,
+            bins,
+            angles=[0]
+        )
+
+        for i in range(gamma.size):
+            self.assertAlmostEqual(0, gamma[i], places=3)
+
+    def test_angles_2D_xy2xy(self):
+
+        x = self.test_data_rotation_1['x']
+        field = self.test_data_rotation_1['field']
+        gamma_exp = self.test_data_rotation_1['gamma']
+        bins = self.test_data_rotation_1['bins']
+        y = np.zeros_like(x)
+
+        # test case 5.)
+        #    data along 45deg axis and calculation for 45deg
+
+        ccos, csin = np.cos(np.pi/4.), np.sin(np.pi/4.)
+
+        xr = [xx * ccos - yy * csin for xx, yy in zip(x, y)]
+        yr = [xx * csin + yy * ccos for xx, yy in zip(x, y)]
+
+        bin_centres, gamma = vario_estimate_unstructured(
+            (xr, yr),
+            field,
+            bins,
+            angles=[np.pi/4.]
+        )
+
+        for i in range(gamma.size):
+            self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3)
+
 if __name__ == "__main__":
     unittest.main()