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

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_geografic_coordinates.py

@@ -0,0 +1,34 @@
+"""
+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`,
+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` will result in a time length-scale of `77.7` days.
+"""
+import gstools as gs
+
+model = gs.Gaussian(
+    latlon=True,
+    time=True,
+    var=1,
+    len_scale=777,
+    anis=0.1,
+    rescale=gs.EARTH_RADIUS,
+)
+
+lat = lon = range(-80, 81)
+time = range(0, 110, 10)
+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/covmodel/base.py

@@ -103,6 +103,12 @@ class CovModel:
         disabled. `rescale` can be set to e.g. earth's radius,
         to have a meaningful `len_scale` parameter.
         Default: False
+    time : :class:`bool`, optional
+        Create a metric spatio-temporal covariance model.
+        Setting this to true will increase `dim` and `field_dim` by 1.
+        `spatial_dim` will be `field_dim - 1`.
+        The time-dimension is appended, meaning the pos tuple is (x,y,z,...,t).
+        Default: False
     var_raw : :class:`float` or :any:`None`, optional
         raw variance of the model which will be multiplied with
         :any:`CovModel.var_factor` to result in the actual variance.
@@ -123,6 +129,7 @@ class CovModel:
     def __init__(
         self,
         dim=3,
+        *,
         var=1.0,
         len_scale=1.0,
         nugget=0.0,
@@ -131,6 +138,7 @@ def __init__(
         integral_scale=None,
         rescale=None,
         latlon=False,
+        time=False,
         var_raw=None,
         hankel_kw=None,
         **opt_arg,
@@ -156,11 +164,13 @@ def __init__(
         self._nugget_bounds = None
         self._anis_bounds = None
         self._opt_arg_bounds = {}
-        # Set latlon first
+        # Set latlon and time first
         self._latlon = bool(latlon)
+        self._time = bool(time)
         # SFT class will be created within dim.setter but needs hankel_kw
         self.hankel_kw = hankel_kw
-        self.dim = dim
+        # using time increases model dimension
+        self.dim = dim + int(self.time)
 
         # optional arguments for the variogram-model
         set_opt_args(self, opt_arg)
@@ -176,7 +186,9 @@ def __init__(
         # set anisotropy and len_scale, disable anisotropy for latlon models
         self._len_scale, anis = set_len_anis(self.dim, len_scale, anis)
         if self.latlon:
+            # keep time anisotropy for metric spatio-temporal model
             self._anis = np.array((self.dim - 1) * [1], dtype=np.double)
+            self._anis[-1] = anis[-1] if self.time else 1.0
             self._angles = np.array(self.dim * [0], dtype=np.double)
         else:
             self._anis = anis
@@ -530,14 +542,14 @@ def isometrize(self, pos):
         """Make a position tuple ready for isotropic operations."""
         pos = np.asarray(pos, dtype=np.double).reshape((self.field_dim, -1))
         if self.latlon:
-            return latlon2pos(pos)
+            return latlon2pos(pos, time=self.time)
         return np.dot(matrix_isometrize(self.dim, self.angles, self.anis), pos)
 
     def anisometrize(self, pos):
         """Bring a position tuple into the anisotropic coordinate-system."""
         pos = np.asarray(pos, dtype=np.double).reshape((self.dim, -1))
         if self.latlon:
-            return pos2latlon(pos)
+            return pos2latlon(pos, time=self.time)
         return np.dot(
             matrix_anisometrize(self.dim, self.angles, self.anis), pos
         )
@@ -862,6 +874,11 @@ def arg_bounds(self):
         res.update(self.opt_arg_bounds)
         return res
 
+    @property
+    def time(self):
+        """:class:`bool`: Whether the model is a metric spatio-temporal one."""
+        return self._time
+
     # geographical coordinates related
 
     @property
@@ -871,8 +888,13 @@ def latlon(self):
 
     @property
     def field_dim(self):
-        """:class:`int`: The field dimension of the model."""
-        return 2 if self.latlon else self.dim
+        """:class:`int`: The (parametric) field dimension of the model (with time)."""
+        return 2 + int(self._time) if self.latlon else self.dim
+
+    @property
+    def spatial_dim(self):
+        """:class:`int`: The spatial field dimension of the model (without time)."""
+        return 2 if self.latlon else self.dim - int(self._time)
 
     # standard parameters
 

src/gstools/covmodel/fit.py

@@ -522,6 +522,7 @@ def logistic_weights(p=0.1, mean=0.7):  # pragma: no cover
     callable
         Weighting function.
     """
+
     # define the callable weights function
     def func(x_data):
         """Callable function for the weights."""

src/gstools/covmodel/tools.py

@@ -498,13 +498,13 @@ def set_dim(model, dim):
             AttributeWarning,
         )
         dim = model.fix_dim()
-        if model.latlon and dim != 3:
+        if model.latlon and dim != (3 + int(model.time)):
             raise ValueError(
                 f"{model.name}: using fixed dimension {model.fix_dim()}, "
-                "which is not compatible with a latlon model."
+                f"which is not compatible with a latlon model (with time={model.time})."
             )
-    # force dim=3 for latlon models
-    dim = 3 if model.latlon else dim
+    # force dim=3 (or 4 when time=True) for latlon models
+    dim = (3 + int(model.time)) if model.latlon else dim
     # set the dimension
     if dim < 1:
         raise ValueError("Only dimensions of d >= 1 are supported.")
@@ -551,6 +551,7 @@ def compare(this, that):
     equal &= np.all(np.isclose(this.angles, that.angles))
     equal &= np.isclose(this.rescale, that.rescale)
     equal &= this.latlon == that.latlon
+    equal &= this.time == that.time
     for opt in this.opt_arg:
         equal &= np.isclose(getattr(this, opt), getattr(that, opt))
     return equal

src/gstools/field/base.py

@@ -678,6 +678,11 @@ def latlon(self):
         """:class:`bool`: Whether the field depends on geographical coords."""
         return False if self.model is None else self.model.latlon
 
+    @property
+    def time(self):
+        """:class:`bool`: Whether the field depends on time."""
+        return False if self.model is None else self.model.time
+
     @property
     def name(self):
         """:class:`str`: The name of the class."""

src/gstools/field/plot.py

@@ -54,7 +54,14 @@ def plot_field(
     if fld.dim == 1:
         return plot_1d(fld.pos, fld[field], fig, ax, **kwargs)
     return plot_nd(
-        fld.pos, fld[field], fld.mesh_type, fig, ax, fld.latlon, **kwargs
+        fld.pos,
+        fld[field],
+        fld.mesh_type,
+        fig,
+        ax,
+        fld.latlon,
+        fld.time,
+        **kwargs,
     )
 
 
@@ -104,6 +111,7 @@ def plot_nd(
     fig=None,
     ax=None,
     latlon=False,
+    time=False,
     resolution=128,
     ax_names=None,
     aspect="quad",
@@ -136,6 +144,11 @@ def plot_nd(
         ValueError will be raised, if a direction was specified.
         Bin edges need to be given in radians in this case.
         Default: False
+    time : :class:`bool`, optional
+        Indicate a metric spatio-temporal covariance model.
+        The time-dimension is assumed to be appended,
+        meaning the pos tuple is (x,y,z,...,t) or (lat, lon, t).
+        Default: False
     resolution : :class:`int`, optional
         Resolution of the imshow plot. The default is 128.
     ax_names : :class:`list` of :class:`str`, optional
@@ -159,14 +172,20 @@ def plot_nd(
     """
     dim = len(pos)
     assert dim > 1
-    assert not latlon or dim == 2
+    assert not latlon or dim == 2 + int(bool(time))
     if dim == 2 and contour_plot:
         return _plot_2d(
             pos, field, mesh_type, fig, ax, latlon, ax_names, **kwargs
         )
-    pos = pos[::-1] if latlon else pos
-    field = field.T if (latlon and mesh_type != "unstructured") else field
-    ax_names = _ax_names(dim, latlon, ax_names)
+    if latlon:
+        # swap lat-lon to lon-lat (x-y)
+        if time:
+            pos = (pos[1], pos[0], pos[2])
+        else:
+            pos = (pos[1], pos[0])
+        if mesh_type != "unstructured":
+            field = np.moveaxis(field, [0, 1], [1, 0])
+    ax_names = _ax_names(dim, latlon, time, ax_names)
     # init planes
     planes = rotation_planes(dim)
     plane_names = [f" {ax_names[p[0]]} - {ax_names[p[1]]}" for p in planes]
@@ -323,15 +342,20 @@ def plot_vec_field(fld, field="field", fig=None, ax=None):  # pragma: no cover
     return ax
 
 
-def _ax_names(dim, latlon=False, ax_names=None):
+def _ax_names(dim, latlon=False, time=False, ax_names=None):
+    t_fac = int(bool(time))
     if ax_names is not None:
         assert len(ax_names) >= dim
         return ax_names[:dim]
-    if dim == 2 and latlon:
-        return ["lon", "lat"]
+    if dim == 2 + t_fac and latlon:
+        return ["lon", "lat"] + t_fac * ["time"]
     if dim <= 3:
-        return ["$x$", "$y$", "$z$"][:dim] + (dim == 1) * ["field"]
-    return [f"$x_{{{i}}}$" for i in range(dim)]
+        return (
+            ["$x$", "$y$", "$z$"][:dim]
+            + t_fac * ["time"]
+            + (dim == 1) * ["field"]
+        )
+    return [f"$x_{{{i}}}$" for i in range(dim - t_fac)] + t_fac * ["time"]
 
 
 def _plot_2d(

src/gstools/tools/geometric.py

@@ -624,60 +624,75 @@ def ang2dir(angles, dtype=np.double, dim=None):
     return vec
 
 
-def latlon2pos(latlon, radius=1.0, dtype=np.double):
+def latlon2pos(latlon, radius=1.0, dtype=np.double, time=False):
     """Convert lat-lon geo coordinates to 3D position tuple.
 
     Parameters
     ----------
     latlon : :class:`list` of :class:`numpy.ndarray`
         latitude and longitude given in degrees.
+        May includes an appended time axis if `time=True`.
     radius : :class:`float`, optional
         Earth radius. Default: `1.0`
     dtype : data-type, optional
         The desired data-type for the array.
         If not given, then the type will be determined as the minimum type
         required to hold the objects in the sequence. Default: None
+    time : bool, optional
+        Whether latlon includes an appended time axis.
+        Default: False
 
     Returns
     -------
     :class:`numpy.ndarray`
         the 3D position array
     """
-    latlon = np.asarray(latlon, dtype=dtype).reshape((2, -1))
+    latlon = np.asarray(latlon, dtype=dtype).reshape((3 if time else 2, -1))
+    if time:
+        timeax = latlon[2]
+        latlon = latlon[:2]
     lat, lon = np.deg2rad(latlon)
-    return np.array(
-        (
-            radius * np.cos(lat) * np.cos(lon),
-            radius * np.cos(lat) * np.sin(lon),
-            radius * np.sin(lat) * np.ones_like(lon),
-        ),
-        dtype=dtype,
+    pos_tuple = (
+        radius * np.cos(lat) * np.cos(lon),
+        radius * np.cos(lat) * np.sin(lon),
+        radius * np.sin(lat) * np.ones_like(lon),
     )
+    if time:
+        return np.array(pos_tuple + (timeax,), dtype=dtype)
+    return np.array(pos_tuple, dtype=dtype)
 
 
-def pos2latlon(pos, radius=1.0, dtype=np.double):
+def pos2latlon(pos, radius=1.0, dtype=np.double, time=False):
     """Convert 3D position tuple from sphere to lat-lon geo coordinates.
 
     Parameters
     ----------
     pos : :class:`list` of :class:`numpy.ndarray`
         The position tuple containing points on a unit-sphere.
+        May includes an appended time axis if `time=True`.
     radius : :class:`float`, optional
         Earth radius. Default: `1.0`
     dtype : data-type, optional
         The desired data-type for the array.
         If not given, then the type will be determined as the minimum type
         required to hold the objects in the sequence. Default: None
+    time : bool, optional
+        Whether latlon includes an appended time axis.
+        Default: False
 
     Returns
     -------
     :class:`numpy.ndarray`
         the 3D position array
     """
-    pos = np.asarray(pos, dtype=dtype).reshape((3, -1))
+    pos = np.asarray(pos, dtype=dtype).reshape((4 if time else 3, -1))
     # prevent numerical errors in arcsin
     lat = np.arcsin(np.maximum(np.minimum(pos[2] / radius, 1.0), -1.0))
     lon = np.arctan2(pos[1], pos[0])
+    if time:
+        timeax = pos[3]
+        lat, lon = np.rad2deg((lat, lon), dtype=dtype)
+        return np.array((lat, lon, timeax), dtype=dtype)
     return np.rad2deg((lat, lon), dtype=dtype)