more docs on DirectionalVariogram

.zenodo.json

@@ -0,0 +1,34 @@
+{
+    "description": "SciKit-Gstat is a scipy-styled analysis module for geostatistics. It includes, two base classes Variogram and DirectionalVariogram. Both have a very similar interface and can compute experimental variograms and model variograms. The module makes use of a rich selection of semi-variance estimators and variogram model functions, while being extensible at the same time.",
+    "license": "MIT",
+    "title": "mmaelicke/scikit-gstat: Geostatistical variogram toolbox",
+    "version": "v.0.2.2",
+    "upload_type": "software",
+    "keywords": [
+        "geostatistics"
+    ],
+    "publication_date": "2018-08-15",
+    "creators": [
+        {
+            "orcid": "0000-0002-0424-2651",
+            "affiliation": "Institute for Water and River Basin Management, Karlsruhe Institute of Technology (KIT), Germany",
+            "name": "M\u00e4licke, Mirko"
+        },
+        {
+            "name": "Schneider, Helge David"
+        }
+    ],
+    "access_right": "open",
+    "related_identifiers": [
+        {
+            "scheme": "url",
+            "identifier": "https://github.com/mmaelicke/scikit-gstat/tree/v.0.1.8",
+            "relation": "isSupplementTo"
+        },
+        {
+            "scheme": "doi",
+            "identifier": "10.5281/zenodo.1345584",
+            "relation": "isVersionOf"
+        }
+    ]
+}

docs/reference/directionalvariogram.rst

@@ -5,4 +5,8 @@ DirectionalVariogram Class
 .. autoclass:: skgstat.DirectionalVariogram
     :members:
 
-    .. automethod:: __init__
+    .. automethod:: __init__
+    .. automethod:: _triangle
+    .. automethod:: _circle
+    .. automethod:: _compass
+    .. automethod:: _direction_mask

docs/technical/direction.rst

@@ -0,0 +1,222 @@
+======================
+Directional Variograms
+======================
+
+General
+=======
+
+With version 0.2.2, directional variograms have been introduced. A
+directional variogram is a variogram where point pairs are only included into
+the semivariance calculation if they fulfill a specified spatial relation.
+This relation is expressed as a *search area* that identifies all
+*directional* points for a given specific point. SciKit-GStat refers to this
+point as *poi* (point of interest). The implementation is done by the
+:class:`DirectionalVariogram <skgstat.DirectionalVariogram>` class.
+
+Understanding Search Area
+=========================
+
+.. note::
+
+    The :class:`DirectionalVariogram <skgstat.DirectionalVariogram>` is
+    in general capable of handling n-dimensional coordinates. The application
+    of directional dependency is, however, only applied to the first two
+    dimensions.
+
+Understanding the search area of a directional is vital for using the
+:class:`DirectionalVariogram <skgstat.DirectionalVariogram>` class. The
+search area is controlled by the
+:func:`directional_model <skgstat.DirectionalVariogram.directional_model>`
+property which determines the shape of the search area. The extend and
+orientation of this area is controlled by the parameters:
+
+- :func:`azimuth <skgstat.DirectionalVariogram.azimuth>`
+- :func:`tolerance <skgstat.DirectionalVariogram.tolerance>`
+- :func:`bandwidth <skgstat.DirectionalVariogram.bandwidth>`
+
+As of this writing, SciKit-GStat supports three different search area shapes:
+
+- :func:`triangle <skgstat.DirectionalVariogram._triangle>` (*default*)
+- :func:`circle <skgstat.DirectionalVariogram._circle>`
+- :func:`compass <skgstat.DirectionalVariogram._compass>`
+
+Additionally, the shape generation is controlled by the
+:func:`tolerance <skgstat.DirectionalVariogram.tolerance>` parameter
+(:func:`triangle <skgstat.DirectionalVariogram._triangle>`,
+:func:`compass <skgstat.DirectionalVariogram._compass>`) and
+:func:`bandwidth <skgstat.DirectionalVariogram.bandwidth>` parameter
+(:func:`triangle <skgstat.DirectionalVariogram._triangle>`,
+:func:`circle <skgstat.DirectionalVariogram._circle>`). The
+:func:`azimuth <skgstat.DirectionalVariogram.azimuth>` is used to rotate the
+search area into a desired direction. An azimuth of 0° is heading East of the
+coordinate plane. Positive values for azimuth rotate the search area
+clockwise, negative values counter-clockwise.
+The :func:`tolerance <skgstat.DirectionalVariogram.tolerance>` specifies how
+far the angle (against 'x-axis') between two points can be off the azimuth to
+be still considered as a directional point pair. Based on this definition,
+two points at a larger distance would generally be allowed to differ more
+from azimuth in terms of coordinate distance. Therefore the
+:func:`bandwidth <skgstat.DirectionalVariogram.bandwidth>` defines a maximum
+coordinate distance, a point can have from the azimuth line.
+The difference between the
+:func:`triangle <skgstat.DirectionalVariogram._triangle>` and the
+:func:`compass <skgstat.DirectionalVariogram._compass>` search area is that
+the triangle uses the bandwidth and the compass does not.
+
+The :class:`DirectionalVariogram <skgstat.DirectionalVariogram>` has a
+function to plot the current search area. As the search area is specific to
+the current poi, it has to be defined as the index of the coordinate to be used.
+The method is called
+:func:`search_area <skgstat.DirectionalVariogram.search_area>`.
+Using random coordinates, the search area shapes are presented below.
+
+.. ipython:: python
+
+    from skgstat import DirectionalVariogram
+    import numpy as np
+    import matplotlib.pyplot as plt
+    plt.style.use('ggplot')
+    np.random.seed(42)
+    coords = np.random.gamma(15, 6, (40, 2))
+    np.random.seed(42)
+    vals = np.random.normal(5,1, 40)
+
+    DV = DirectionalVariogram(coords, vals,
+        azimuth=0,
+        tolerance=45,
+        directional_model='triangle')
+
+    @savefig dv1.png width=6in
+    DV.search_area(poi=3)
+
+The model can easily be changed, using the
+:func:`set_directional_model <skgstat.DirectionalVariogram.set_directional_model>`
+function:
+
+.. ipython:: python
+
+    fig, axes = plt.subplots(1, 2, figsize=(8, 4))
+
+    DV.set_directional_model('circle')
+    DV.search_area(poi=3, ax=axes[0])
+
+    DV.set_directional_model('compass')
+    DV.search_area(poi=3, ax=axes[1])
+
+    @savefig dv2.png width=8in
+    fig.show()
+
+
+Local Reference System
+======================
+
+In order to apply different search area shapes and rotate them considering
+the given azimuth, a few preprocessing steps are necessary. This can lead to
+some calculation overhead, as in the case of a
+:func:`compass <skgstat.DirectionalVariogram._compass>` model.
+The selection of point pairs being directional is implemented by transforming
+the coordinates into a local reference system iteratively for each coordinate.
+For multidimensional coordinates, only the first two dimensions are used.
+They are shifted to make the current point of interest the origin of the
+local reference system. Then all other points are rotated until the azimuth
+overlays the local x-axis. This makes the definition of different shapes way
+easier.
+In this local system, the bandwidth can easily be applied to the transformed
+y-axis. The
+:func:`set_directional_model <skgstat.DirectionalVariogram.set_directional_model>`
+can also set a custom function as search area shape, that accepts the current
+local reference system and returns the search area for the given poi.
+The search area has to be returned as a shapely Polygon. Unfortunately, the
+:func:`tolerance <skgstat.DirectionalVariogram.tolerance>` and
+:func:`bandwidth <skgstat.DirectionalVariogram.bandwidth>` parameter are not
+passed yet.
+
+.. note::
+
+    For the next release, it is planned to pass all necessary parameters to
+    the directional model function. This should greatly improve the
+    definition of custom shapes. Until the implementation, the parameters
+    have to be injected directly.
+
+The following example will illustrate the rotation of the local reference
+system.
+
+.. ipython:: python
+
+    from matplotlib.patches import FancyArrowPatch as arrow
+    np.random.seed(42)
+    c = np.random.gamma(10, 6, (9, 2))
+    mid = np.array([[np.mean(c[:,0]), np.mean(c[:,1])]])
+    coords = np.append(mid, c, axis=0) - mid
+
+    DV = DirectionalVariogram(coords, vals[:10],
+        azimuth=45, tolerance=45)
+
+    loc = DV.local_reference_system(poi=0)
+
+    fig, ax = plt.subplots(1, 1, figsize=(6, 6))
+    ax.scatter(coords[:,0], coords[:,1], 15, c='b')
+    ax.scatter(loc[:,0], loc[:,1], 20, c='r')
+    ax.scatter([0], [0], 40, c='y')
+
+    for u,v in zip(coords[1:], loc[1:]):
+        arrowstyle="Simple,head_width=6,head_length=12,tail_width=1"
+        a = arrow(u, v,  color='grey', linestyle='--',
+            connectionstyle="arc3, rad=.3", arrowstyle=arrowstyle)
+        ax.add_patch(a)
+
+    @savefig transform.png width=6in
+    fig.show()
+
+After moving and shifting, any type of Geometry could be generated and passed
+as the search area.
+
+.. ipython:: python
+
+    from shapely.geometry import Polygon
+
+    def M(loc):
+        xmax = np.max(loc[:,0])
+        ymax = np.max(loc[:,1])
+        return Polygon([
+            (0, 0),
+            (0, ymax * 0.6),
+            (0.05*xmax, ymax * 0.6),
+            (xmax * 0.3, ymax * 0.3),
+            (0.55 * xmax, 0.6 * ymax),
+            (0.6 * xmax, 0.6 * ymax),
+            (0.6 * xmax, 0),
+            (0.55 * xmax, 0),
+            (0.55 * xmax, 0.55 * ymax),
+            (xmax * 0.325, ymax * 0.275),
+            (xmax * 0.275, ymax * 0.275),
+            (0.05 * xmax, 0.55 * ymax),
+            (0.05 * xmax, 0),
+            (0, 0)
+        ])
+
+    DV.set_directional_model(M)
+
+    @savefig custom.png width=6in
+    DV.search_area(poi=0)
+
+
+Directional variograms
+======================
+
+In principle, the :class:`DirectionalVariogram <skgstat.DirectionalVariogram>`
+can be used just like the :class:`Variogram <skgstat.Variogram>` base class.
+In fact :class:`DirectionalVariogram <skgstat.DirectionalVariogram>` inherits
+most of the behaviour. All the functionality described in the previous
+sections is added to the basic :class:`Variogram <skgstat.Variogram>`.
+All other methods and attributes can be used in the same way.
+
+.. warning::
+
+    In order to implement the directional dependency, some methods have been
+    rewritten in :class:`DirectionalVariogram <skgstat.DirectionalVariogram>`.
+    Thus the following methods do **not** show the same behaviour:
+
+    - :func:`DirectionalVariogram.bins <skgstat.DirectionalVariogram.bins>`
+    - :func:`DirectionalVariogram._calc_groups <skgstat.DirectionalVariogram._calc_groups>`
+

docs/technical/technical.rst

@@ -12,3 +12,4 @@ make sense in every situation.
     :caption: Contents:
 
     fitting
+    direction

skgstat/DirectionalVariogram.py

@@ -278,10 +278,6 @@ def __init__(self,
 
     @property
     def azimuth(self):
-        return self._azimuth
-
-    @azimuth.setter
-    def azimuth(self, angle):
         """Direction azimuth
 
         Main direction for te selection of points in the formation of point
@@ -298,6 +294,10 @@ def azimuth(self, angle):
         ValueError : in case angle < -180° or angle > 180
 
         """
+        return self._azimuth
+
+    @azimuth.setter
+    def azimuth(self, angle):
         if angle < -180 or angle > 180:
             raise ValueError('The azimuth is an angle in degree and has to '
                              'meet -180 <= angle <= 180')
@@ -309,26 +309,26 @@ def azimuth(self, angle):
 
     @property
     def tolerance(self):
-        return self._tolerance
-
-    @tolerance.setter
-    def tolerance(self, angle):
         """Azimuth tolerance
 
-        Tolerance angle of how far a point can be off the azimuth for being
-        still counted as directional. A tolerance angle will be applied to
-        the azimuth angle symmetrically.
+         Tolerance angle of how far a point can be off the azimuth for being
+         still counted as directional. A tolerance angle will be applied to
+         the azimuth angle symmetrically.
 
-        Parameters
-        ----------
-        angle : float
-            New tolerance angle in **degree**. Has to meet 0 <= angle <= 360.
+         Parameters
+         ----------
+         angle : float
+             New tolerance angle in **degree**. Has to meet 0 <= angle <= 360.
 
-        Raises
-        ------
-        ValueError : in case angle < 0 or angle > 360
+         Raises
+         ------
+         ValueError : in case angle < 0 or angle > 360
 
-        """
+         """
+        return self._tolerance
+
+    @tolerance.setter
+    def tolerance(self, angle):
         if angle < 0 or angle > 360:
             raise ValueError('The tolerance is an angle in degree and has to '
                              'meet 0 <= angle <= 360')
@@ -340,13 +340,6 @@ def tolerance(self, angle):
 
     @property
     def bandwidth(self):
-        if self._bandwidth is None:
-            return 0
-        else:
-            return self._bandwidth
-
-    @bandwidth.setter
-    def bandwidth(self, width):
         """Tolerance bandwidth
 
         New bandwidth parameter. As the tolerance from azimuth is given as an
@@ -364,6 +357,13 @@ def bandwidth(self, width):
         ValueError : in case width is negative
 
         """
+        if self._bandwidth is None:
+            return 0
+        else:
+            return self._bandwidth
+
+    @bandwidth.setter
+    def bandwidth(self, width):
         # check if quantiles is given
         if isinstance(width, str):
             # TODO document and handle more exceptions
@@ -601,11 +601,14 @@ def _triangle(self, local_ref):
 
         Notes
         -----
-               C
-              /|\
-           a / | \ a
-            /__|h_\
-           A   c   B
+
+        .. code-block:: text
+
+                 C
+                /|\
+             a / | \ a
+              /__|h_\
+             A   c   B
 
         The point of interest is C and c is the bandwidth. The angle at C
         (gamma) is the tolerance. From this, a and then h can be calculated.
@@ -698,7 +701,7 @@ def _circle(self, local_ref):
         return Polygon(half_circle)
 
     def _compass(self, local_ref):
-        """Compass direction Search Area
+        r"""Compass direction Search Area
 
         Construct a search area for building directional dependent point
         pairs. The compass search area will **not** be bounded by the