CovModel: prevent rotation between spatial and temporal dims

src/gstools/covmodel/base.py

@@ -28,6 +28,7 @@
     set_arg_bounds,
     set_dim,
     set_len_anis,
+    set_model_angles,
     set_opt_args,
     spectral_rad_pdf,
 )
@@ -39,7 +40,6 @@
     matrix_isometrize,
     pos2latlon,
     rotated_main_axes,
-    set_angles,
 )
 
 __all__ = ["CovModel"]
@@ -194,16 +194,15 @@ def __init__(
         # set parameters
         self.rescale = rescale
         self._nugget = float(nugget)
+
         # 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.temporal else 1.0
-            self._angles = np.array(self.dim * [0], dtype=np.double)
-        else:
-            self._anis = anis
-            self._angles = set_angles(self.dim, angles)
+        self._len_scale, self._anis = set_len_anis(
+            self.dim, len_scale, anis, self.latlon, self.temporal
+        )
+        self._angles = set_model_angles(
+            self.dim, angles, self.latlon, self.temporal
+        )
+
         # set var at last, because of the var_factor (to be right initialized)
         if var_raw is None:
             self._var = None
@@ -1004,12 +1003,9 @@ def anis(self):
 
     @anis.setter
     def anis(self, anis):
-        if self.latlon:
-            self._anis = np.array((self.dim - 1) * [1], dtype=np.double)
-        else:
-            self._len_scale, self._anis = set_len_anis(
-                self.dim, self.len_scale, anis
-            )
+        self._len_scale, self._anis = set_len_anis(
+            self.dim, self.len_scale, anis, self.latlon, self.temporal
+        )
         self.check_arg_bounds()
 
     @property
@@ -1019,10 +1015,9 @@ def angles(self):
 
     @angles.setter
     def angles(self, angles):
-        if self.latlon:
-            self._angles = np.array(self.dim * [0], dtype=np.double)
-        else:
-            self._angles = set_angles(self.dim, angles)
+        self._angles = set_model_angles(
+            self.dim, angles, self.latlon, self.temporal
+        )
         self.check_arg_bounds()
 
     @property

src/gstools/covmodel/tools.py

@@ -30,14 +30,15 @@
 from scipy import special as sps
 from scipy.optimize import root
 
-from gstools.tools.geometric import set_angles, set_anis
+from gstools.tools.geometric import no_of_angles, set_angles, set_anis
 from gstools.tools.misc import list_format
 
 __all__ = [
     "AttributeWarning",
     "rad_fac",
     "set_opt_args",
     "set_len_anis",
+    "set_model_angles",
     "check_bounds",
     "check_arg_in_bounds",
     "default_arg_from_bounds",
@@ -183,7 +184,7 @@ def set_opt_args(model, opt_arg):
         setattr(model, opt_name, float(opt_arg[opt_name]))
 
 
-def set_len_anis(dim, len_scale, anis):
+def set_len_anis(dim, len_scale, anis, latlon=False, temporal=False):
     """Set the length scale and anisotropy factors for the given dimension.
 
     Parameters
@@ -194,6 +195,13 @@ def set_len_anis(dim, len_scale, anis):
         the length scale of the SRF in x direction or in x- (y-, ...) direction
     anis : :class:`float` or :class:`list`
         the anisotropy of length scales along the transversal axes
+    latlon : :class:`bool`, optional
+        Whether the model is describing 2D fields on earths surface described
+        by latitude and longitude.
+        Default: False
+    temporal : :class:`bool`, optional
+        Whether a time-dimension is appended.
+        Default: False
 
     Returns
     -------
@@ -235,9 +243,47 @@ def set_len_anis(dim, len_scale, anis):
             raise ValueError(
                 "anisotropy-ratios needs to be > 0, " + "got: " + str(out_anis)
             )
+    # no anisotropy for latlon (only when temporal, but then only in time-dimension)
+    if latlon:
+        out_anis[: dim - (2 if temporal else 1)] = 1.0
     return out_len_scale, out_anis
 
 
+def set_model_angles(dim, angles, latlon=False, temporal=False):
+    """Set the model angles for the given dimension.
+
+    Parameters
+    ----------
+    dim : :class:`int`
+        spatial dimension
+    angles : :class:`float` or :class:`list`
+        the angles of the SRF
+    latlon : :class:`bool`, optional
+        Whether the model is describing 2D fields on earths surface described
+        by latitude and longitude.
+        Default: False
+    temporal : :class:`bool`, optional
+        Whether a time-dimension is appended.
+        Default: False
+
+    Returns
+    -------
+    angles : :class:`float`
+        the angles fitting to the dimension
+
+    Notes
+    -----
+        If too few angles are given, they are filled up with `0`.
+    """
+    if latlon:
+        return np.array(no_of_angles(dim) * [0], dtype=np.double)
+    out_angles = set_angles(dim, angles)
+    if temporal:
+        # no rotation between spatial dimensions and temporal dimension
+        out_angles[no_of_angles(dim - 1) :] = 0.0
+    return out_angles
+
+
 def check_bounds(bounds):
     """
     Check if given bounds are valid.
@@ -522,7 +568,9 @@ def set_dim(model, dim):
             model.dim, model._len_scale, model._anis
         )
     if model._angles is not None:
-        model._angles = set_angles(model.dim, model._angles)
+        model._angles = set_model_angles(
+            model.dim, model._angles, model.latlon, model.temporal
+        )
     model.check_arg_bounds()