Merge pull request #50 from mmaelicke/direction

Direction

docs/reference/directionalvariogram.rst

@@ -6,7 +6,11 @@ DirectionalVariogram Class
     :members:
 
     .. automethod:: __init__
+    .. automethod:: azimuth
+    .. automethod:: tolerance
+    .. automethod:: bandwidth
+    .. automethod:: pair_field
+    .. automethod:: _calc_direction_mask_data
     .. automethod:: _triangle
-    .. automethod:: _circle
     .. automethod:: _compass
     .. automethod:: _direction_mask

docs/technical/direction.rst

@@ -64,11 +64,9 @@ The difference between the
 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.
+function to plot the effect of the search area. The method is called
+:func:`pair_field <skgstat.DirectionalVariogram.pair_field>`. Using
+random coordinates, the visualization is shown below.
 
 .. ipython:: python
     :okwarning:
@@ -88,8 +86,8 @@ Using random coordinates, the search area shapes are presented below.
         directional_model='triangle')
 
     @savefig dv1.png width=6in
-    DV.search_area(poi=3)
-
+    DV.pair_field(plt.gca())
+    
 The model can easily be changed, using the
 :func:`set_directional_model <skgstat.DirectionalVariogram.set_directional_model>`
 function:
@@ -99,112 +97,19 @@ function:
 
     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])
+    DV.set_directional_model('triangle')
+    DV.pair_field(plt.gca())
 
     @savefig dv2.png width=8in
+    DV.pair_field(plt.gca())
     fig.show()
 
+    DV.set_directional_model('compass')
 
-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
-    :okwarning:
-
-    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
+    @savefig dv3.png width=8in
+    DV.pair_field(plt.gca())
     fig.show()
 
-After moving and shifting, any type of Geometry could be generated and passed
-as the search area.
-
-.. ipython:: python
-    :okwarning:
-
-    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
 ======================
 

docs/userguide/variogram.rst

@@ -569,8 +569,9 @@ The Matérn model takes an additional smoothness paramter, that can
 change the shape of the function in between an exponential 
 model shape and a Gaussian one. 
 
-.. iypthon:: python
-
+.. ipython:: python
+    :okwarning:
+   
     xi = np.linspace(0, 100, 100)
 
     # plot a exponential and a gaussian
@@ -593,6 +594,202 @@ When direction matters
 What is 'direction'?
 --------------------
 
+The classic approach to calculate a variogram is based on the 
+assumption that covariance between observations can be related to 
+their separating distance. For this, point pairs of all observation 
+points are formed and it is assumed that they can be formed without any restriction.
+The only paramter to be influenced is a limiting distance, beyond which 
+a point pair does not make sense anymore. 
+
+This assumption might not always hold. Especially in landscapes, processes do 
+not occur randomly, but in an organized manner. This organization is often 
+directed, which can lead to stronger covariance in one direction than another.
+Therefore, another step has to be introduced before lag classes are formed.
+
+The *direction* of a variogram is then a orientation, which two points need. 
+If they are not oriented in the specified way, they will be ignored while calculating 
+a semi-variance value for a given lag class. Usually, you will specify a 
+orientation, which is called :func:`azimuth <skgstat.DirectionalVariogram.azimuth>`, 
+and a :func:`tolerance <skgstat.DirectionalVariogram.tolerance>`, which is an 
+offset from the given azimuth, at which a point pair will still be accepted.
+
+Defining orientiation
+---------------------
+
+One has to decide how orientation of two points is determined. In scikit-gstat,
+orientation between two observation points is only defined in :math:`\mathbb{R}^2`.
+We define the orientation as the **angle between the vector connecting two observation points 
+with the x-axis**.
+
+Thus, also the :func:`azimuth <skgstat.DirectionalVariogram.azimuth>` is defined as an 
+angle of the azimutal vector to the x-axis, with an 
+:func:`tolerance <skgstat.DirectionalVariogram.tolerance>` in degrees added to the 
+exact azimutal orientation clockwise and counter clockwise.
+
+The angle :math:`\Phi` between two vetors ``u,v`` is given like:
+
+.. math::
+
+    \Phi = cos^{-1}\left(\frac{u \circ v}{||u|| \cdot ||v||}\right)
+
+.. ipython:: python
+    :okwarning:
+
+    from matplotlib.patches import FancyArrowPatch as farrow
+    fig, ax = plt.subplots(1, 1, figsize=(6,4))
+    ax.arrow(0,0,2,1,color='k')
+    ax.arrow(-.1,0,3.1,0,color='k')
+    ax.set_xlim(-.1, 3)
+    ax.set_ylim(-.1,2.)
+    ax.scatter([0,2], [0,1], 50, c='r')
+    ax.annotate('A (0, 0)', (.0, .26), fontsize=14)
+    ax.annotate('B (2, 1)', (2.05,1.05), fontsize=14)
+    arrowstyle="Simple,head_width=6,head_length=12,tail_width=1"
+    ar = farrow([1.5,0], [1.25, 0.625],  color='r', connectionstyle="arc3, rad=.2", arrowstyle=arrowstyle)
+    ax.add_patch(ar)
+    @savefig sample_orientation_of_2_1.png width=6in
+    ax.annotate('26.5°', (1.5, 0.25), fontsize=14, color='r')
+
+The described definition of orientation is illustrated in the figure above. 
+There are two observation points, :math:`A (0,0)` and :math:`B (2, 1)`. To decide
+wether to account for them when calculating the semi-variance at their separating 
+distance lag, their orientation is used. Only if the direction of the varigram includes
+this orientation, the points are used. Imagine the azimuth and tolerance would be 
+``45°``, then anything between ``0°`` (East) and ``90°`` orientation would be included.
+The given example shows the orientation angle :math:`\Phi = 26.5°`, which means the 
+vector :math:`\overrightarrow{AB}` is included.
+
+Calculating orientations
+------------------------
+
+SciKit-GStat implements a slightly adaped version of the formula given in the 
+last section. It makes use of symmetric search areas (tolerance is applied clockwise 
+and counter clockwise) und therefore any calculated angle might be the result 
+of calculating the orientation of :math:`\overrightarrow{AB}` or 
+:math:`\overrightarrow{BA}`. Mathematically, these two vectors have two different 
+angles, but they are always both taken into account or omitted for a variagram 
+at the same time. Thus, it does not make a difference for variography. 
+However, it does make a difference when you try to use the orientation angles 
+directly as the containing matrix can contain the inverse angles.
+
+This can be demonstrated by an easy example. Let ``c`` be a set of points mirrored 
+along the x-axis.
+
+.. ipython:: python
+    :okwarning:
+
+    c = np.array([[0,0], [2,1], [1,2], [2, -1], [1, -2]])
+    east = np.array([1,0])
+
+We can plug these two arrays into the the formula above:
+
+.. ipython:: python
+    :okwarning:
+
+    u = c[1:]   # omit the first one
+    angles = np.degrees(np.arccos(u.dot(east) / np.sqrt(np.sum(u**2, axis=1))))
+    angles.round(1)
+
+You can see, that the both points and their mirrored counterpart have the same 
+angle to the x-axis, just like expected. This can be visualized by the plot below:
+
+.. ipython:: python
+    :okwarning:
+
+    fig, ax = plt.subplots(1, 1, figsize=(6,4))
+    ax.set_xlim(-.1, 2.25)
+    ax.set_ylim(-2.1,2.1)
+    ax.arrow(-.1,0,3.1,0,color='k')
+    for i,p in enumerate(u):
+        ax.arrow(0,0,p[0],p[1],color='r')
+        ax.annotate('%.1f°' % angles[i], (p[0] / 2, p[1] / 2 ), fontsize=14, color='r')
+    @savefig sample_orientation_of_multiple_points.png width=6in
+    ax.scatter(c[:,0], c[:,1], 50, c='r')
+
+The main difference to the internal structure storing the orientation angles for a 
+:class:`DirectionalVariogram <skgstat.DirectionalVariogram>` instance will store different
+angles.
+To use the class on only five points, we need to prevent the class from fitting, as 
+fitting on only 5 points will not work. But this does not affect the orientation calculations.
+Therefore, the :func:`fit <skgstat.DirectionalVariogram.fit>` mehtod is overwritten.
+
+.. iypthon:: python
+    :okwarning:
+
+    class TestCls(DirectionalVariogram):
+        def fit(*args, **kwargs):
+            pass
+
+    DV = TestCls(c, np.random.normal(0,1,len(c))
+    DV._calc_direction_mask_data()
+    np.degrees(DV._angles + np.pi)[:len(c) - 1]
+
+The first two points (with positive y-coordinate) show the same result. The other two, 
+with negative y-coordinates, are also calculated counter clockwise:
+
+.. ipython:: python
+    :okwarning:
+
+    360 - np.degrees(DV._angles + np.pi)[[2,3]]
+
+The :class:`DirectionalVariogram <skgstat.DirectionalVariogram>` class has a plotting 
+function to show a network graph of all point pairs that are oriented in the 
+variogram direction. But first we need to increase the tolerance as half tolerance 
+(``45° / 2 = 22.5°`` clockwise and counter clockwise) is smaller than both orientations.
+
+.. ipython:: python
+    :okwarning:
+
+    DV.tolerance = 90 
+    @savefig sample_pair_field_plot.png width=8in
+    DV.pair_field()
+
+Directional variogram
+---------------------
+
+
+.. ipython:: python
+    :okwarning:
+
+    field = np.loadtxt('data/aniso_x2.txt')
+    np.random.seed(1312)
+    coords = np.random.randint(100, size=(300,2))
+    vals = [field[_[0], _[1]] for _ in coords]
+
+The next step is to create two different variogram instances, which share the same 
+parameters, but use a different azimuth angle. One oriented to North and the 
+second one oriented to East.
+
+.. ipython:: python
+    :okwarning:
+
+    Vnorth = DirectionalVariogram(coords, vals, azimuth=90, tolerance=90, maxlag=80, n_lags=20)
+    Veast = DirectionalVariogram(coords, vals, azimuth=0, tolerance=90, maxlag=80, n_lags=20)
+    pd.DataFrame({'north':Vnorth.describe(), 'east': Veast.describe()})
+
+You can see, how the two are differing in effective range and also sill, only 
+caused by the orientation. Let's look at the experimental variogram:
+
+.. ipython:: python
 
-Space-time variography
-======================
+    fix, ax = plt.subplots(1,1,figsize=(8,6))
+    ax.plot(Vnorth.bins, Vnorth.experimental, '.--r', label='North-South')
+    ax.plot(Veast.bins, Veast.experimental, '.--b', label='East-West')
+    ax.set_xlabel('lag [m]')
+    ax.set_ylabel('semi-variance (matheron)')
+    @savefig expermiental_direcional_varigram_comparison.png width=8in
+    plt.legend(loc='upper left')
+
+The shape of both experimental variograms is very similar on the first 40 meters 
+of distance. Within this range, the apparent anisotropy is not pronounced. 
+The East-West oriented variograms also have an effective range of only about 40 meters,
+which means that in this direction the observations become statistically independent 
+at larger distances.
+For the North-South variogram the effective range is way bigger and the variogram 
+plot reveals much larger correlation lengths in that direction. The spatial 
+dependency is thus directed in North-South direction.
+
+
+To perform Kriging, you would now transform the data, especially in North-West 
+direction, unitl both variograms look the same within the effective range. 
+Finally, the Kriging result is back-transformed into the original coordinate system.

requirements.txt

@@ -3,5 +3,4 @@ numpy
 pandas
 matplotlib
 numba
-shapely
-scikit-learn
+scikit-learn