Merge pull request #308 from GeoStat-Framework/add_time_feature

Add 'temporal', 'spatial_dim' and 'geo_scale' attributes to Covmodel

Author: MuellerSeb

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
+- `CovModel`s 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
 
@@ -25,7 +53,7 @@ All notable changes to **GSTools** will be documented in this file.
 - better support for custom generators [#250](https://github.com/GeoStat-Framework/GSTools/pull/250) [#259](https://github.com/GeoStat-Framework/GSTools/pull/259)
 - add `valid_value_types` class variable to all field classes [#250](https://github.com/GeoStat-Framework/GSTools/pull/250)
 - PyKrige: fix passed variogram in case of latlon models [#254](https://github.com/GeoStat-Framework/GSTools/pull/254)
-- add bounds checks for optional arguments of CovModel when resetting by class attribute [#255](https://github.com/GeoStat-Framework/GSTools/pull/255)
+- add bounds checks for optional arguments of `CovModel` when resetting by class attribute [#255](https://github.com/GeoStat-Framework/GSTools/pull/255)
 - minor coverage improvements [#255](https://github.com/GeoStat-Framework/GSTools/pull/255)
 - documentation: readability improvements [#257](https://github.com/GeoStat-Framework/GSTools/pull/257)
 
@@ -181,7 +209,7 @@ See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197)
 - added new `len_rescaled` attribute to the `CovModel` class, which is the rescaled `len_scale`: `len_rescaled = len_scale / rescale`
 - new method `default_rescale` to provide default rescale factor (can be overridden)
 - remove `doctest` calls
-- docstring updates in CovModel and derived models
+- docstring updates in `CovModel` and derived models
 - updated all models to use the `cor` routine and make use of the `rescale` argument (See: [#90](https://github.com/GeoStat-Framework/GSTools/issues/90))
 - TPL models got a separate base class to not repeat code
 - added **new models** (See: [#88](https://github.com/GeoStat-Framework/GSTools/issues/88)):
@@ -208,7 +236,7 @@ See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197)
 #### Arbitrary dimensions ([#112](https://github.com/GeoStat-Framework/GSTools/issues/112))
 - allow arbitrary dimensions in all routines (CovModel, Krige, SRF, variogram)
 - anisotropy and rotation following a generalization of tait-bryan angles
-- CovModel provides `isometrize` and `anisometrize` routines to convert points
+- `CovModel` provides `isometrize` and `anisometrize` routines to convert points
 
 #### New Class for Conditioned Random Fields ([#130](https://github.com/GeoStat-Framework/GSTools/issues/130))
 - **THIS BREAKS BACKWARD COMPATIBILITY**
@@ -232,7 +260,7 @@ See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197)
 
 ### Changes
 - drop support for Python 3.5 [#146](https://github.com/GeoStat-Framework/GSTools/pull/146)
-- added a finit limit for shape-parameters in some CovModels [#147](https://github.com/GeoStat-Framework/GSTools/pull/147)
+- added a finit limit for shape-parameters in some `CovModel`s [#147](https://github.com/GeoStat-Framework/GSTools/pull/147)
 - drop usage of `pos2xyz` and `xyz2pos`
 - remove structured option from generators (structured pos need to be converted first)
 - explicitly assert dim=2,3 when generating vector fields
@@ -248,7 +276,7 @@ See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197)
 - typo in keyword argument for vario_estimate_structured [#80](https://github.com/GeoStat-Framework/GSTools/issues/80)
 - isotropic rotation of SRF was not possible [#100](https://github.com/GeoStat-Framework/GSTools/issues/100)
 - `CovModel.opt_arg` now sorted [#103](https://github.com/GeoStat-Framework/GSTools/issues/103)
-- CovModel.fit: check if weights are given as a string (numpy comparison error) [#111](https://github.com/GeoStat-Framework/GSTools/issues/111)
+- `CovModel.fit`: check if weights are given as a string (numpy comparison error) [#111](https://github.com/GeoStat-Framework/GSTools/issues/111)
 - several pylint fixes ([#159](https://github.com/GeoStat-Framework/GSTools/pull/159))
 
 
@@ -266,7 +294,7 @@ See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197)
 ### Enhancements
 - different variogram estimator functions can now be used #51
 - the TPLGaussian and TPLExponential now have analytical spectra #67
-- added property ``is_isotropic`` to CovModel #67
+- added property `is_isotropic` to `CovModel` #67
 - reworked the whole krige sub-module to provide multiple kriging methods #67
   - Simple
   - Ordinary
@@ -279,7 +307,7 @@ See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197)
 
 ### Changes
 - Python versions 2.7 and 3.4 are no longer supported #40 #43
-- CovModel: in 3D the input of anisotropy is now treated slightly different: #67
+- `CovModel`: in 3D the input of anisotropy is now treated slightly different: #67
   - single given anisotropy value [e] is converted to [1, e] (it was [e, e] before)
   - two given length-scales [l_1, l_2] are converted to [l_1, l_2, l_2] (it was [l_1, l_2, l_1] before)
 
@@ -297,31 +325,31 @@ See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197)
 
 ### Bugfixes
 - define spectral_density instead of spectrum in covariance models since Cov-base derives spectrum. See: [commit 00f2747](https://github.com/GeoStat-Framework/GSTools/commit/00f2747fd0503ff8806f2eebfba36acff813416b)
-- better boundaries for CovModel parameters. See: https://github.com/GeoStat-Framework/GSTools/issues/37
+- better boundaries for `CovModel` parameters. See: https://github.com/GeoStat-Framework/GSTools/issues/37
 
 
 ## [1.1.0] - Reverberating Red - 2019-10-01
 
 ### Enhancements
 - by using Cython for all the heavy computations, we could achieve quite some speed ups and reduce the memory consumption significantly #16
 - parallel computation in Cython is now supported with the help of OpenMP and the performance increase is nearly linear with increasing cores #16
-- new submodule ``krige`` providing simple (known mean) and ordinary (estimated mean) kriging working analogous to the srf class
-- interface to pykrige to use the gstools CovModel with the pykrige routines (https://github.com/bsmurphy/PyKrige/issues/124)
-- the srf class now provides a ``plot`` and a ``vtk_export`` routine
+- new submodule `krige` providing simple (known mean) and ordinary (estimated mean) kriging working analogous to the srf class
+- interface to pykrige to use the gstools `CovModel` with the pykrige routines (https://github.com/bsmurphy/PyKrige/issues/124)
+- the srf class now provides a `plot` and a `vtk_export` routine
 - incompressible flow fields can now be generated #14
 - new submodule providing several field transformations like: Zinn&Harvey, log-normal, bimodal, ... #13
 - Python 3.4 and 3.7 wheel support #19
 - field can now be generated directly on meshes from [meshio](https://github.com/nschloe/meshio) and [ogs5py](https://github.com/GeoStat-Framework/ogs5py), see: [commit f4a3439](https://github.com/GeoStat-Framework/GSTools/commit/f4a3439400b81d8d9db81a5f7fbf6435f603cf05)
-- the srf and kriging classes now store the last ``pos``, ``mesh_type`` and ``field`` values to keep them accessible, see: [commit 29f7f1b](https://github.com/GeoStat-Framework/GSTools/commit/29f7f1b029866379ce881f44765f72534d757fae)
+- the srf and kriging classes now store the last `pos`, `mesh_type` and `field` values to keep them accessible, see: [commit 29f7f1b](https://github.com/GeoStat-Framework/GSTools/commit/29f7f1b029866379ce881f44765f72534d757fae)
 - tutorials on all important features of GSTools have been written for you guys #20
 - a new interface to pyvista is provided to export fields to python vtk representation, which can be used for plotting, exploring and exporting fields #29
 
 ### Changes
 - the license was changed from GPL to LGPL in order to promote the use of this library #25
 - the rotation angles are now interpreted in positive direction (counter clock wise)
-- the ``force_moments`` keyword was removed from the SRF call method, it is now in provided as a field transformation #13
+- the `force_moments` keyword was removed from the SRF call method, it is now in provided as a field transformation #13
 - drop support of python implementations of the variogram estimators #18
-- the ``variogram_normed`` method was removed from the ``CovModel`` class due to redundance [commit 25b1647](https://github.com/GeoStat-Framework/GSTools/commit/25b164722ac6744ebc7e03f3c0bf1c30be1eba89)
+- the `variogram_normed` method was removed from the `CovModel` class due to redundance [commit 25b1647](https://github.com/GeoStat-Framework/GSTools/commit/25b164722ac6744ebc7e03f3c0bf1c30be1eba89)
 - the position vector of 1D fields does not have to be provided in a list-like object with length 1 [commit a6f5be8](https://github.com/GeoStat-Framework/GSTools/commit/a6f5be8bfd2db1f002e7889ecb8e9a037ea08886)
 
 ### Bugfixes
@@ -353,7 +381,7 @@ See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197)
 
 ### Changes
 - release is not downwards compatible with release v0.4.0
-- SRF creation has been adapted for the CovModel
+- SRF creation has been adapted for the `CovModel`
 - a tuple `pos` is now used instead of `x`, `y`, and `z` for the axes
 - renamed `estimate_unstructured` and `estimate_structured` to `vario_estimate_unstructured` and `vario_estimate_structured` for less ambiguity
 

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

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

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)

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)
 

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)
 
 ###############################################################################

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::

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)
 

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)