Source code for gstools.covmodel.base

"""
GStools subpackage providing the base class for covariance models.

.. currentmodule:: gstools.covmodel.base

The following classes are provided

.. autosummary::
   CovModel
"""
# pylint: disable=C0103, R0201, E1101, C0302, W0613
import copy

import numpy as np
from hankel import SymmetricFourierTransform as SFT
from scipy.integrate import quad as integral

from gstools.covmodel import plot
from gstools.covmodel.fit import fit_variogram
from gstools.covmodel.tools import (
    _init_subclass,
    check_arg_bounds,
    check_bounds,
    compare,
    default_arg_from_bounds,
    model_repr,
    percentile_scale,
    set_arg_bounds,
    set_dim,
    set_len_anis,
    set_model_angles,
    set_opt_args,
    spectral_rad_pdf,
)
from gstools.tools import RADIAN_SCALE
from gstools.tools.geometric import (
    great_circle_to_chordal,
    latlon2pos,
    matrix_anisometrize,
    matrix_isometrize,
    pos2latlon,
    rotated_main_axes,
)

__all__ = ["CovModel"]

# default arguments for hankel.SymmetricFourierTransform
HANKEL_DEFAULT = {"a": -1, "b": 1, "N": 200, "h": 0.001, "alt": True}


[docs]class CovModel: r"""Base class for the GSTools covariance models. Parameters ---------- dim : :class:`int`, optional dimension of the model. Includes the temporal dimension if temporal is true. To specify only the spatial dimension in that case, use `spatial_dim`. Default: ``3`` var : :class:`float`, optional variance of the model (the nugget is not included in "this" variance) Default: ``1.0`` len_scale : :class:`float` or :class:`list`, optional length scale of the model. If a single value is given, the same length-scale will be used for every direction. If multiple values (for main and transversal directions) are given, `anis` will be recalculated accordingly. If only two values are given in 3D, the latter one will be used for both transversal directions. Default: ``1.0`` nugget : :class:`float`, optional nugget of the model. Default: ``0.0`` anis : :class:`float` or :class:`list`, optional anisotropy ratios in the transversal directions [e_y, e_z]. * e_y = l_y / l_x * e_z = l_z / l_x If only one value is given in 3D, e_y will be set to 1. This value will be ignored, if multiple len_scales are given. Default: ``1.0`` angles : :class:`float` or :class:`list`, optional angles of rotation (given in rad): * in 2D: given as rotation around z-axis * in 3D: given by yaw, pitch, and roll (known as Tait–Bryan angles) Default: ``0.0`` integral_scale : :class:`float` or :class:`list` or :any:`None`, optional If given, ``len_scale`` will be ignored and recalculated, so that the integral scale of the model matches the given one. Default: :any:`None` rescale : :class:`float` or :any:`None`, optional Optional rescaling factor to divide the length scale with. This could be used for unit conversion or rescaling the length scale to coincide with e.g. the integral scale. Will be set by each model individually. Default: :any:`None` latlon : :class:`bool`, optional Whether the model is describing 2D fields on earths surface described by latitude and longitude. When using this, the model will internally use the associated 'Yadrenko' model to represent a valid model. This means, the spatial distance :math:`r` will be replaced by :math:`2\sin(\alpha/2)`, where :math:`\alpha` is the great-circle distance, which is equal to the spatial distance of two points in 3D. As a consequence, `dim` will be set to `3` and anisotropy will be disabled. `geo_scale` can be set to e.g. earth's radius, to have a meaningful `len_scale` parameter. Default: False geo_scale : :class:`float`, optional Geographic unit scaling in case of latlon coordinates to get a meaningful length scale unit. By default, len_scale is assumed to be in radians with latlon=True. Can be set to :any:`KM_SCALE` to have len_scale in km or :any:`DEGREE_SCALE` to have len_scale in degrees. Default: :any:`RADIAN_SCALE` temporal : :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 spatial_dim : :class:`int`, optional spatial dimension of the model. If given, the model dimension will be determined from this spatial dimension and the possible temporal dimension if temporal is ture. Default: None 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. If given, ``var`` will be ignored. (This is just for models that override :any:`CovModel.var_factor`) Default: :any:`None` hankel_kw: :class:`dict` or :any:`None`, optional Modify the init-arguments of :any:`hankel.SymmetricFourierTransform` used for the spectrum calculation. Use with caution (Better: Don't!). ``None`` is equivalent to ``{"a": -1, "b": 1, "N": 1000, "h": 0.001}``. Default: :any:`None` **opt_arg Optional arguments are covered by these keyword arguments. If present, they are described in the section `Other Parameters`. """ def __init__( self, dim=3, var=1.0, len_scale=1.0, nugget=0.0, anis=1.0, angles=0.0, *, integral_scale=None, rescale=None, latlon=False, geo_scale=RADIAN_SCALE, temporal=False, spatial_dim=None, var_raw=None, hankel_kw=None, **opt_arg, ): # assert, that we use a subclass # this is the case, if __init_subclass__ is called, which creates # the "variogram"... so we check for that if not hasattr(self, "variogram"): raise TypeError("Don't instantiate 'CovModel' directly!") # prepare dim setting self._dim = None self._hankel_kw = None self._sft = None # prepare parameters (they are checked in dim setting) self._rescale = None self._len_scale = None self._anis = None self._angles = None # prepare parameters boundaries self._var_bounds = None self._len_scale_bounds = None self._nugget_bounds = None self._anis_bounds = None self._opt_arg_bounds = {} # Set latlon and temporal first self._latlon = bool(latlon) self._temporal = bool(temporal) self._geo_scale = abs(float(geo_scale)) # SFT class will be created within dim.setter but needs hankel_kw self.hankel_kw = hankel_kw # using time increases model dimension given by "spatial_dim" self.dim = ( dim if spatial_dim is None else spatial_dim + int(self.temporal) ) # optional arguments for the variogram-model set_opt_args(self, opt_arg) # set standard boundaries for variance, len_scale, nugget and opt_arg bounds = self.default_arg_bounds() bounds.update(self.default_opt_arg_bounds()) self.set_arg_bounds(check_args=False, **bounds) # set parameters self.rescale = rescale self._nugget = float(nugget) # set anisotropy and len_scale, disable anisotropy for latlon models self._len_scale, self._anis = set_len_anis( self.dim, len_scale, anis, self.latlon ) 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 self.var = var else: self._var = float(var_raw) self._integral_scale = None self.integral_scale = integral_scale # set var again, if int_scale affects var_factor if var_raw is None: self._var = None self.var = var else: self._var = float(var_raw) # final check for parameter bounds self.check_arg_bounds() # additional checks for the optional arguments (provided by user) self.check_opt_arg() # precision for printing self._prec = 3 # one of these functions needs to be overridden def __init_subclass__(cls): """Initialize gstools covariance model.""" _init_subclass(cls) # modify the docstrings: class docstring gets attributes added if cls.__doc__ is None: cls.__doc__ = "User defined GSTools Covariance-Model." cls.__doc__ += CovModel.__doc__[45:] # overridden functions get standard doc if no new doc was created ign = ["__", "variogram", "covariance", "cor"] for att, attr_cls in cls.__dict__.items(): if any(att.startswith(i) for i in ign) or att not in dir(CovModel): continue attr_doc = getattr(CovModel, att).__doc__ if attr_cls.__doc__ is None: attr_cls.__doc__ = attr_doc # special variogram functions
[docs] def vario_axis(self, r, axis=0): r"""Variogram along axis of anisotropy.""" if axis == 0: return self.variogram(r) return self.variogram(np.abs(r) / self.anis[axis - 1])
[docs] def cov_axis(self, r, axis=0): r"""Covariance along axis of anisotropy.""" if axis == 0: return self.covariance(r) return self.covariance(np.abs(r) / self.anis[axis - 1])
[docs] def cor_axis(self, r, axis=0): r"""Correlation along axis of anisotropy.""" if axis == 0: return self.correlation(r) return self.correlation(np.abs(r) / self.anis[axis - 1])
[docs] def vario_yadrenko(self, zeta): r"""Yadrenko variogram for great-circle distance from latlon-pos.""" return self.variogram(great_circle_to_chordal(zeta, self.geo_scale))
[docs] def cov_yadrenko(self, zeta): r"""Yadrenko covariance for great-circle distance from latlon-pos.""" return self.covariance(great_circle_to_chordal(zeta, self.geo_scale))
[docs] def cor_yadrenko(self, zeta): r"""Yadrenko correlation for great-circle distance from latlon-pos.""" return self.correlation(great_circle_to_chordal(zeta, self.geo_scale))
[docs] def vario_spatial(self, pos): r"""Spatial variogram respecting anisotropy and rotation.""" return self.variogram(self._get_iso_rad(pos))
[docs] def cov_spatial(self, pos): r"""Spatial covariance respecting anisotropy and rotation.""" return self.covariance(self._get_iso_rad(pos))
[docs] def cor_spatial(self, pos): r"""Spatial correlation respecting anisotropy and rotation.""" return self.correlation(self._get_iso_rad(pos))
[docs] def vario_nugget(self, r): """Isotropic variogram of the model respecting the nugget at r=0.""" r = np.asarray(np.abs(r), dtype=np.double) r_gz = np.logical_not(np.isclose(r, 0)) res = np.empty_like(r, dtype=np.double) res[r_gz] = self.variogram(r[r_gz]) res[np.logical_not(r_gz)] = 0.0 return res
[docs] def cov_nugget(self, r): """Isotropic covariance of the model respecting the nugget at r=0.""" r = np.asarray(np.abs(r), dtype=np.double) r_gz = np.logical_not(np.isclose(r, 0)) res = np.empty_like(r, dtype=np.double) res[r_gz] = self.covariance(r[r_gz]) res[np.logical_not(r_gz)] = self.sill return res
[docs] def plot(self, func="variogram", **kwargs): # pragma: no cover """ Plot a function of a the CovModel. Parameters ---------- func : :class:`str`, optional Function to be plotted. Could be one of: * "variogram" * "covariance" * "correlation" * "vario_spatial" * "cov_spatial" * "cor_spatial" * "vario_yadrenko" * "cov_yadrenko" * "cor_yadrenko" * "vario_axis" * "cov_axis" * "cor_axis" * "spectrum" * "spectral_density" * "spectral_rad_pdf" **kwargs Keyword arguments forwarded to the plotting function `"plot_" + func` in :py:mod:`gstools.covmodel.plot`. See Also -------- gstools.covmodel.plot """ routine = getattr(plot, "plot_" + func) return routine(self, **kwargs)
# pykrige functions
[docs] def pykrige_vario(self, args=None, r=0): # pragma: no cover """Isotropic variogram of the model for pykrige.""" if self.latlon: return self.vario_yadrenko(np.deg2rad(r)) return self.variogram(r)
@property def pykrige_anis(self): """2D anisotropy ratio for pykrige.""" if self.dim == 2: return 1 / self.anis[0] return 1.0 # pragma: no cover @property def pykrige_anis_y(self): """3D anisotropy ratio in y direction for pykrige.""" if self.dim >= 2: return 1 / self.anis[0] return 1.0 # pragma: no cover @property def pykrige_anis_z(self): """3D anisotropy ratio in z direction for pykrige.""" if self.dim == 3: return 1 / self.anis[1] return 1.0 # pragma: no cover @property def pykrige_angle(self): """2D rotation angle for pykrige.""" if self.dim == 2: return self.angles[0] / np.pi * 180 return 0.0 # pragma: no cover @property def pykrige_angle_z(self): """3D rotation angle around z for pykrige.""" if self.dim >= 2: return self.angles[0] / np.pi * 180 return 0.0 # pragma: no cover @property def pykrige_angle_y(self): """3D rotation angle around y for pykrige.""" if self.dim == 3: return self.angles[1] / np.pi * 180 return 0.0 # pragma: no cover @property def pykrige_angle_x(self): """3D rotation angle around x for pykrige.""" if self.dim == 3: return self.angles[2] / np.pi * 180 return 0.0 # pragma: no cover @property def pykrige_kwargs(self): """Keyword arguments for pykrige routines.""" kwargs = { "variogram_model": "custom", "variogram_parameters": [], "variogram_function": self.pykrige_vario, } if self.dim == 1: add_kwargs = {} elif self.dim == 2: add_kwargs = { "anisotropy_scaling": self.pykrige_anis, "anisotropy_angle": self.pykrige_angle, } else: add_kwargs = { "anisotropy_scaling_y": self.pykrige_anis_y, "anisotropy_scaling_z": self.pykrige_anis_z, "anisotropy_angle_x": self.pykrige_angle_x, "anisotropy_angle_y": self.pykrige_angle_y, "anisotropy_angle_z": self.pykrige_angle_z, } kwargs.update(add_kwargs) return kwargs # methods for optional/default arguments (can be overridden)
[docs] def default_opt_arg(self): """Provide default optional arguments by the user. Should be given as a dictionary when overridden. """ return { opt: default_arg_from_bounds(bnd) for (opt, bnd) in self.default_opt_arg_bounds().items() }
[docs] def default_opt_arg_bounds(self): """Provide default boundaries for optional arguments.""" res = {} for opt in self.opt_arg: res[opt] = [-np.inf, np.inf] return res
[docs] def check_opt_arg(self): """Run checks for the optional arguments. This is in addition to the bound-checks Notes ----- * You can use this to raise a ValueError/warning * Any return value will be ignored * This method will only be run once, when the class is initialized """
[docs] def check_dim(self, dim): """Check the given dimension.""" return True
[docs] def fix_dim(self): """Set a fix dimension for the model.""" return None
[docs] def var_factor(self): """Factor for the variance.""" return 1.0
[docs] def default_rescale(self): """Provide default rescaling factor.""" return 1.0
# calculation of different scales
[docs] def calc_integral_scale(self): """Calculate the integral scale of the isotrope model.""" self._integral_scale = integral(self.correlation, 0, np.inf)[0] return self._integral_scale
[docs] def percentile_scale(self, per=0.9): """Calculate the percentile scale of the isotrope model. This is the distance, where the given percentile of the variance is reached by the variogram """ return percentile_scale(self, per)
# spectrum methods (can be overridden for speedup)
[docs] def spectrum(self, k): r""" Spectrum of the covariance model. This is given by: .. math:: S(\mathbf{k}) = \left(\frac{1}{2\pi}\right)^n \int C(r) e^{i \mathbf{k}\cdot\mathbf{r}} d^n\mathbf{r} Internally, this is calculated by the hankel transformation: .. math:: S(k) = \left(\frac{1}{2\pi}\right)^n \cdot \frac{(2\pi)^{n/2}}{k^{n/2-1}} \int_0^\infty r^{n/2} C(r) J_{n/2-1}(kr) dr Where :math:`C(r)` is the covariance function of the model. Parameters ---------- k : :class:`float` Radius of the phase: :math:`k=\left\Vert\mathbf{k}\right\Vert` """ return self.spectral_density(k) * self.var
[docs] def spectral_density(self, k): r""" Spectral density of the covariance model. This is given by: .. math:: \tilde{S}(k) = \frac{S(k)}{\sigma^2} Where :math:`S(k)` is the spectrum of the covariance model. Parameters ---------- k : :class:`float` Radius of the phase: :math:`k=\left\Vert\mathbf{k}\right\Vert` """ k = np.asarray(np.abs(k), dtype=np.double) return self._sft.transform(self.correlation, k, ret_err=False)
[docs] def spectral_rad_pdf(self, r): """Radial spectral density of the model.""" return spectral_rad_pdf(self, r)
[docs] def ln_spectral_rad_pdf(self, r): """Log radial spectral density of the model.""" with np.errstate(divide="ignore"): return np.log(self.spectral_rad_pdf(r))
def _has_cdf(self): """State if a cdf is defined with 'spectral_rad_cdf'.""" return hasattr(self, "spectral_rad_cdf") def _has_ppf(self): """State if a ppf is defined with 'spectral_rad_ppf'.""" return hasattr(self, "spectral_rad_ppf") # spatial routines
[docs] 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, radius=self.geo_scale, temporal=self.temporal, time_scale=self.anis[-1], ) return np.dot(matrix_isometrize(self.dim, self.angles, self.anis), pos)
[docs] 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, radius=self.geo_scale, temporal=self.temporal, time_scale=self.anis[-1], ) return np.dot( matrix_anisometrize(self.dim, self.angles, self.anis), pos )
[docs] def main_axes(self): """Axes of the rotated coordinate-system.""" return rotated_main_axes(self.dim, self.angles)
def _get_iso_rad(self, pos): """Isometrized radians.""" pos = np.asarray(pos, dtype=np.double).reshape((self.dim, -1)) iso = np.dot(matrix_isometrize(self.dim, self.angles, self.anis), pos) return np.linalg.norm(iso, axis=0) # fitting routine
[docs] def fit_variogram( self, x_data, y_data, anis=True, sill=None, init_guess="default", weights=None, method="trf", loss="soft_l1", max_eval=None, return_r2=False, curve_fit_kwargs=None, **para_select, ): """ Fitting the variogram-model to an empirical variogram. Parameters ---------- x_data : :class:`numpy.ndarray` The bin-centers of the empirical variogram. y_data : :class:`numpy.ndarray` The measured variogram If multiple are given, they are interpreted as the directional variograms along the main axis of the associated rotated coordinate system. Anisotropy ratios will be estimated in that case. anis : :class:`bool`, optional In case of a directional variogram, you can control anisotropy by this argument. Deselect the parameter from fitting, by setting it "False". You could also pass a fixed value to be set in the model. Then the anisotropy ratios wont be altered during fitting. Default: True sill : :class:`float` or :class:`bool`, optional Here you can provide a fixed sill for the variogram. It needs to be in a fitting range for the var and nugget bounds. If variance or nugget are not selected for estimation, the nugget will be recalculated to fulfill: * sill = var + nugget * if the variance is bigger than the sill, nugget will bet set to its lower bound and the variance will be set to the fitting partial sill. If variance is deselected, it needs to be less than the sill, otherwise a ValueError comes up. Same for nugget. If sill=False, it will be deselected from estimation and set to the current sill of the model. Then, the procedure above is applied. Default: None init_guess : :class:`str` or :class:`dict`, optional Initial guess for the estimation. Either: * "default": using the default values of the covariance model ("len_scale" will be mean of given bin centers; "var" and "nugget" will be mean of given variogram values (if in given bounds)) * "current": using the current values of the covariance model * dict: dictionary with parameter names and given value (separate "default" can bet set to "default" or "current" for unspecified values to get same behavior as given above ("default" by default)) Example: ``{"len_scale": 10, "default": "current"}`` Default: "default" weights : :class:`str`, :class:`numpy.ndarray`, :class:`callable`, optional Weights applied to each point in the estimation. Either: * 'inv': inverse distance ``1 / (x_data + 1)`` * list: weights given per bin * callable: function applied to x_data If callable, it must take a 1-d ndarray. Then ``weights = f(x_data)``. Default: None method : {'trf', 'dogbox'}, optional Algorithm to perform minimization. * 'trf' : Trust Region Reflective algorithm, particularly suitable for large sparse problems with bounds. Generally robust method. * 'dogbox' : dogleg algorithm with rectangular trust regions, typical use case is small problems with bounds. Not recommended for problems with rank-deficient Jacobian. Default: 'trf' loss : :class:`str` or :class:`callable`, optional Determines the loss function in scipys curve_fit. The following keyword values are allowed: * 'linear' (default) : ``rho(z) = z``. Gives a standard least-squares problem. * 'soft_l1' : ``rho(z) = 2 * ((1 + z)**0.5 - 1)``. The smooth approximation of l1 (absolute value) loss. Usually a good choice for robust least squares. * 'huber' : ``rho(z) = z if z <= 1 else 2*z**0.5 - 1``. Works similarly to 'soft_l1'. * 'cauchy' : ``rho(z) = ln(1 + z)``. Severely weakens outliers influence, but may cause difficulties in optimization process. * 'arctan' : ``rho(z) = arctan(z)``. Limits a maximum loss on a single residual, has properties similar to 'cauchy'. If callable, it must take a 1-d ndarray ``z=f**2`` and return an array_like with shape (3, m) where row 0 contains function values, row 1 contains first derivatives and row 2 contains second derivatives. Default: 'soft_l1' max_eval : :class:`int` or :any:`None`, optional Maximum number of function evaluations before the termination. If None (default), the value is chosen automatically: 100 * n. return_r2 : :class:`bool`, optional Whether to return the r2 score of the estimation. Default: False curve_fit_kwargs : :class:`dict`, optional Other keyword arguments passed to scipys curve_fit. Default: None **para_select You can deselect parameters from fitting, by setting them "False" using their names as keywords. You could also pass fixed values for each parameter. Then these values will be applied and the involved parameters wont be fitted. By default, all parameters are fitted. Returns ------- fit_para : :class:`dict` Dictionary with the fitted parameter values pcov : :class:`numpy.ndarray` The estimated covariance of `popt` from :any:`scipy.optimize.curve_fit`. To compute one standard deviation errors on the parameters use ``perr = np.sqrt(np.diag(pcov))``. r2_score : :class:`float`, optional r2 score of the curve fitting results. Only if return_r2 is True. Notes ----- You can set the bounds for each parameter by accessing :any:`CovModel.set_arg_bounds`. The fitted parameters will be instantly set in the model. """ return fit_variogram( model=self, x_data=x_data, y_data=y_data, anis=anis, sill=sill, init_guess=init_guess, weights=weights, method=method, loss=loss, max_eval=max_eval, return_r2=return_r2, curve_fit_kwargs=curve_fit_kwargs, **para_select, )
# bounds setting and checks
[docs] def default_arg_bounds(self): """Provide default boundaries for arguments. Given as a dictionary. """ res = { "var": (0.0, np.inf, "oo"), "len_scale": (0.0, np.inf, "oo"), "nugget": (0.0, np.inf, "co"), "anis": (0.0, np.inf, "oo"), } return res
[docs] def set_arg_bounds(self, check_args=True, **kwargs): r"""Set bounds for the parameters of the model. Parameters ---------- check_args : bool, optional Whether to check if the arguments are in their valid bounds. In case not, a proper default value will be determined. Default: True **kwargs Parameter name as keyword ("var", "len_scale", "nugget", <opt_arg>) and a list of 2 or 3 values: ``[a, b]`` or ``[a, b, <type>]`` where <type> is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ return set_arg_bounds(self, check_args, **kwargs)
[docs] def check_arg_bounds(self): """Check arguments to be within their given bounds.""" return check_arg_bounds(self)
# bounds properties @property def var_bounds(self): """:class:`list`: Bounds for the variance. Notes ----- Is a list of 2 or 3 values: ``[a, b]`` or ``[a, b, <type>]`` where <type> is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ return self._var_bounds @var_bounds.setter def var_bounds(self, bounds): if not check_bounds(bounds): raise ValueError( f"Given bounds for 'var' are not valid, got: {bounds}" ) self._var_bounds = bounds @property def len_scale_bounds(self): """:class:`list`: Bounds for the length scale. Notes ----- Is a list of 2 or 3 values: ``[a, b]`` or ``[a, b, <type>]`` where <type> is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ return self._len_scale_bounds @len_scale_bounds.setter def len_scale_bounds(self, bounds): if not check_bounds(bounds): raise ValueError( f"Given bounds for 'len_scale' are not valid, got: {bounds}" ) self._len_scale_bounds = bounds @property def nugget_bounds(self): """:class:`list`: Bounds for the nugget. Notes ----- Is a list of 2 or 3 values: ``[a, b]`` or ``[a, b, <type>]`` where <type> is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ return self._nugget_bounds @nugget_bounds.setter def nugget_bounds(self, bounds): if not check_bounds(bounds): raise ValueError( f"Given bounds for 'nugget' are not valid, got: {bounds}" ) self._nugget_bounds = bounds @property def anis_bounds(self): """:class:`list`: Bounds for the nugget. Notes ----- Is a list of 2 or 3 values: ``[a, b]`` or ``[a, b, <type>]`` where <type> is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ return self._anis_bounds @anis_bounds.setter def anis_bounds(self, bounds): if not check_bounds(bounds): raise ValueError( f"Given bounds for 'anis' are not valid, got: {bounds}" ) self._anis_bounds = bounds @property def opt_arg_bounds(self): """:class:`dict`: Bounds for the optional arguments. Notes ----- Keys are the opt-arg names and values are lists of 2 or 3 values: ``[a, b]`` or ``[a, b, <type>]`` where <type> is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ return self._opt_arg_bounds @property def arg_bounds(self): """:class:`dict`: Bounds for all parameters. Notes ----- Keys are the arg names and values are lists of 2 or 3 values: ``[a, b]`` or ``[a, b, <type>]`` where <type> is one of ``"oo"``, ``"cc"``, ``"oc"`` or ``"co"`` to define if the bounds are open ("o") or closed ("c"). """ res = { "var": self.var_bounds, "len_scale": self.len_scale_bounds, "nugget": self.nugget_bounds, "anis": self.anis_bounds, } res.update(self.opt_arg_bounds) return res @property def temporal(self): """:class:`bool`: Whether the model is a metric spatio-temporal one.""" return self._temporal # geographical coordinates related @property def latlon(self): """:class:`bool`: Whether the model depends on geographical coords.""" return self._latlon @property def geo_scale(self): """:class:`float`: Geographic scaling for geographical coords.""" return self._geo_scale @property def field_dim(self): """:class:`int`: The (parametric) field dimension of the model (with time).""" return 2 + int(self.temporal) 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.temporal) # standard parameters @property def dim(self): """:class:`int`: The dimension of the model.""" return self._dim @dim.setter def dim(self, dim): set_dim(self, dim) @property def var(self): """:class:`float`: The variance of the model.""" return self._var * self.var_factor() @var.setter def var(self, var): self._var = float(var) / self.var_factor() self.check_arg_bounds() @property def var_raw(self): """:class:`float`: The raw variance of the model without factor. (See. CovModel.var_factor) """ return self._var @var_raw.setter def var_raw(self, var_raw): self._var = float(var_raw) self.check_arg_bounds() @property def nugget(self): """:class:`float`: The nugget of the model.""" return self._nugget @nugget.setter def nugget(self, nugget): self._nugget = float(nugget) self.check_arg_bounds() @property def len_scale(self): """:class:`float`: The main length scale of the model.""" return self._len_scale @len_scale.setter def len_scale(self, len_scale): self._len_scale, anis = set_len_anis( self.dim, len_scale, self.anis, self.latlon ) if self.latlon: self._anis = np.array((self.dim - 1) * [1], dtype=np.double) else: self._anis = anis self.check_arg_bounds() @property def rescale(self): """:class:`float`: Rescale factor for the length scale of the model.""" return self._rescale @rescale.setter def rescale(self, rescale): rescale = self.default_rescale() if rescale is None else rescale self._rescale = abs(float(rescale)) @property def len_rescaled(self): """:class:`float`: The rescaled main length scale of the model.""" return self._len_scale / self._rescale @property def anis(self): """:class:`numpy.ndarray`: The anisotropy factors of the model.""" return self._anis @anis.setter def anis(self, anis): self._len_scale, self._anis = set_len_anis( self.dim, self.len_scale, anis, self.latlon ) self.check_arg_bounds() @property def angles(self): """:class:`numpy.ndarray`: Rotation angles (in rad) of the model.""" return self._angles @angles.setter def angles(self, angles): self._angles = set_model_angles( self.dim, angles, self.latlon, self.temporal ) self.check_arg_bounds() @property def integral_scale(self): """:class:`float`: The main integral scale of the model. Raises ------ ValueError If integral scale is not setable. """ self._integral_scale = self.calc_integral_scale() return self._integral_scale @integral_scale.setter def integral_scale(self, integral_scale): if integral_scale is not None: # format int-scale right self.len_scale = integral_scale integral_scale = self.len_scale # reset len_scale self.len_scale = 1.0 int_tmp = self.calc_integral_scale() self.len_scale = integral_scale / int_tmp if not np.isclose(self.integral_scale, integral_scale, rtol=1e-3): raise ValueError( f"{self.name}: Integral scale could not be set correctly! " "Please just provide a 'len_scale'!" ) @property def hankel_kw(self): """:class:`dict`: :any:`hankel.SymmetricFourierTransform` kwargs.""" return self._hankel_kw @hankel_kw.setter def hankel_kw(self, hankel_kw): if self._hankel_kw is None or hankel_kw is None: self._hankel_kw = copy.copy(HANKEL_DEFAULT) if hankel_kw is not None: self._hankel_kw.update(hankel_kw) if self.dim is not None: self._sft = SFT(ndim=self.dim, **self.hankel_kw) @property def dist_func(self): """:class:`tuple` of :any:`callable`: pdf, cdf and ppf. Spectral distribution info from the model. """ pdf = self.spectral_rad_pdf cdf = None ppf = None if self.has_cdf: cdf = self.spectral_rad_cdf if self.has_ppf: ppf = self.spectral_rad_ppf return pdf, cdf, ppf @property def has_cdf(self): """:class:`bool`: State if a cdf is defined by the user.""" return self._has_cdf() @property def has_ppf(self): """:class:`bool`: State if a ppf is defined by the user.""" return self._has_ppf() @property def sill(self): """:class:`float`: The sill of the variogram. Notes ----- This is calculated by: * ``sill = variance + nugget`` """ return self.var + self.nugget @property def arg(self): """:class:`list` of :class:`str`: Names of all arguments.""" return ["var", "len_scale", "nugget", "anis", "angles"] + self._opt_arg @property def arg_list(self): """:class:`list` of :class:`float`: Values of all arguments.""" alist = [self.var, self.len_scale, self.nugget, self.anis, self.angles] for opt in self.opt_arg: alist.append(getattr(self, opt)) return alist @property def iso_arg(self): """:class:`list` of :class:`str`: Names of isotropic arguments.""" return ["var", "len_scale", "nugget"] + self._opt_arg @property def iso_arg_list(self): """:class:`list` of :class:`float`: Values of isotropic arguments.""" alist = [self.var, self.len_scale, self.nugget] for opt in self.opt_arg: alist.append(getattr(self, opt)) return alist @property def opt_arg(self): """:class:`list` of :class:`str`: Names of the optional arguments.""" return self._opt_arg @property def len_scale_vec(self): """:class:`numpy.ndarray`: The length scales in each direction. Notes ----- This is calculated by: * ``len_scale_vec[0] = len_scale`` * ``len_scale_vec[1] = len_scale*anis[0]`` * ``len_scale_vec[2] = len_scale*anis[1]`` """ res = np.zeros(self.dim, dtype=np.double) res[0] = self.len_scale for i in range(1, self._dim): res[i] = self.len_scale * self.anis[i - 1] return res @property def integral_scale_vec(self): """:class:`numpy.ndarray`: The integral scales in each direction. Notes ----- This is calculated by: * ``integral_scale_vec[0] = integral_scale`` * ``integral_scale_vec[1] = integral_scale*anis[0]`` * ``integral_scale_vec[2] = integral_scale*anis[1]`` """ res = np.zeros(self.dim, dtype=np.double) res[0] = self.integral_scale for i in range(1, self.dim): res[i] = self.integral_scale * self.anis[i - 1] return res @property def name(self): """:class:`str`: The name of the CovModel class.""" return self.__class__.__name__ @property def do_rotation(self): """:any:`bool`: State if a rotation is performed.""" return not np.all(np.isclose(self.angles, 0.0)) @property def is_isotropic(self): """:any:`bool`: State if a model is isotropic.""" return np.all(np.isclose(self.anis, 1.0)) def __eq__(self, other): """Compare CovModels.""" if not isinstance(other, CovModel): return False return compare(self, other) def __setattr__(self, name, value): """Set an attribute.""" super().__setattr__(name, value) # if an optional variogram argument was given, check bounds if hasattr(self, "_opt_arg") and name in self._opt_arg: self.check_arg_bounds() def __repr__(self): """Return String representation.""" return model_repr(self)