From 80746e23e786a8f472f3192c4ed9155e1e6d44d3 Mon Sep 17 00:00:00 2001
From: TobiasGlaubach <tobias.glaubach@gmail.com>
Date: Mon, 27 Apr 2020 23:00:11 +0200
Subject: [PATCH 01/15] implemented a prototype for the unstructured function
 with angles.

---
 gstools/variogram/estimator.pyx | 76 +++++++++++++++++++++++++++++++--
 1 file changed, 73 insertions(+), 3 deletions(-)

diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx
index 8762b946..75a1c3fa 100644
--- a/gstools/variogram/estimator.pyx
+++ b/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,6 +153,16 @@ 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,
@@ -117,6 +170,8 @@ def unstructured(
     const double[:] x,
     const double[:] y=None,
     const double[:] z=None,
+    const double[:] angles=None,
+    const double angles_tol_rad=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)

From d99e33be27046dc9c7872d82b576eb4e94cde198 Mon Sep 17 00:00:00 2001
From: Tobi Glaubach <tobias.glaubach@gmail.com>
Date: Sat, 2 May 2020 19:29:26 +0200
Subject: [PATCH 02/15] minor fix with variable name

---
 gstools/variogram/estimator.pyx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx
index 75a1c3fa..2dc7e8d0 100644
--- a/gstools/variogram/estimator.pyx
+++ b/gstools/variogram/estimator.pyx
@@ -171,7 +171,7 @@ def unstructured(
     const double[:] y=None,
     const double[:] z=None,
     const double[:] angles=None,
-    const double angles_tol_rad=0.0872665,
+    const double angles_tol=0.0872665,
     str estimator_type='m'
 ):
     if x.shape[0] != f.shape[0]:

From 8973b18098c5f01ec33e81534f4fc30f72944e23 Mon Sep 17 00:00:00 2001
From: Tobi Glaubach <tobias.glaubach@gmail.com>
Date: Sat, 2 May 2020 19:30:34 +0200
Subject: [PATCH 03/15] added docstring and arguments for angle estimation

---
 gstools/variogram/variogram.py | 26 +++++++++++++++++++++++++-
 1 file changed, 25 insertions(+), 1 deletion(-)

diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py
index ab87562f..f915e8b4 100644
--- a/gstools/variogram/variogram.py
+++ b/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 : :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
         ),
     )
 

From b7355f845c1e029e4d6d6e8052c1dd1c90cb3306 Mon Sep 17 00:00:00 2001
From: Tobi Glaubach <tobias.glaubach@gmail.com>
Date: Sat, 2 May 2020 19:31:00 +0200
Subject: [PATCH 04/15] added a working version which tests the basic
 functionality

---
 tests/test_variogram_unstructured.py | 95 ++++++++++++++++++++++++++++
 1 file changed, 95 insertions(+)

diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py
index c5a59a85..b1e6b5f0 100644
--- a/tests/test_variogram_unstructured.py
+++ b/tests/test_variogram_unstructured.py
@@ -217,6 +217,101 @@ def test_assertions(self):
             ValueError, vario_estimate_unstructured, [x], field_e, bins
         )
 
+    def test_angles(self):
+        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])
+        
+        bins = np.arange(0, 10)
+        
+        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)
+
 
+        # 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)
+            
+        # 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)
+            
+        # 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)
+            
+        # 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()

From 8cbdc2b3713b8406734041dda9252a38509d9752 Mon Sep 17 00:00:00 2001
From: Tobi Glaubach <tobias.glaubach@gmail.com>
Date: Sat, 9 May 2020 12:11:01 +0200
Subject: [PATCH 05/15] docstring adaption

---
 gstools/variogram/variogram.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py
index f915e8b4..99c8c1ea 100644
--- a/gstools/variogram/variogram.py
+++ b/gstools/variogram/variogram.py
@@ -84,7 +84,7 @@ def vario_estimate_unstructured(
         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 : :float
+    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)

From 4f64b85a7002bec1c8b73911058383ef7ac93d30 Mon Sep 17 00:00:00 2001
From: Tobi Glaubach <tobias.glaubach@gmail.com>
Date: Sat, 9 May 2020 12:11:31 +0200
Subject: [PATCH 06/15] changed for automatic testcase data generation and
 split up the test cases

---
 tests/test_variogram_unstructured.py | 90 +++++++++++++++++++++-------
 1 file changed, 68 insertions(+), 22 deletions(-)

diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py
index b1e6b5f0..e31471e9 100644
--- a/tests/test_variogram_unstructured.py
+++ b/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,24 +243,12 @@ def test_assertions(self):
             ValueError, vario_estimate_unstructured, [x], field_e, bins
         )
 
-    def test_angles(self):
-        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])
-        
-        bins = np.arange(0, 10)
-        
+    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.)
@@ -250,6 +264,13 @@ def test_angles(self):
         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
@@ -264,6 +285,14 @@ def test_angles(self):
         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
             
@@ -276,7 +305,16 @@ def test_angles(self):
         
         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
             
@@ -295,6 +333,14 @@ def test_angles(self):
         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
             

From 16b517f351409493ab17252d3a11f470e59e0c0e Mon Sep 17 00:00:00 2001
From: Tobi Glaubach <tobias.glaubach@gmail.com>
Date: Thu, 14 May 2020 10:40:23 +0200
Subject: [PATCH 07/15] changed default angle tolerance to 25deg

---
 gstools/variogram/estimator.pyx | 2 +-
 gstools/variogram/variogram.py  | 5 +++--
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx
index 2dc7e8d0..d826b73c 100644
--- a/gstools/variogram/estimator.pyx
+++ b/gstools/variogram/estimator.pyx
@@ -171,7 +171,7 @@ def unstructured(
     const double[:] y=None,
     const double[:] z=None,
     const double[:] angles=None,
-    const double angles_tol=0.0872665,
+    const double angles_tol=0.436332,
     str estimator_type='m'
 ):
     if x.shape[0] != f.shape[0]:
diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py
index 99c8c1ea..2702d27b 100644
--- a/gstools/variogram/variogram.py
+++ b/gstools/variogram/variogram.py
@@ -38,7 +38,7 @@ def vario_estimate_unstructured(
     field,
     bin_edges,
     angles=None,
-    angles_tol=0.0872665,
+    angles_tol=0.436332,
     sampling_size=None,
     sampling_seed=None,
     estimator="matheron",
@@ -87,7 +87,8 @@ def vario_estimate_unstructured(
     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)
+        the directional vector given by angles) 
+        Default: 25°≈0.436332
     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

From 360151a68f096fcb54d58956d075aeba86a9e6a2 Mon Sep 17 00:00:00 2001
From: Tobi Glaubach <tobias.glaubach@gmail.com>
Date: Tue, 2 Jun 2020 06:55:58 +0200
Subject: [PATCH 08/15] added option to also return the counts (number of
 pairs) from unstructured

---
 gstools/variogram/estimator.pyx |  2 +-
 gstools/variogram/variogram.py  | 29 ++++++++++++++++++++++-------
 2 files changed, 23 insertions(+), 8 deletions(-)

diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx
index d826b73c..fbcd06bd 100644
--- a/gstools/variogram/estimator.pyx
+++ b/gstools/variogram/estimator.pyx
@@ -234,7 +234,7 @@ def unstructured(
                         variogram[i] += estimator_func(f[k] - f[j])
 
     normalization_func(variogram, counts)
-    return np.asarray(variogram)
+    return np.asarray(variogram), np.asarray(counts)
 
 
 def structured(const double[:,:,:] f, str estimator_type='m'):
diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py
index 2702d27b..c25dc6e2 100644
--- a/gstools/variogram/variogram.py
+++ b/gstools/variogram/variogram.py
@@ -42,6 +42,7 @@ def vario_estimate_unstructured(
     sampling_size=None,
     sampling_seed=None,
     estimator="matheron",
+    return_counts=False
 ):
     r"""
     Estimates the variogram on a unstructured grid.
@@ -104,11 +105,16 @@ def vario_estimate_unstructured(
             * "cressie": an estimator more robust to outliers
 
         Default: "matheron"
-
+    return_counts: class:`bool`, optional
+        if set to true, this function will also return the number of data 
+        points found at each lag distance as a third return value
+        Default: False
     Returns
     -------
     :class:`tuple` of :class:`numpy.ndarray`
-        the estimated variogram and the bin centers
+        1. the bin centers
+        2. the estimated variogram values at bin centers
+        (optional) 3. the number of points found at each bin center (see argument return_counts)
     """
     # TODO check_mesh
     field = np.array(field, ndmin=1, dtype=np.double)
@@ -140,12 +146,21 @@ def vario_estimate_unstructured(
 
     cython_estimator = _set_estimator(estimator)
 
-    return (
-        bin_centres,
-        unstructured(
+    estimates, counts = unstructured(
             field, bin_edges, x, y, z, angles, angles_tol, estimator_type=cython_estimator
-        ),
-    )
+        )
+
+    if return_counts:
+        return (
+            bin_centres,
+            estimates,
+            counts
+        )
+    else:
+        return (
+            bin_centres,
+            estimates
+        )
 
 
 def vario_estimate_structured(field, direction="x", estimator="matheron"):

From 86b2fdcee23d4505172c02dd46b6aa4a54a3a4da Mon Sep 17 00:00:00 2001
From: Tobi Glaubach <tobias.glaubach@gmail.com>
Date: Mon, 8 Jun 2020 15:22:01 +0200
Subject: [PATCH 09/15] implemented a meaningful test for 2d variogram
 estimation

---
 tests/test_variogram_unstructured.py | 68 ++++++++++++++++++----------
 1 file changed, 45 insertions(+), 23 deletions(-)

diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py
index e31471e9..1baef635 100644
--- a/tests/test_variogram_unstructured.py
+++ b/tests/test_variogram_unstructured.py
@@ -335,29 +335,51 @@ def test_angles_2D_xy2x(self):
             
     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)
+        seed = gs.random.MasterRNG(19970221)
+        rng = np.random.RandomState(seed())
+        rng = np.random
+        x = rng.randint(0, 100, size=3000)
+        y = rng.randint(0, 100, size=3000)
+
+        model = gs.Exponential(dim=2, var=1, len_scale=[12, 3], angles=np.pi / 8)
+        model_maj = gs.Exponential(dim=1, var=1, len_scale=[12], angles=np.pi / 8)
+        model_min = gs.Exponential(dim=1, var=1, len_scale=[3], angles=np.pi / 8 + np.pi/2)
+
+        srf = gs.SRF(model, seed=20170519)
+        field = srf((x, y))
+
+
+        bins = np.arange(0, 50, 2.5)
+        angle_mask = 22.5
+        angle_tol = 22.5
+
+        bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured(
+                (x, y),
+                field,
+                bins,
+                angles=[np.deg2rad(angle_mask)],
+                angles_tol=np.deg2rad(angle_tol)
+            )
+
+        bin_centers_min, gamma_min = gs.vario_estimate_unstructured(
+                (x, y),
+                field,
+                bins,
+                angles=[np.deg2rad(angle_mask + 90.0)],
+                angles_tol=np.deg2rad(angle_tol)
+            )
+
+        gamma_maj_real = model_maj.variogram(bin_centers_maj)
+        gamma_min_real = model_min.variogram(bin_centers_min)
+
+        # we have no real way of testing values, but we can test some basic properties which definitelly need to be true
+        # test that the major estimate aligns better with major real than minor real
+        self.assertTrue(np.sum((gamma_maj_real - gamma_maj)**2) < np.sum((gamma_min_real - gamma_maj)**2))
+        # test that the minor estimate aligns better with minor real than major real
+        self.assertTrue(np.sum((gamma_min_real - gamma_min)**2) < np.sum((gamma_maj_real - gamma_min)**2))
+        # test that both variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance
+        self.assertTrue((np.mean(gamma_min[-5:]) - np.var(field)) / np.var(field) < 0.1)
+        self.assertTrue((np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1)
 
-        # 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()

From 980946491ab08548cbd2356deead1209655d2f3e Mon Sep 17 00:00:00 2001
From: Tobi Glaubach <tobias.glaubach@gmail.com>
Date: Tue, 9 Jun 2020 14:34:21 +0200
Subject: [PATCH 10/15] =?UTF-8?q?bugfix=20for=203d=20case=20when=20elevati?=
 =?UTF-8?q?on=20is=2090=C2=B0=20or=20270=C2=B0?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 gstools/variogram/estimator.pyx | 32 +++++++++++++++++++++++++++++---
 1 file changed, 29 insertions(+), 3 deletions(-)

diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx
index fbcd06bd..e4845430 100644
--- a/gstools/variogram/estimator.pyx
+++ b/gstools/variogram/estimator.pyx
@@ -70,9 +70,26 @@ cdef inline bint _angle_test_2d(
     cdef double dx = x[i] - x[j]
     cdef double dy = y[i] - y[j]
 
+    # azimuth
     cdef double phi = atan2(dy,dx)
     return fabs(phi - angles[0]) <= angles_tol
 
+cdef inline bint _angle_test_3d_ignore_azim(
+    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]
+
+    # elevation
+    cdef double theta = atan2(dz,sqrt(dx + dy))
+    return fabs(angles[1] - theta) <= angles_tol
 
 cdef inline bint _angle_test_3d(
     const double[:] x,
@@ -87,9 +104,12 @@ cdef inline bint _angle_test_3d(
     cdef double dy = y[i] - y[j]
     cdef double dz = z[i] - z[j]
 
-    cdef double theta = atan2(dz,sqrt(dx + dy))
+    # azimuth
     cdef double phi = atan2(dy,dx)
-    return fabs(theta - angles[0]) <= angles_tol and fabs(phi - angles[1]) <= angles_tol
+    # elevation
+    cdef double theta = atan2(dz,sqrt(dx + dy))
+
+    return sqrt((angles[0] - phi)**2 + (angles[1] - theta)**2) <= angles_tol
 
 cdef inline double estimator_matheron(const double f_diff) nogil:
     return f_diff * f_diff
@@ -189,7 +209,13 @@ def unstructured(
             raise ValueError('len(z) = {0} != len(f) = {1} '.
                              format(z.shape[0], f.shape[0]))
         distance = _distance_3d
-        angle_test = _angle_test_3d
+
+        # check if elevation is perpendicular to original axis
+        if angles is not None and fabs(angles[1] - 0.5*np.pi) < 1e-12 or fabs(angles[1] - 1.5*np.pi) < 1e-12:
+            # if so it does not matter how we rotated before (at least when dealing with two points)
+            angle_test = _angle_test_3d_ignore_azim
+        else:
+            angle_test = _angle_test_3d
     # 2d
     elif y is not None:
         if y.shape[0] != f.shape[0]:

From 60ca648904f8ffad62545c3dd9d28184a472153a Mon Sep 17 00:00:00 2001
From: Tobi Glaubach <tobias.glaubach@gmail.com>
Date: Tue, 9 Jun 2020 14:34:43 +0200
Subject: [PATCH 11/15] implemented some basic 3d test cases

---
 tests/test_variogram_unstructured.py | 133 ++++++++++++++++++++++++++-
 1 file changed, 130 insertions(+), 3 deletions(-)

diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py
index 1baef635..c219ec22 100644
--- a/tests/test_variogram_unstructured.py
+++ b/tests/test_variogram_unstructured.py
@@ -333,7 +333,7 @@ def test_angles_2D_xy2x(self):
         for i in range(gamma.size):
             self.assertAlmostEqual(0, gamma[i], places=3)
             
-    def test_angles_2D_xy2xy(self):
+    def test_angles_2D_estim(self):
 
         seed = gs.random.MasterRNG(19970221)
         rng = np.random.RandomState(seed())
@@ -342,8 +342,8 @@ def test_angles_2D_xy2xy(self):
         y = rng.randint(0, 100, size=3000)
 
         model = gs.Exponential(dim=2, var=1, len_scale=[12, 3], angles=np.pi / 8)
-        model_maj = gs.Exponential(dim=1, var=1, len_scale=[12], angles=np.pi / 8)
-        model_min = gs.Exponential(dim=1, var=1, len_scale=[3], angles=np.pi / 8 + np.pi/2)
+        model_maj = gs.Exponential(dim=1, var=1, len_scale=[12])
+        model_min = gs.Exponential(dim=1, var=1, len_scale=[3])
 
         srf = gs.SRF(model, seed=20170519)
         field = srf((x, y))
@@ -381,5 +381,132 @@ def test_angles_2D_xy2xy(self):
         self.assertTrue((np.mean(gamma_min[-5:]) - np.var(field)) / np.var(field) < 0.1)
         self.assertTrue((np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1)
 
+    def test_angles_line_3D(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']            
+
+        def test_xyz(x, y, z, angles, gamma_exp):
+            bin_centres, gamma = vario_estimate_unstructured(
+                (x, y, z),
+                field,
+                bins,
+                angles=angles
+            )
+
+            if np.ndim(gamma_exp) == 0:
+                gamma_exp = np.ones_like(gamma) * gamma_exp
+
+            for i in range(gamma.size):
+                self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3)
+
+        # all along x axis and calculation for x axis
+        test_xyz(x, np.zeros_like(x), np.zeros_like(x), [0], gamma_exp)
+        # all along y axis and calculation for x axis
+        test_xyz(np.zeros_like(x), x, np.zeros_like(x), [0], 0)
+        # all along z axis and calculation for x axis
+        test_xyz(np.zeros_like(x), np.zeros_like(x), x, [0], 0)
+
+        angles_rot_azim_90 = [np.pi/2]
+        # all along x axis and calculation for y axis
+        test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90, 0)
+        # all along y axis and calculation for y axis
+        test_xyz(np.zeros_like(x), x, np.zeros_like(x), angles_rot_azim_90, gamma_exp)
+        # all along z axis and calculation for y axis
+        test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90, 0)
+
+        # for elevation it is important to check, that IF elevation is 90° or 270° it does 
+        # not matter how we rotated before, since any rotation around z (in XY plane) 
+        # followed by a rotation around x' (in YZ' plane) by 90° will result in the same
+        # coordinates, (when the structure is two points with zero extend)
+
+        # test with [0, 90]
+        angles_rot_azim_90_elev_90 = [0, np.pi/2]
+        # all along x axis and calculation for z axis
+        test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90_elev_90, 0)
+        # all along y axis and calculation for z axis
+        test_xyz(np.zeros_like(x), x, np.zeros_like(x), angles_rot_azim_90_elev_90, 0)
+        # all along z axis and calculation for z axis
+        test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90_elev_90, gamma_exp)
+
+        # test with [90, 90]     
+        angles_rot_azim_90_elev_90 = [np.pi/2, np.pi/2]
+        # all along x axis and calculation for z axis
+        test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90_elev_90, 0)
+        # all along y axis and calculation for z axis
+        test_xyz(np.zeros_like(x), x, np.zeros_like(x), angles_rot_azim_90_elev_90, 0)
+        # all along z axis and calculation for z axis
+        test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90_elev_90, gamma_exp)
+
+    def test_angles_3D_estim(self):
+
+        seed = gs.random.MasterRNG(19970221)
+        rng = np.random.RandomState(seed())
+        rng = np.random
+        x = rng.randint(0, 50, size=1000)
+        y = rng.randint(0, 50, size=1000)
+        z = rng.randint(0, 50, size=1000)
+
+        model = gs.Exponential(dim=3, var=1, len_scale=[10, 5, 2], angles=[np.pi / 8, np.pi / 16])
+        model_maj = gs.Exponential(dim=1, var=1, len_scale=[10])
+        model_min1 = gs.Exponential(dim=1, var=1, len_scale=[5])
+        model_min2 = gs.Exponential(dim=1, var=1, len_scale=[2])
+
+        srf = gs.SRF(model, seed=20170519)
+        field = srf((x, y, z))
+
+        bins = np.arange(0, 25, 5)
+        angle_mask = [np.pi / 8, np.pi / 16]
+        angle_tol = 22.5
+
+        # in x'
+        bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured(
+                (x, y, z),
+                field,
+                bins,
+                angles=[angle_mask[0], angle_mask[1]],
+                angles_tol=np.deg2rad(angle_tol)
+            )
+
+        # in y'
+        bin_centers_min1, gamma_min1 = gs.vario_estimate_unstructured(
+                (x, y, z),
+                field,
+                bins,
+                angles=[angle_mask[0] + 0.5*np.pi, angle_mask[1]],
+                angles_tol=np.deg2rad(angle_tol)
+            )
+
+        # in z'
+        bin_centers_min2, gamma_min2 = gs.vario_estimate_unstructured(
+                (x, y, z),
+                field,
+                bins,
+                angles=[angle_mask[0], angle_mask[1] + 0.5*np.pi],
+                angles_tol=np.deg2rad(angle_tol)
+            )
+
+        gamma_maj_real = model_maj.variogram(bin_centers_maj)
+        gamma_min1_real = model_min1.variogram(bin_centers_min1)
+        gamma_min2_real = model_min2.variogram(bin_centers_min2)
+
+        # we have no real way of testing values, but we can test some basic properties which definitelly need to be true
+        # test that the major estimate aligns better with major real than the minor reals
+        self.assertTrue(np.sum((gamma_maj_real - gamma_maj)**2) < np.sum((gamma_min1_real - gamma_maj)**2))
+        self.assertTrue(np.sum((gamma_maj_real - gamma_maj)**2) < np.sum((gamma_min2_real - gamma_maj)**2))
+        # test that the minor1 estimate aligns better with minor1 real than major real and minor2 real
+        self.assertTrue(np.sum((gamma_min1_real - gamma_min1)**2) < np.sum((gamma_maj_real - gamma_min1)**2))
+        self.assertTrue(np.sum((gamma_min1_real - gamma_min1)**2) < np.sum((gamma_min2_real - gamma_min1)**2))
+        # test that the minor2 estimate aligns better with minor2 real than major real and minor1 real
+        self.assertTrue(np.sum((gamma_min2_real - gamma_min2)**2) < np.sum((gamma_maj_real - gamma_min2)**2))
+        self.assertTrue(np.sum((gamma_min2_real - gamma_min2)**2) < np.sum((gamma_min1_real - gamma_min2)**2))
+        # test that all variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance
+        self.assertTrue((np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1)
+        self.assertTrue((np.mean(gamma_min1[-5:]) - np.var(field)) / np.var(field) < 0.1)
+        self.assertTrue((np.mean(gamma_min2[-5:]) - np.var(field)) / np.var(field) < 0.1)
+        
+
 if __name__ == "__main__":
     unittest.main()

From b05569feea0cdeaa1140259614ab6513d2c95e74 Mon Sep 17 00:00:00 2001
From: MuellerSeb <mueller.seb@posteo.de>
Date: Thu, 5 Nov 2020 23:07:08 +0100
Subject: [PATCH 12/15] vario: cleanup cython routines; use greate-circle for
 tolerance in 3D; check both directions between point pairs

---
 gstools/variogram/estimator.pyx | 59 +++++++++++++++------------------
 1 file changed, 26 insertions(+), 33 deletions(-)

diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx
index e4845430..85781a54 100644
--- a/gstools/variogram/estimator.pyx
+++ b/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, atan2
+from libc.math cimport fabs, sqrt, atan2, acos, asin, sin, cos
 cimport numpy as np
 
 
@@ -69,27 +69,13 @@ cdef inline bint _angle_test_2d(
 ) nogil:
     cdef double dx = x[i] - x[j]
     cdef double dy = y[i] - y[j]
-
     # azimuth
-    cdef double phi = atan2(dy,dx)
-    return fabs(phi - angles[0]) <= angles_tol
-
-cdef inline bint _angle_test_3d_ignore_azim(
-    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]
-
-    # elevation
-    cdef double theta = atan2(dz,sqrt(dx + dy))
-    return fabs(angles[1] - theta) <= angles_tol
+    cdef double phi1 = atan2(dy,dx)
+    cdef double phi2 = atan2(-dy,-dx)
+    # check both directions (+/-)
+    cdef bint dir1 = fabs(phi1 - angles[0]) <= angles_tol
+    cdef bint dir2 = fabs(phi2 - angles[0]) <= angles_tol
+    return dir1 or dir2
 
 cdef inline bint _angle_test_3d(
     const double[:] x,
@@ -103,13 +89,26 @@ cdef inline bint _angle_test_3d(
     cdef double dx = x[i] - x[j]
     cdef double dy = y[i] - y[j]
     cdef double dz = z[i] - z[j]
-
+    cdef double dr = sqrt(dx**2 + dy**2 + dz**2)
     # azimuth
-    cdef double phi = atan2(dy,dx)
+    cdef double phi1 = atan2(dy, dx)
+    cdef double phi2 = atan2(-dy, -dx)
     # elevation
-    cdef double theta = atan2(dz,sqrt(dx + dy))
-
-    return sqrt((angles[0] - phi)**2 + (angles[1] - theta)**2) <= angles_tol
+    cdef double theta1 = acos(dz / dr)
+    cdef double theta2 = acos(-dz / dr)
+    # definitions for great-circle distance (for tolerance check)
+    cdef double dx1 = sin(theta1) * cos(phi1) - sin(angles[1]) * cos(angles[0])
+    cdef double dy1 = sin(theta1) * sin(phi1) - sin(angles[1]) * sin(angles[0])
+    cdef double dz1 = cos(theta1) - cos(angles[1])
+    cdef double dx2 = sin(theta2) * cos(phi2) - sin(angles[1]) * cos(angles[0])
+    cdef double dy2 = sin(theta2) * sin(phi2) - sin(angles[1]) * sin(angles[0])
+    cdef double dz2 = cos(theta2) - cos(angles[1])
+    cdef double dist1 = 2.0 * asin(sqrt(dx1**2 + dy1**2 + dz1**2) * 0.5)
+    cdef double dist2 = 2.0 * asin(sqrt(dx2**2 + dy2**2 + dz2**2) * 0.5)
+    # check both directions (+/-)
+    cdef bint dir1 = dist1 <= angles_tol
+    cdef bint dir2 = dist2 <= angles_tol
+    return dir1 or dir2
 
 cdef inline double estimator_matheron(const double f_diff) nogil:
     return f_diff * f_diff
@@ -209,13 +208,7 @@ def unstructured(
             raise ValueError('len(z) = {0} != len(f) = {1} '.
                              format(z.shape[0], f.shape[0]))
         distance = _distance_3d
-
-        # check if elevation is perpendicular to original axis
-        if angles is not None and fabs(angles[1] - 0.5*np.pi) < 1e-12 or fabs(angles[1] - 1.5*np.pi) < 1e-12:
-            # if so it does not matter how we rotated before (at least when dealing with two points)
-            angle_test = _angle_test_3d_ignore_azim
-        else:
-            angle_test = _angle_test_3d
+        angle_test = _angle_test_3d
     # 2d
     elif y is not None:
         if y.shape[0] != f.shape[0]:

From 536858c8fa4055d77907ff08de1a74bc81a3c1c4 Mon Sep 17 00:00:00 2001
From: MuellerSeb <mueller.seb@posteo.de>
Date: Thu, 5 Nov 2020 23:08:35 +0100
Subject: [PATCH 13/15] vario: doc update; correct intervals for angles;
 general formatting of angles array

---
 gstools/variogram/variogram.py | 98 +++++++++++++++++++---------------
 1 file changed, 54 insertions(+), 44 deletions(-)

diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py
index c25dc6e2..9ce2a499 100644
--- a/gstools/variogram/variogram.py
+++ b/gstools/variogram/variogram.py
@@ -42,7 +42,7 @@ def vario_estimate_unstructured(
     sampling_size=None,
     sampling_seed=None,
     estimator="matheron",
-    return_counts=False
+    return_counts=False,
 ):
     r"""
     Estimates the variogram on a unstructured grid.
@@ -53,7 +53,8 @@ def vario_estimate_unstructured(
        \gamma(r_k) = \frac{1}{2 N(r_k)} \sum_{i=1}^{N(r_k)} (z(\mathbf x_i) -
        z(\mathbf x_i'))^2 \; ,
 
-    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}` being the bins.
+    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
+    being the bins.
 
     Or if the estimator "cressie" was chosen:
 
@@ -62,13 +63,10 @@ def vario_estimate_unstructured(
        \left|z(\mathbf x_i) - z(\mathbf x_i')\right|^{0.5}\right)^4}
        {0.457 + 0.494 / N(r_k) + 0.045 / N^2(r_k)} \; ,
 
-    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}` being the bins.
+    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
+    being the bins.
     The Cressie estimator is more robust to outliers.
 
-    Notes
-    -----
-    Internally uses double precision and also returns doubles.
-
     Parameters
     ----------
     pos : :class:`list`
@@ -80,15 +78,15 @@ def vario_estimate_unstructured(
         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 
+        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)
+        for 3d supply two angles which are azimuth φ (ccw from +x in xy plane)
+        and inclination θ (cw from +z)
     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) 
+        the directional vector given by angles)
         Default: 25°≈0.436332
     sampling_size : :class:`int` or :any:`None`, optional
         for large input data, this method can take a long
@@ -106,32 +104,42 @@ def vario_estimate_unstructured(
 
         Default: "matheron"
     return_counts: class:`bool`, optional
-        if set to true, this function will also return the number of data 
+        if set to true, this function will also return the number of data
         points found at each lag distance as a third return value
         Default: False
+
     Returns
     -------
     :class:`tuple` of :class:`numpy.ndarray`
         1. the bin centers
         2. the estimated variogram values at bin centers
-        (optional) 3. the number of points found at each bin center (see argument return_counts)
+        3. (optional) the number of points found at each bin center
+           (see argument return_counts)
+
+    Notes
+    -----
+    Internally uses double precision and also returns doubles.
     """
     # TODO check_mesh
     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])
-            
+        # there are at most (dim-1) angles
+        angles = np.array(angles, ndmin=1, dtype=np.double).ravel()[
+            : (dim - 1)
+        ]
+        angles = np.pad(
+            angles, (0, dim - angles.size - 1), "constant", constant_values=0.0
+        )
+        # correct intervalls for angles (prepared for n-d)
+        if dim >= 2:
+            angles[0] = angles[0] % (2 * np.pi)  # azimuth in [0, 2pi)
+        if dim >= 3:
+            angles[1:] = angles[1:] % np.pi  # inclinations in [0, pi)
 
     if sampling_size is not None and sampling_size < len(field):
         sampled_idx = np.random.RandomState(sampling_seed).choice(
@@ -147,20 +155,20 @@ def vario_estimate_unstructured(
     cython_estimator = _set_estimator(estimator)
 
     estimates, counts = unstructured(
-            field, bin_edges, x, y, z, angles, angles_tol, estimator_type=cython_estimator
-        )
+        field,
+        bin_edges,
+        x,
+        y,
+        z,
+        angles,
+        angles_tol,
+        estimator_type=cython_estimator,
+    )
 
     if return_counts:
-        return (
-            bin_centres,
-            estimates,
-            counts
-        )
+        return bin_centres, estimates, counts
     else:
-        return (
-            bin_centres,
-            estimates
-        )
+        return bin_centres, estimates
 
 
 def vario_estimate_structured(field, direction="x", estimator="matheron"):
@@ -173,7 +181,8 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"):
        \gamma(r_k) = \frac{1}{2 N(r_k)} \sum_{i=1}^{N(r_k)} (z(\mathbf x_i) -
        z(\mathbf x_i'))^2 \; ,
 
-    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}` being the bins.
+    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
+    being the bins.
 
     Or if the estimator "cressie" was chosen:
 
@@ -182,17 +191,10 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"):
        \left|z(\mathbf x_i) - z(\mathbf x_i')\right|^{0.5}\right)^4}
        {0.457 + 0.494 / N(r_k) + 0.045 / N^2(r_k)} \; ,
 
-    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}` being the bins.
+    with :math:`r_k \leq \| \mathbf x_i - \mathbf x_i' \| < r_{k+1}`
+    being the bins.
     The Cressie estimator is more robust to outliers.
 
-    Warnings
-    --------
-    It is assumed that the field is defined on an equidistant Cartesian grid.
-
-    Notes
-    -----
-    Internally uses double precision and also returns doubles.
-
     Parameters
     ----------
     field : :class:`numpy.ndarray`
@@ -211,6 +213,14 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"):
     -------
     :class:`numpy.ndarray`
         the estimated variogram along the given direction.
+
+    Warnings
+    --------
+    It is assumed that the field is defined on an equidistant Cartesian grid.
+
+    Notes
+    -----
+    Internally uses double precision and also returns doubles.
     """
     try:
         mask = np.array(field.mask, dtype=np.int32)

From a8ef9f8c9a0785b7d29e3174d07b121c64d41f96 Mon Sep 17 00:00:00 2001
From: MuellerSeb <mueller.seb@posteo.de>
Date: Fri, 6 Nov 2020 01:40:05 +0100
Subject: [PATCH 14/15] vario: better handling of angle ranges

---
 gstools/variogram/estimator.pyx | 10 +++++-----
 gstools/variogram/variogram.py  | 18 ++++++++----------
 2 files changed, 13 insertions(+), 15 deletions(-)

diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx
index 85781a54..27ed1fad 100644
--- a/gstools/variogram/estimator.pyx
+++ b/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, atan2, acos, asin, sin, cos
+from libc.math cimport fabs, sqrt, atan2, acos, asin, sin, cos, M_PI
 cimport numpy as np
 
 
@@ -70,8 +70,8 @@ cdef inline bint _angle_test_2d(
     cdef double dx = x[i] - x[j]
     cdef double dy = y[i] - y[j]
     # azimuth
-    cdef double phi1 = atan2(dy,dx)
-    cdef double phi2 = atan2(-dy,-dx)
+    cdef double phi1 = atan2(dy,dx) % (2.0 * M_PI)
+    cdef double phi2 = atan2(-dy,-dx) % (2.0 * M_PI)
     # check both directions (+/-)
     cdef bint dir1 = fabs(phi1 - angles[0]) <= angles_tol
     cdef bint dir2 = fabs(phi2 - angles[0]) <= angles_tol
@@ -91,8 +91,8 @@ cdef inline bint _angle_test_3d(
     cdef double dz = z[i] - z[j]
     cdef double dr = sqrt(dx**2 + dy**2 + dz**2)
     # azimuth
-    cdef double phi1 = atan2(dy, dx)
-    cdef double phi2 = atan2(-dy, -dx)
+    cdef double phi1 = atan2(dy, dx) % (2.0 * M_PI)
+    cdef double phi2 = atan2(-dy, -dx) % (2.0 * M_PI)
     # elevation
     cdef double theta1 = acos(dz / dr)
     cdef double theta2 = acos(-dz / dr)
diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py
index 9ce2a499..d1cdf4b1 100644
--- a/gstools/variogram/variogram.py
+++ b/gstools/variogram/variogram.py
@@ -127,19 +127,17 @@ def vario_estimate_unstructured(
 
     bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0
 
-    if angles is not None:
+    if angles is not None and dim > 1:
         # there are at most (dim-1) angles
-        angles = np.array(angles, ndmin=1, dtype=np.double).ravel()[
-            : (dim - 1)
-        ]
+        angles = np.array(angles, ndmin=1, dtype=np.double)  # type convert
+        angles = angles.ravel()[: (dim - 1)]  # cutoff
         angles = np.pad(
             angles, (0, dim - angles.size - 1), "constant", constant_values=0.0
-        )
-        # correct intervalls for angles (prepared for n-d)
-        if dim >= 2:
-            angles[0] = angles[0] % (2 * np.pi)  # azimuth in [0, 2pi)
-        if dim >= 3:
-            angles[1:] = angles[1:] % np.pi  # inclinations in [0, pi)
+        )  # fill with 0 if too less given
+        # correct intervalls for angles (only need upper hemispheres)
+        angles = angles % np.pi
+    elif dim == 1:
+        angles = None
 
     if sampling_size is not None and sampling_size < len(field):
         sampled_idx = np.random.RandomState(sampling_seed).choice(

From b8955e9f55829a0cb6a37cc706690d6760dbfb39 Mon Sep 17 00:00:00 2001
From: MuellerSeb <mueller.seb@posteo.de>
Date: Fri, 6 Nov 2020 12:25:20 +0100
Subject: [PATCH 15/15] vario: fix wrong assumption about hemisphere for angles

---
 gstools/variogram/variogram.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py
index d1cdf4b1..1dcf6e40 100644
--- a/gstools/variogram/variogram.py
+++ b/gstools/variogram/variogram.py
@@ -134,8 +134,9 @@ def vario_estimate_unstructured(
         angles = np.pad(
             angles, (0, dim - angles.size - 1), "constant", constant_values=0.0
         )  # fill with 0 if too less given
-        # correct intervalls for angles (only need upper hemispheres)
-        angles = angles % np.pi
+        # correct intervalls for angles
+        angles[0] = angles[0] % (2 * np.pi)
+        angles[1:] = angles[1:] % np.pi
     elif dim == 1:
         angles = None