Merge pull request #131 from GeoStat-Framework/binning

Automatic Binning

README.md

@@ -177,26 +177,25 @@ y = np.random.RandomState(20011012).rand(1000) * 100.
 model = gs.Exponential(dim=2, var=2, len_scale=8)
 srf = gs.SRF(model, mean=0, seed=19970221)
 field = srf((x, y))
-# estimate the variogram of the field with 40 bins
-bins = np.arange(40)
-bin_center, gamma = gs.vario_estimate((x, y), field, bins)
+# estimate the variogram of the field
+bin_center, gamma = gs.vario_estimate((x, y), field)
 # fit the variogram with a stable model. (no nugget fitted)
 fit_model = gs.Stable(dim=2)
 fit_model.fit_variogram(bin_center, gamma, nugget=False)
 # output
-ax = fit_model.plot(x_max=40)
-ax.plot(bin_center, gamma)
+ax = fit_model.plot(x_max=bin_center[-1])
+ax.scatter(bin_center, gamma)
 print(fit_model)
 ```
 
 Which gives:
 
 ```python
-Stable(dim=2, var=1.92, len_scale=8.15, nugget=0.0, anis=[1.], angles=[0.], alpha=1.05)
+Stable(dim=2, var=1.85, len_scale=7.42, nugget=0.0, anis=[1.0], angles=[0.0], alpha=1.09)
 ```
 
 <p align="center">
-<img src="https://raw.githubusercontent.com/GeoStat-Framework/GSTools/master/docs/source/pics/exp_vario_fit.png" alt="Variogram" width="600px"/>
+<img src="https://github.com/GeoStat-Framework/GeoStat-Framework.github.io/raw/master/img/GS_vario_est.png" alt="Variogram" width="600px"/>
 </p>
 
 
@@ -325,6 +324,7 @@ in memory for immediate 3D plotting in Python.
 - [hankel >= 1.0.2](https://github.com/steven-murray/hankel)
 - [emcee >= 3.0.0](https://github.com/dfm/emcee)
 - [pyevtk >= 1.1.1](https://github.com/pyscience-projects/pyevtk)
+- [meshio>=4.0.3, <5.0](https://github.com/nschloe/meshio)
 
 ### Optional
 
@@ -339,7 +339,7 @@ You can contact us via <[email protected]>.
 
 ## License
 
-[LGPLv3][license_link] © 2018-2020
+[LGPLv3][license_link] © 2018-2021
 
 [pip_link]: https://pypi.org/project/gstools
 [conda_link]: https://docs.conda.io/en/latest/miniconda.html

docs/source/index.rst

@@ -205,24 +205,23 @@ model again.
     model = gs.Exponential(dim=2, var=2, len_scale=8)
     srf = gs.SRF(model, mean=0, seed=19970221)
     field = srf((x, y))
-    # estimate the variogram of the field with 40 bins
-    bins = np.arange(40)
-    bin_center, gamma = gs.vario_estimate((x, y), field, bins)
+    # estimate the variogram of the field
+    bin_center, gamma = gs.vario_estimate((x, y), field)
     # fit the variogram with a stable model. (no nugget fitted)
     fit_model = gs.Stable(dim=2)
     fit_model.fit_variogram(bin_center, gamma, nugget=False)
     # output
-    ax = fit_model.plot(x_max=40)
-    ax.plot(bin_center, gamma)
+    ax = fit_model.plot(x_max=bin_center[-1])
+    ax.scatter(bin_center, gamma)
     print(fit_model)
 
 Which gives:
 
 .. code-block:: python
 
-    Stable(dim=2, var=1.92, len_scale=8.15, nugget=0.0, anis=[1.], angles=[0.], alpha=1.05)
+    Stable(dim=2, var=1.85, len_scale=7.42, nugget=0.0, anis=[1.0], angles=[0.0], alpha=1.09)
 
-.. image:: https://raw.githubusercontent.com/GeoStat-Framework/GSTools/master/docs/source/pics/exp_vario_fit.png
+.. image:: https://raw.githubusercontent.com/GeoStat-Framework/GeoStat-Framework.github.io/master/img/GS_vario_est.png
    :width: 400px
    :align: center
 
@@ -361,6 +360,7 @@ Requirements
 - `hankel >= 1.0.2 <https://github.com/steven-murray/hankel>`_
 - `emcee >= 3.0.0 <https://github.com/dfm/emcee>`_
 - `pyevtk >= 1.1.1 <https://github.com/pyscience-projects/pyevtk>`_
+- `meshio>=4.0.3, <5.0 <https://github.com/nschloe/meshio>`_
 
 
 Optional

examples/00_misc/04_herten.py

@@ -145,7 +145,7 @@ def generate_transmissivity():
 # results reproducible, we can also set a seed.
 
 
-bins = np.linspace(0, 10, 50)
+bins = gs.standard_bins(pos=(x_u, y_u), max_dist=10)
 bin_center, gamma = gs.vario_estimate(
     (x_u, y_u),
     herten_log_trans.reshape(-1),

examples/03_variogram/05_auto_fit_variogram.py

@@ -0,0 +1,35 @@
+"""
+Fit Variogram with automatic binning
+------------------------------------
+"""
+import numpy as np
+import gstools as gs
+
+###############################################################################
+# Generate a synthetic field with an exponential model.
+
+x = np.random.RandomState(19970221).rand(1000) * 100.0
+y = np.random.RandomState(20011012).rand(1000) * 100.0
+model = gs.Exponential(dim=2, var=2, len_scale=8)
+srf = gs.SRF(model, mean=0, seed=19970221)
+field = srf((x, y))
+print(field.var())
+###############################################################################
+# Estimate the variogram of the field with automatic binning.
+
+bin_center, gamma = gs.vario_estimate((x, y), field)
+print("estimated bin number:", len(bin_center))
+print("maximal bin distance:", bin_center[-1])
+
+###############################################################################
+# Fit the variogram with a stable model (no nugget fitted).
+
+fit_model = gs.Stable(dim=2)
+fit_model.fit_variogram(bin_center, gamma, nugget=False)
+print(fit_model)
+
+###############################################################################
+# Plot the fitting result.
+
+ax = fit_model.plot(x_max=bin_center[-1])
+ax.scatter(bin_center, gamma)

examples/03_variogram/06_auto_bin_latlon.py

@@ -0,0 +1,76 @@
+"""
+Automatic binning with lat-lon data
+-----------------------------------
+
+In this example we demonstrate automatic binning for a tiny data set
+containing temperature records from Germany
+(See the detailed DWD example for more information on the data).
+
+We use a data set from 20 meteo-stations choosen randomly.
+"""
+import numpy as np
+import gstools as gs
+
+# lat, lon, temperature
+data = np.array(
+    [
+        [52.9336, 8.237, 15.7],
+        [48.6159, 13.0506, 13.9],
+        [52.4853, 7.9126, 15.1],
+        [50.7446, 9.345, 17.0],
+        [52.9437, 12.8518, 21.9],
+        [53.8633, 8.1275, 11.9],
+        [47.8342, 10.8667, 11.4],
+        [51.0881, 12.9326, 17.2],
+        [48.406, 11.3117, 12.9],
+        [49.7273, 8.1164, 17.2],
+        [49.4691, 11.8546, 13.4],
+        [48.0197, 12.2925, 13.9],
+        [50.4237, 7.4202, 18.1],
+        [53.0316, 13.9908, 21.3],
+        [53.8412, 13.6846, 21.3],
+        [54.6792, 13.4343, 17.4],
+        [49.9694, 9.9114, 18.6],
+        [51.3745, 11.292, 20.2],
+        [47.8774, 11.3643, 12.7],
+        [50.5908, 12.7139, 15.8],
+    ]
+)
+pos = data.T[:2]  # lat, lon
+field = data.T[2]  # temperature
+
+###############################################################################
+# Since the overall range of these meteo-stations is too low, we can use the
+# data-variance as additional information during the fit of the variogram.
+
+emp_v = gs.vario_estimate(pos, field, latlon=True)
+sph = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS)
+sph.fit_variogram(*emp_v, sill=np.var(field))
+ax = sph.plot(x_max=2 * np.max(emp_v[0]))
+ax.scatter(*emp_v, label="Empirical variogram")
+ax.legend()
+print(sph)
+
+###############################################################################
+# As we can see, the variogram fitting was successful and providing the data
+# variance helped finding the right length-scale.
+#
+# Now, we'll use this covariance model to interpolate the given data with
+# ordinary kriging.
+
+# enclosing box for data points
+grid_lat = np.linspace(np.min(pos[0]), np.max(pos[0]))
+grid_lon = np.linspace(np.min(pos[1]), np.max(pos[1]))
+# ordinary kriging
+krige = gs.krige.Ordinary(sph, pos, field)
+krige((grid_lat, grid_lon), mesh_type="structured")
+ax = krige.plot()
+# plotting lat on y-axis and lon on x-axis
+ax.scatter(pos[1], pos[0], 50, c=field, edgecolors="k", label="input")
+ax.legend()
+
+###############################################################################
+# Looks good, doesn't it?
+#
+# This example shows, that setting up variogram estimation and kriging routines
+# is straight forward with GSTools. ;-)

examples/08_geo_coordinates/00_field_generation.py

@@ -19,7 +19,7 @@
 model = gs.Gaussian(latlon=True, var=1, len_scale=777, rescale=gs.EARTH_RADIUS)
 
 lat = lon = range(-80, 81)
-srf = gs.SRF(model, seed=12345)
+srf = gs.SRF(model, seed=1234)
 field = srf.structured((lat, lon))
 srf.plot()
 
@@ -53,7 +53,7 @@
 # .. note::
 #
 #    Note, that the estimated variogram coincides with the yadrenko variogram,
-#    which means it depends on the great-circle distance.
+#    which means it depends on the great-circle distance given in radians.
 #
 #    Keep that in mind when defining bins: The range is at most
 #    :math:`\pi\approx 3.14`, which corresponds to the half globe.

examples/08_geo_coordinates/01_dwd_krige.py

@@ -79,8 +79,7 @@ def get_dwd_temperature():
 # As the maximal bin distance we choose 8 degrees, which corresponds to a
 # chordal length of about 900 km.
 
-bin_max = np.deg2rad(8)
-bins = np.linspace(0, bin_max, 20)
+bins = gs.standard_bins((lat, lon), max_dist=np.deg2rad(8), latlon=True)
 bin_c, vario = gs.vario_estimate((lat, lon), temp, bins, latlon=True)
 
 ###############################################################################
@@ -100,7 +99,7 @@ def get_dwd_temperature():
 
 model = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS)
 model.fit_variogram(bin_c, vario, nugget=False)
-ax = model.plot("vario_yadrenko", x_max=bin_max)
+ax = model.plot("vario_yadrenko", x_max=bins[-1])
 ax.scatter(bin_c, vario)
 print(model)
 

gstools/__init__.py

@@ -94,13 +94,14 @@
 
 Variogram Estimation
 ^^^^^^^^^^^^^^^^^^^^
-Estimate the variogram of a given field
+Estimate the variogram of a given field with these routines
 
 .. currentmodule:: gstools.variogram
 
 .. autosummary::
    vario_estimate
    vario_estimate_axis
+   standard_bins
 
 Misc
 ====
@@ -129,6 +130,7 @@
     vario_estimate_axis,
     vario_estimate_structured,
     vario_estimate_unstructured,
+    standard_bins,
 )
 from gstools.covmodel import (
     CovModel,
@@ -184,6 +186,7 @@
     "vario_estimate_axis",
     "vario_estimate_structured",
     "vario_estimate_unstructured",
+    "standard_bins",
 ]
 
 __all__ += [

gstools/tools/geometric.py

@@ -26,6 +26,7 @@
    ang2dir
    latlon2pos
    pos2latlon
+   chordal_to_great_circle
 """
 # pylint: disable=C0103
 
@@ -51,6 +52,7 @@
     "ang2dir",
     "latlon2pos",
     "pos2latlon",
+    "chordal_to_great_circle",
 ]
 
 
@@ -641,3 +643,24 @@ def pos2latlon(pos, radius=1.0, dtype=np.double):
     lat = np.arcsin(np.maximum(np.minimum(pos[2] / radius, 1.0), -1.0))
     lon = np.arctan2(pos[1], pos[0])
     return np.rad2deg((lat, lon), dtype=dtype)
+
+
+def chordal_to_great_circle(dist):
+    """
+    Calculate great circle distance corresponding to given chordal distance.
+
+    Parameters
+    ----------
+    dist : array_like
+        Chordal distance of two points on the unit-sphere.
+
+    Returns
+    -------
+    :class:`numpy.ndarray`
+        Great circle distance corresponding to given chordal distance.
+
+    Notes
+    -----
+    If given values are not in [0, 1], they will be truncated.
+    """
+    return 2 * np.arcsin(np.maximum(np.minimum(np.divide(dist, 2), 1), 0))

gstools/variogram/__init__.py

@@ -11,6 +11,12 @@
    vario_estimate
    vario_estimate_axis
 
+Binning
+^^^^^^^
+
+.. autosummary::
+   standard_bins
+
 ----
 """
 
@@ -20,10 +26,12 @@
     vario_estimate_structured,
     vario_estimate_unstructured,
 )
+from gstools.variogram.binning import standard_bins
 
 __all__ = [
     "vario_estimate",
     "vario_estimate_axis",
     "vario_estimate_unstructured",
     "vario_estimate_structured",
+    "standard_bins",
 ]