GeoStat-Framework · MuellerSeb · Jun 15, 2023 · Jun 6, 2023 · Jun 6, 2023 · Jun 6, 2023
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,34 @@
 
 All notable changes to **GSTools** will be documented in this file.
 
+## [Unreleased] - ? - 2023-?
+
+### Enhancements
+- added `temporal` flag to CovModel to explicitly specify spatio-temporal models [#308](https://github.com/GeoStat-Framework/GSTools/pull/308)
+  - rotation between spatial and temporal dimension will be ignored
+  - added `spatial_dim` to CovModel to explicitly set spatial dimension for spatio-temporal models
+    - if not using `spatial_dim`, the provided `dim` needs to include the possible temporal dimension
+    - `spatial_dim` is always one less than `field_dim` for spatio-temporal models
+  - also works with `latlon=True` to have a spatio-temporal model with geographic coordinates
+  - all plotting routines respect this
+  - the `Field` class now has a `temporal` attribute which forwards the model attribute
+  - automatic variogram fitting in kriging classes for `temporal=True` and `latlon=True` will raise an error
+- added `geo_scale` to CovModel to have a more consistent way to set the units of the model length scale for geographic coordinates [#308](https://github.com/GeoStat-Framework/GSTools/pull/308)
+  - no need to use `rescale` for this anymore (was rather a hack)
+  - added `gs.KM_SCALE` which is the same as `gs.EARTH_RADIUS` for kilometer scaling
+  - added `gs.DEGREE_SCALE` for great circle distance in degrees
+  - added `gs.RADIAN_SCALE` for great circle distance in radians (default and previous behavior)
+  - yadrenko variogram respects this and assumes the great circle distances is given in the respective unit
+  - `vario_estimate` also has `geo_scale` now to control the units of the bins
+- `vario_estimate` now forwards additional kwargs to `standard_bins` (`bin_no`, `max_dist`) [#308](https://github.com/GeoStat-Framework/GSTools/pull/308)
+
+### Changes
+- CovModels expect special arguments by keyword now [#308](https://github.com/GeoStat-Framework/GSTools/pull/308)
+- always use f-strings internally [#283](https://github.com/GeoStat-Framework/GSTools/pull/283)
+
+### Bugfixes
+- latex equations were not rendered correctly in docs [#290](https://github.com/GeoStat-Framework/GSTools/pull/290)
+
 
 ## [1.4.1] - Sassy Sapphire - 2022-11
 

diff --git a/README.md b/README.md
@@ -147,7 +147,7 @@ import cartopy.crs as ccrs
 import gstools as gs
 # define a structured field by latitude and longitude
 lat = lon = range(-80, 81)
-model = gs.Gaussian(latlon=True, len_scale=777, rescale=gs.EARTH_RADIUS)
+model = gs.Gaussian(latlon=True, len_scale=777, geo_scale=gs.KM_SCALE)
 srf = gs.SRF(model, seed=12345)
 field = srf.structured((lat, lon))
 # Orthographic plotting with cartopy

diff --git a/docs/source/index.rst b/docs/source/index.rst
@@ -200,7 +200,7 @@ This works perfectly well with `cartopy <https://scitools.org.uk/cartopy/docs/la
     import gstools as gs
     # define a structured field by latitude and longitude
     lat = lon = range(-80, 81)
-    model = gs.Gaussian(latlon=True, len_scale=777, rescale=gs.EARTH_RADIUS)
+    model = gs.Gaussian(latlon=True, len_scale=777, geo_scale=gs.KM_SCALE)
     srf = gs.SRF(model, seed=12345)
     field = srf.structured((lat, lon))
     # Orthographic plotting with cartopy

diff --git a/examples/03_variogram/06_auto_bin_latlon.py b/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.KM_SCALE)
+sph = gs.Spherical(latlon=True, geo_scale=gs.KM_SCALE)
 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)

diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py
@@ -8,15 +8,19 @@
 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 a kilometer scale provided by :any:`KM_SCALE`
+as ``geo_scale`` to have a meaningful length scale in km.
+By default the length scale would be given in radians (:any:`RADIAN_SCALE`).
+A third option is a length scale in degrees (:any:`DEGREE_SCALE`).
 
 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.KM_SCALE)
 
 lat = lon = range(-80, 81)
 srf = gs.SRF(model, seed=1234)
@@ -32,7 +36,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 +45,12 @@
     mesh_type="structured",
     sampling_size=2000,
     sampling_seed=12345,
+    geo_scale=gs.KM_SCALE,
 )
 
-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)
 

diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/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.KM_SCALE, 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.KM_SCALE)
+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)
 
 ###############################################################################

diff --git a/examples/08_geo_coordinates/README.rst b/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.KM_SCALE)
 
 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::

diff --git a/examples/09_spatio_temporal/01_precip_1d.py b/examples/09_spatio_temporal/01_precip_1d.py
@@ -26,13 +26,13 @@
 # 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(
+    temporal=True, spatial_dim=1, var=1, len_scale=5, anis=st_anis
+)
 # create a spatial random field instance
 srf = gs.SRF(model, seed=seed)
 

diff --git a/examples/09_spatio_temporal/02_precip_2d.py b/examples/09_spatio_temporal/02_precip_2d.py
@@ -27,13 +27,13 @@
 # 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(
+    temporal=True, spatial_dim=2, var=1, len_scale=5, anis=st_anis
+)
 # create a spatial random field instance
 srf = gs.SRF(model, seed=seed)
 

diff --git a/examples/09_spatio_temporal/03_geographic_coordinates.py b/examples/09_spatio_temporal/03_geographic_coordinates.py
@@ -0,0 +1,37 @@
+"""
+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 ``temporal=True``,
+to get the associated spatio-temporal Yadrenko model.
+
+In addition, we will use a kilometer scale provided by :any:`KM_SCALE`
+as ``geo_scale`` to have a meaningful length scale in km.
+By default the length scale would be given in radians (:any:`RADIAN_SCALE`).
+A third option is a length scale in degrees (:any:`DEGREE_SCALE`).
+
+To generate the field, we simply pass ``(lat, lon, time)`` as the position tuple
+to the :any:`SRF` class.
+
+We will set a spatial length-scale of `1000` and a time length-scale of `100` days.
+"""
+import numpy as np
+
+import gstools as gs
+
+model = gs.Matern(
+    latlon=True,
+    temporal=True,
+    var=1,
+    len_scale=[1000, 100],
+    geo_scale=gs.KM_SCALE,
+)
+
+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()
diff --git a/examples/09_spatio_temporal/README.rst b/examples/09_spatio_temporal/README.rst
@@ -5,28 +5,47 @@ Spatio-Temporal modelling can provide insights into time dependent processes
 like rainfall, air temperature or crop yield.
 
 GSTools provides the metric spatio-temporal model for all covariance models
-by enhancing the spatial model dimension with a time dimension to result in
-the spatio-temporal dimension ``st_dim`` and setting a
-spatio-temporal anisotropy ratio with ``st_anis``:
+by setting ``temporal=True``, which enhances the spatial model dimension with
+a time dimension to result in the spatio-temporal dimension.
+Since the model dimension is then higher than the spatial dimension, you can use
+the ``spatial_dim`` argument to explicitly set the spatial dimension.
+Doing that and setting a spatio-temporal anisotropy ratio looks like this:
 
 .. code-block:: python
 
     import gstools as gs
     dim = 3  # spatial dimension
-    st_dim = dim + 1
     st_anis = 0.4
-    st_model = gs.Exponential(dim=st_dim, anis=st_anis)
+    st_model = gs.Exponential(temporal=True, spatial_dim=dim, anis=st_anis)
 
-Since it is given in the name "spatio-temporal",
-we will always treat the time as last dimension.
-This enables us to have spatial anisotropy and rotation defined as in
+Since it is given in the name "spatio-temporal", time is always treated as last dimension.
+You could also use ``dim`` to specify the dimension but note that it needs to include
+the temporal dimension.
+
+There are now three different dimension attributes giving information about (i) the
+model dimension (``dim``), (ii) the field dimension (``field_dim``, including time) and
+(iii) the spatial dimension (``spatial_dim`` always 1 less than ``field_dim`` for temporal models).
+Model and field dimension can differ in case of geographic coordinates where the model dimension is 3,
+but the field or parametric dimension is 2.
+If the model is spatio-temporal with geographic coordinates, the model dimension is 4,
+the field dimension is 3 and the spatial dimension is 2.
+
+In the case above we get:
+
+.. code-block:: python
+
+    st_model.dim == 4
+    st_model.field_dim == 4
+    st_model.spatial_dim == 3
+
+This formulation enables us to have spatial anisotropy and rotation defined as in
 non-temporal models, without altering the behavior in the time dimension:
 
 .. code-block:: python
 
     anis = [0.4, 0.2]  # spatial anisotropy in 3D
     angles = [0.5, 0.4, 0.3]  # spatial rotation in 3D
-    st_model = gs.Exponential(dim=st_dim, anis=anis+[st_anis], angles=angles)
+    st_model = gs.Exponential(temporal=True, spatial_dim=dim, anis=anis+[st_anis], angles=angles)
 
 In order to generate spatio-temporal position tuples, GSTools provides a
 convenient function :any:`generate_st_grid`. The output can be used for

diff --git a/pyproject.toml b/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

diff --git a/src/gstools/__init__.py b/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",