Merge 12012c6 into f60f2cb

examples/03_variogram/06_auto_bin_latlon.py

@@ -44,10 +44,10 @@
 # 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)
+emp_v = gs.vario_estimate(pos, field, latlon=True, geo_scale=gs.EARTH_RADIUS)
+sph = gs.Spherical(latlon=True, geo_scale=gs.EARTH_RADIUS)
 sph.fit_variogram(*emp_v, sill=np.var(field))
-ax = sph.plot(x_max=2 * np.max(emp_v[0]))
+ax = sph.plot("vario_yadrenko", x_max=2 * np.max(emp_v[0]))
 ax.scatter(*emp_v, label="Empirical variogram")
 ax.legend()
 print(sph)

examples/08_geo_coordinates/00_field_generation.py

@@ -8,15 +8,17 @@
 First we setup a model, with ``latlon=True``, to get the associated
 Yadrenko model.
 
-In addition, we will use the earth radius provided by :any:`EARTH_RADIUS`,
-to have a meaningful length scale in km.
+In addition, we will use the earth radius provided by :any:`EARTH_RADIUS`
+as ``geo_scale`` to have a meaningful length scale in km.
 
 To generate the field, we simply pass ``(lat, lon)`` as the position tuple
 to the :any:`SRF` class.
 """
+import numpy as np
+
 import gstools as gs
 
-model = gs.Gaussian(latlon=True, var=1, len_scale=777, rescale=gs.EARTH_RADIUS)
+model = gs.Gaussian(latlon=True, len_scale=777, geo_scale=gs.EARTH_RADIUS)
 
 lat = lon = range(-80, 81)
 srf = gs.SRF(model, seed=1234)
@@ -32,7 +34,7 @@
 #
 # As we will see, everthing went well... phew!
 
-bin_edges = [0.01 * i for i in range(30)]
+bin_edges = np.linspace(0, 777 * 3, 30)
 bin_center, emp_vario = gs.vario_estimate(
     (lat, lon),
     field,
@@ -41,11 +43,12 @@
     mesh_type="structured",
     sampling_size=2000,
     sampling_seed=12345,
+    geo_scale=gs.EARTH_RADIUS,
 )
 
-ax = model.plot("vario_yadrenko", x_max=0.3)
+ax = model.plot("vario_yadrenko", x_max=max(bin_center))
 model.fit_variogram(bin_center, emp_vario, nugget=False)
-model.plot("vario_yadrenko", ax=ax, label="fitted", x_max=0.3)
+model.plot("vario_yadrenko", ax=ax, label="fitted", x_max=max(bin_center))
 ax.scatter(bin_center, emp_vario, color="k")
 print(model)
 

examples/08_geo_coordinates/01_dwd_krige.py

@@ -76,17 +76,17 @@ def get_dwd_temperature(date="2020-06-09 12:00:00"):
 
 ###############################################################################
 # First we will estimate the variogram of our temperature data.
-# As the maximal bin distance we choose 8 degrees, which corresponds to a
-# chordal length of about 900 km.
+# As the maximal bin distance we choose 900 km.
 
-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)
+bin_center, vario = gs.vario_estimate(
+    (lat, lon), temp, latlon=True, geo_scale=gs.EARTH_RADIUS, max_dist=900
+)
 
 ###############################################################################
 # Now we can use this estimated variogram to fit a model to it.
 # Here we will use a :any:`Spherical` model. We select the ``latlon`` option
 # to use the `Yadrenko` variant of the model to gain a valid model for lat-lon
-# coordinates and we rescale it to the earth-radius. Otherwise the length
+# coordinates and we set the ``geo_scale`` to the earth-radius. Otherwise the length
 # scale would be given in radians representing the great-circle distance.
 #
 # We deselect the nugget from fitting and plot the result afterwards.
@@ -97,10 +97,10 @@ def get_dwd_temperature(date="2020-06-09 12:00:00"):
 #    still holds the ordinary routine that is not respecting the great-circle
 #    distance.
 
-model = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS)
-model.fit_variogram(bin_c, vario, nugget=False)
-ax = model.plot("vario_yadrenko", x_max=bins[-1])
-ax.scatter(bin_c, vario)
+model = gs.Spherical(latlon=True, geo_scale=gs.EARTH_RADIUS)
+model.fit_variogram(bin_center, vario, nugget=False)
+ax = model.plot("vario_yadrenko", x_max=max(bin_center))
+ax.scatter(bin_center, vario)
 print(model)
 
 ###############################################################################

examples/08_geo_coordinates/README.rst

@@ -22,35 +22,36 @@ in your desired model (see :any:`CovModel`):
 By doing so, the model will use the associated `Yadrenko` model on a sphere
 (see `here <https://onlinelibrary.wiley.com/doi/abs/10.1002/sta4.84>`_).
 The `len_scale` is given in radians to scale the arc-length.
-In order to have a more meaningful length scale, one can use the ``rescale``
+In order to have a more meaningful length scale, one can use the ``geo_scale``
 argument:
 
 .. code-block:: python
 
     import gstools as gs
 
-    model = gs.Gaussian(latlon=True, var=2, len_scale=500, rescale=gs.EARTH_RADIUS)
+    model = gs.Gaussian(latlon=True, var=2, len_scale=500, geo_scale=gs.EARTH_RADIUS)
 
 Then ``len_scale`` can be interpreted as given in km.
 
 A `Yadrenko` model :math:`C` is derived from a valid
 isotropic covariance model in 3D :math:`C_{3D}` by the following relation:
 
 .. math::
-   C(\zeta)=C_{3D}\left(2 \cdot \sin\left(\frac{\zeta}{2}\right)\right)
+   C(\zeta)=C_{3D}\left(2r \cdot \sin\left(\frac{\zeta}{2r}\right)\right)
 
 Where :math:`\zeta` is the
-`great-circle distance <https://en.wikipedia.org/wiki/Great-circle_distance>`_.
+`great-circle distance <https://en.wikipedia.org/wiki/Great-circle_distance>`_
+and :math:`r` is the ``geo_scale``.
 
 .. note::
 
    ``lat`` and ``lon`` are given in degree, whereas the great-circle distance
-   :math:`zeta` is given in radians.
+   :math:`zeta` is given in units of the ``geo_scale``.
 
-Note, that :math:`2 \cdot \sin(\frac{\zeta}{2})` is the
+Note, that :math:`2r \cdot \sin(\frac{\zeta}{2r})` is the
 `chordal distance <https://en.wikipedia.org/wiki/Chord_(geometry)>`_
-of two points on a sphere, which means we simply think of the earth surface
-as a sphere, that is cut out of the surrounding three dimensional space,
+of two points on a sphere with radius :math:`r`, which means we simply think of the
+earth surface as a sphere, that is cut out of the surrounding three dimensional space,
 when using the `Yadrenko` model.
 
 .. note::

examples/09_spatio_temporal/01_precip_1d.py

@@ -26,13 +26,11 @@
 # half daily timesteps over three months
 t = np.arange(0.0, 90.0, 0.5)
 
-# total spatio-temporal dimension
-st_dim = 1 + 1
 # space-time anisotropy ratio given in units d / km
 st_anis = 0.4
 
 # an exponential variogram with a corr. lengths of 2d and 5km
-model = gs.Exponential(dim=st_dim, var=1.0, len_scale=5.0, anis=st_anis)
+model = gs.Exponential(dim=1, time=True, var=1.0, len_scale=5.0, anis=st_anis)
 # create a spatial random field instance
 srf = gs.SRF(model, seed=seed)
 

examples/09_spatio_temporal/02_precip_2d.py

@@ -27,13 +27,11 @@
 # half daily timesteps over three months
 t = np.arange(0.0, 90.0, 0.5)
 
-# total spatio-temporal dimension
-st_dim = 2 + 1
 # space-time anisotropy ratio given in units d / km
 st_anis = 0.4
 
 # an exponential variogram with a corr. lengths of 5km, 5km, and 2d
-model = gs.Exponential(dim=st_dim, var=1.0, len_scale=5.0, anis=st_anis)
+model = gs.Exponential(dim=2, time=True, var=1.0, len_scale=5.0, anis=st_anis)
 # create a spatial random field instance
 srf = gs.SRF(model, seed=seed)
 

examples/09_spatio_temporal/03_geographic_coordinates.py

@@ -0,0 +1,36 @@
+"""
+Working with spatio-temporal lat-lon fields
+-------------------------------------------
+
+In this example, we demonstrate how to generate a spatio-temporal
+random field on geographical coordinates.
+
+First we setup a model, with ``latlon=True`` and ``time=True``,
+to get the associated spatio-temporal Yadrenko model.
+
+In addition, we will use the earth radius provided by :any:`EARTH_RADIUS`
+as ``geo_scale`` to have a meaningful length scale in km.
+
+To generate the field, we simply pass ``(lat, lon, time)`` as the position tuple
+to the :any:`SRF` class.
+
+The anisotropy factor of `0.1` (days/km) will result in a time length-scale of `100` days.
+"""
+import numpy as np
+
+import gstools as gs
+
+model = gs.Matern(
+    latlon=True,
+    time=True,
+    var=1,
+    len_scale=1000,
+    anis=0.1,
+    geo_scale=gs.EARTH_RADIUS,
+)
+
+lat = lon = np.linspace(-80, 81, 50)
+time = np.linspace(0, 777, 50)
+srf = gs.SRF(model, seed=1234)
+field = srf.structured((lat, lon, time))
+srf.plot()

pyproject.toml

@@ -147,9 +147,9 @@ target-version = [
     max-args = 20
     max-locals = 50
     max-branches = 30
-    max-statements = 80
+    max-statements = 85
     max-attributes = 25
-    max-public-methods = 75
+    max-public-methods = 80
 
 [tool.cibuildwheel]
 # Switch to using build

src/gstools/__init__.py

@@ -125,6 +125,9 @@
 
 .. autosummary::
    EARTH_RADIUS
+   KM_SCALE
+   DEGREE_SCALE
+   RADIAN_SCALE
 """
 # Hooray!
 from gstools import (
@@ -161,7 +164,10 @@
 from gstools.field import SRF, CondSRF
 from gstools.krige import Krige
 from gstools.tools import (
+    DEGREE_SCALE,
     EARTH_RADIUS,
+    KM_SCALE,
+    RADIAN_SCALE,
     generate_grid,
     generate_st_grid,
     rotated_main_axes,
@@ -188,7 +194,7 @@
 
 __all__ = ["__version__"]
 __all__ += ["covmodel", "field", "variogram", "krige", "random", "tools"]
-__all__ += ["transform", "normalizer"]
+__all__ += ["transform", "normalizer", "config"]
 __all__ += [
     "CovModel",
     "Gaussian",
@@ -226,6 +232,9 @@
     "generate_grid",
     "generate_st_grid",
     "EARTH_RADIUS",
+    "KM_SCALE",
+    "DEGREE_SCALE",
+    "RADIAN_SCALE",
     "vtk_export",
     "vtk_export_structured",
     "vtk_export_unstructured",