diff --git a/CHANGELOG.md b/CHANGELOG.md index 61aa289f..733d75e3 100755 --- 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 +- `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,7 +325,7 @@ 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 @@ -305,23 +333,23 @@ See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197) ### 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 diff --git a/README.md b/README.md index 2740aa9f..69f1f3ec 100644 --- 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 index 86ec0671..df583778 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -200,7 +200,7 @@ This works perfectly well with `cartopy `_). 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. @@ -37,20 +37,21 @@ 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 `_. +`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 `_ -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 index 4671b2f7..4b4c6b8a 100644 --- 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 index 049225d3..81c78964 100644 --- 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 new file mode 100644 index 00000000..56684a94 --- /dev/null +++ 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 index 07aa5faa..3cb06b9e 100644 --- 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 index 41090458..3a6c7ec0 100644 --- 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 index 60dd1e33..6d62d558 100644 --- 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", diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 686768fa..be530075 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -28,16 +28,18 @@ 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, - set_angles, ) __all__ = ["CovModel"] @@ -52,7 +54,10 @@ class CovModel: Parameters ---------- dim : :class:`int`, optional - dimension of the model. Default: ``3`` + 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`` @@ -100,9 +105,27 @@ class CovModel: :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. `rescale` can be set to e.g. earth's radius, + 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. @@ -128,9 +151,13 @@ def __init__( 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, @@ -156,11 +183,16 @@ def __init__( self._nugget_bounds = None self._anis_bounds = None self._opt_arg_bounds = {} - # Set latlon first + # 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 - self.dim = dim + # 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) @@ -173,14 +205,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: - self._anis = np.array((self.dim - 1) * [1], dtype=np.double) - 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._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 @@ -242,15 +275,15 @@ def cor_axis(self, r, axis=0): def vario_yadrenko(self, zeta): r"""Yadrenko variogram for great-circle distance from latlon-pos.""" - return self.variogram(2 * np.sin(zeta / 2)) + return self.variogram(great_circle_to_chordal(zeta, self.geo_scale)) def cov_yadrenko(self, zeta): r"""Yadrenko covariance for great-circle distance from latlon-pos.""" - return self.covariance(2 * np.sin(zeta / 2)) + return self.covariance(great_circle_to_chordal(zeta, self.geo_scale)) def cor_yadrenko(self, zeta): r"""Yadrenko correlation for great-circle distance from latlon-pos.""" - return self.correlation(2 * np.sin(zeta / 2)) + return self.correlation(great_circle_to_chordal(zeta, self.geo_scale)) def vario_spatial(self, pos): r"""Spatial variogram respecting anisotropy and rotation.""" @@ -530,14 +563,24 @@ 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, + radius=self.geo_scale, + temporal=self.temporal, + time_scale=self.anis[-1], + ) 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, + radius=self.geo_scale, + temporal=self.temporal, + time_scale=self.anis[-1], + ) return np.dot( matrix_anisometrize(self.dim, self.angles, self.anis), pos ) @@ -548,7 +591,9 @@ def main_axes(self): def _get_iso_rad(self, pos): """Isometrized radians.""" - return np.linalg.norm(self.isometrize(pos), axis=0) + 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 @@ -862,6 +907,11 @@ def arg_bounds(self): 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 @@ -869,10 +919,20 @@ 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 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.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 @@ -925,7 +985,9 @@ def len_scale(self): @len_scale.setter def len_scale(self, len_scale): - self._len_scale, anis = set_len_anis(self.dim, len_scale, self.anis) + 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: @@ -954,12 +1016,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.check_arg_bounds() @property @@ -969,10 +1028,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 diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index 2ff5398b..dc2d5a3a 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -13,7 +13,7 @@ from scipy.optimize import curve_fit from gstools.covmodel.tools import check_arg_in_bounds, default_arg_from_bounds -from gstools.tools.geometric import set_anis +from gstools.tools.geometric import great_circle_to_chordal, set_anis __all__ = ["fit_variogram"] @@ -341,7 +341,7 @@ def _check_vario(model, x_data, y_data): ) if model.latlon: # convert to yadrenko model - x_data = 2 * np.sin(x_data / 2) + x_data = great_circle_to_chordal(x_data, model.geo_scale) return x_data, y_data, is_dir_vario @@ -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.""" diff --git a/src/gstools/covmodel/plot.py b/src/gstools/covmodel/plot.py index efcc2630..334ccbe2 100644 --- a/src/gstools/covmodel/plot.py +++ b/src/gstools/covmodel/plot.py @@ -52,12 +52,14 @@ # plotting routines ####################################################### -def _plot_spatial(dim, pos, field, fig, ax, latlon, **kwargs): +def _plot_spatial(dim, pos, field, fig, ax, temporal, **kwargs): from gstools.field.plot import plot_1d, plot_nd if dim == 1: - return plot_1d(pos, field, fig, ax, **kwargs) - return plot_nd(pos, field, "structured", fig, ax, latlon, **kwargs) + return plot_1d(pos, field, fig, ax, temporal, **kwargs) + return plot_nd( + pos, field, "structured", fig, ax, temporal=temporal, **kwargs + ) def plot_vario_spatial( @@ -70,7 +72,9 @@ def plot_vario_spatial( pos = [x_s] * model.dim shp = tuple(len(p) for p in pos) fld = model.vario_spatial(generate_grid(pos)).reshape(shp) - return _plot_spatial(model.dim, pos, fld, fig, ax, model.latlon, **kwargs) + return _plot_spatial( + model.dim, pos, fld, fig, ax, model.temporal, **kwargs + ) def plot_cov_spatial( @@ -83,7 +87,9 @@ def plot_cov_spatial( pos = [x_s] * model.dim shp = tuple(len(p) for p in pos) fld = model.cov_spatial(generate_grid(pos)).reshape(shp) - return _plot_spatial(model.dim, pos, fld, fig, ax, model.latlon, **kwargs) + return _plot_spatial( + model.dim, pos, fld, fig, ax, model.temporal, **kwargs + ) def plot_cor_spatial( @@ -96,7 +102,9 @@ def plot_cor_spatial( pos = [x_s] * model.dim shp = tuple(len(p) for p in pos) fld = model.cor_spatial(generate_grid(pos)).reshape(shp) - return _plot_spatial(model.dim, pos, fld, fig, ax, model.latlon, **kwargs) + return _plot_spatial( + model.dim, pos, fld, fig, ax, model.temporal, **kwargs + ) def plot_variogram( @@ -150,7 +158,7 @@ def plot_vario_yadrenko( """Plot Yadrenko variogram of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = min(3 * model.len_rescaled, np.pi) + x_max = min(3 * model.len_scale, model.geo_scale * np.pi) x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} Yadrenko variogram") ax.plot(x_s, model.vario_yadrenko(x_s), **kwargs) @@ -165,7 +173,7 @@ def plot_cov_yadrenko( """Plot Yadrenko covariance of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = min(3 * model.len_rescaled, np.pi) + x_max = min(3 * model.len_scale, model.geo_scale * np.pi) x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} Yadrenko covariance") ax.plot(x_s, model.cov_yadrenko(x_s), **kwargs) @@ -180,7 +188,7 @@ def plot_cor_yadrenko( """Plot Yadrenko correlation function of a given CovModel.""" fig, ax = get_fig_ax(fig, ax) if x_max is None: - x_max = min(3 * model.len_rescaled, np.pi) + x_max = min(3 * model.len_scale, model.geo_scale * np.pi) x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} Yadrenko correlation") ax.plot(x_s, model.cor_yadrenko(x_s), **kwargs) diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 98ed3b8a..dddeb441 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -30,7 +30,7 @@ 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__ = [ @@ -38,6 +38,7 @@ "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): """Set the length scale and anisotropy factors for the given dimension. Parameters @@ -194,6 +195,10 @@ 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. In this case there is no spatial anisotropy. + Default: False Returns ------- @@ -233,11 +238,49 @@ def set_len_anis(dim, len_scale, anis): for ani in out_anis: if not ani > 0.0: raise ValueError( - "anisotropy-ratios needs to be > 0, " + "got: " + str(out_anis) + f"anisotropy-ratios needs to be > 0, got: {out_anis}" ) + # no spatial anisotropy for latlon + if latlon: + out_anis[:2] = 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. @@ -378,9 +421,7 @@ def percentile_scale(model, per=0.9): """ # check the given percentile if not 0.0 < per < 1.0: - raise ValueError( - "percentile needs to be within (0, 1), got: " + str(per) - ) + raise ValueError(f"percentile needs to be within (0, 1), got: {per}") # define a curve, that has its root at the wanted point def curve(x): @@ -494,17 +535,17 @@ def set_dim(model, dim): # check if a fixed dimension should be used if model.fix_dim() is not None and model.fix_dim() != dim: warnings.warn( - model.name + ": using fixed dimension " + str(model.fix_dim()), + f"{model.name}: using fixed dimension {model.fix_dim()}", AttributeWarning, ) dim = model.fix_dim() - if model.latlon and dim != 3: + if model.latlon and dim != (3 + int(model.temporal)): 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 temporal={model.temporal})." ) - # force dim=3 for latlon models - dim = 3 if model.latlon else dim + # force dim=3 (or 4 when temporal=True) for latlon models + dim = (3 + int(model.temporal)) if model.latlon else dim # set the dimension if dim < 1: raise ValueError("Only dimensions of d >= 1 are supported.") @@ -522,7 +563,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() @@ -551,6 +594,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.temporal == that.temporal for opt in this.opt_arg: equal &= np.isclose(getattr(this, opt), getattr(that, opt)) return equal @@ -568,21 +612,35 @@ def model_repr(model): # pragma: no cover m = model p = model._prec opt_str = "" + t_str = ", temporal=True" if m.temporal else "" if not np.isclose(m.rescale, m.default_rescale()): opt_str += f", rescale={m.rescale:.{p}}" for opt in m.opt_arg: opt_str += f", {opt}={getattr(m, opt):.{p}}" - # only print anis and angles if model is anisotropic or rotated - ani_str = "" if m.is_isotropic else f", anis={list_format(m.anis, p)}" - ang_str = f", angles={list_format(m.angles, p)}" if m.do_rotation else "" if m.latlon: + ani_str = ( + "" + if m.is_isotropic or not m.temporal + else f", anis={m.anis[-1]:.{p}}" + ) + r_str = ( + "" + if np.isclose(m.geo_scale, 1) + else f", geo_scale={m.geo_scale:.{p}}" + ) repr_str = ( - f"{m.name}(latlon={m.latlon}, var={m.var:.{p}}, " - f"len_scale={m.len_scale:.{p}}, nugget={m.nugget:.{p}}{opt_str})" + f"{m.name}(latlon={m.latlon}{t_str}, var={m.var:.{p}}, " + f"len_scale={m.len_scale:.{p}}, nugget={m.nugget:.{p}}" + f"{ani_str}{r_str}{opt_str})" ) else: + # only print anis and angles if model is anisotropic or rotated + ani_str = "" if m.is_isotropic else f", anis={list_format(m.anis, p)}" + ang_str = ( + f", angles={list_format(m.angles, p)}" if m.do_rotation else "" + ) repr_str = ( - f"{m.name}(dim={m.dim}, var={m.var:.{p}}, " + f"{m.name}(dim={m.spatial_dim}{t_str}, var={m.var:.{p}}, " f"len_scale={m.len_scale:.{p}}, nugget={m.nugget:.{p}}" f"{ani_str}{ang_str}{opt_str})" ) diff --git a/src/gstools/field/base.py b/src/gstools/field/base.py index 903f3893..bb514141 100755 --- a/src/gstools/field/base.py +++ b/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 temporal(self): + """:class:`bool`: Whether the field depends on time.""" + return False if self.model is None else self.model.temporal + @property def name(self): """:class:`str`: The name of the class.""" diff --git a/src/gstools/field/plot.py b/src/gstools/field/plot.py index 528cdcc3..d06a22a1 100644 --- a/src/gstools/field/plot.py +++ b/src/gstools/field/plot.py @@ -52,13 +52,22 @@ def plot_field( Forwarded to the plotting routine. """ if fld.dim == 1: - return plot_1d(fld.pos, fld[field], fig, ax, **kwargs) + return plot_1d(fld.pos, fld[field], fig, ax, fld.temporal, **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.temporal, + **kwargs, ) -def plot_1d(pos, field, fig=None, ax=None, ax_names=None): # pragma: no cover +def plot_1d( + pos, field, fig=None, ax=None, temporal=False, ax_names=None +): # pragma: no cover """ Plot a 1D field. @@ -69,6 +78,11 @@ def plot_1d(pos, field, fig=None, ax=None, ax_names=None): # pragma: no cover or the axes descriptions (for mesh_type='structured') field : :class:`numpy.ndarray` Field values. + temporal : :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 fig : :class:`Figure` or :any:`None`, optional Figure to plot the axes on. If `None`, a new one will be created. Default: `None` @@ -88,7 +102,7 @@ def plot_1d(pos, field, fig=None, ax=None, ax_names=None): # pragma: no cover x = pos[0] x = x.flatten() arg = np.argsort(x) - ax_names = _ax_names(1, ax_names=ax_names) + ax_names = _ax_names(1, temporal=temporal, ax_names=ax_names) ax.plot(x[arg], field.ravel()[arg]) ax.set_xlabel(ax_names[0]) ax.set_ylabel(ax_names[1]) @@ -104,6 +118,7 @@ def plot_nd( fig=None, ax=None, latlon=False, + temporal=False, resolution=128, ax_names=None, aspect="quad", @@ -136,6 +151,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 + temporal : :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 +179,28 @@ def plot_nd( """ dim = len(pos) assert dim > 1 - assert not latlon or dim == 2 + assert not latlon or dim == 2 + int(bool(temporal)) if dim == 2 and contour_plot: return _plot_2d( - pos, field, mesh_type, fig, ax, latlon, ax_names, **kwargs + pos, + field, + mesh_type, + fig, + ax, + latlon, + temporal, + 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 temporal: + 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, temporal, 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 +357,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, temporal=False, ax_names=None): + t_fac = int(bool(temporal)) 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 <= 3: - return ["$x$", "$y$", "$z$"][:dim] + (dim == 1) * ["field"] - return [f"$x_{{{i}}}$" for i in range(dim)] + if dim == 2 + t_fac and latlon: + return ["lon", "lat"] + t_fac * ["time"] + if dim - t_fac <= 3: + return ( + ["$x$", "$y$", "$z$"][: dim - t_fac] + + t_fac * ["time"] + + (dim == 1) * ["field"] + ) + return [f"$x_{{{i}}}$" for i in range(dim - t_fac)] + t_fac * ["time"] def _plot_2d( @@ -341,6 +380,7 @@ def _plot_2d( fig=None, ax=None, latlon=False, + temporal=False, ax_names=None, levels=64, antialias=True, @@ -348,7 +388,7 @@ def _plot_2d( """Plot a 2d field with a contour plot.""" fig, ax = get_fig_ax(fig, ax) title = f"Field 2D {mesh_type}: {field.shape}" - ax_names = _ax_names(2, latlon, ax_names=ax_names) + ax_names = _ax_names(2, latlon, temporal, ax_names=ax_names) x, y = pos[::-1] if latlon else pos if mesh_type == "unstructured": cont = ax.tricontourf(x, y, field.ravel(), levels=levels) diff --git a/src/gstools/krige/base.py b/src/gstools/krige/base.py index 77c6832e..336f4cc0 100755 --- a/src/gstools/krige/base.py +++ b/src/gstools/krige/base.py @@ -113,8 +113,12 @@ class Krige(Field): Default: False fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False @@ -496,8 +500,12 @@ def set_condition( Default: False fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -522,12 +530,18 @@ def set_condition( self.normalizer.fit(self.cond_val - self.cond_trend) if fit_variogram: # fitting model to empirical variogram of data # normalize field + if self.model.latlon and self.model.temporal: + msg = "Krige: can't fit variogram for spatio-temporal latlon data." + raise ValueError(msg) field = self.normalizer.normalize(self.cond_val - self.cond_trend) field -= self.cond_mean sill = np.var(field) if self.model.is_isotropic: emp_vario = vario_estimate( - self.cond_pos, field, latlon=self.model.latlon + self.cond_pos, + field, + latlon=self.model.latlon, + geo_scale=self.model.geo_scale, ) else: axes = rotated_main_axes(self.model.dim, self.model.angles) diff --git a/src/gstools/krige/methods.py b/src/gstools/krige/methods.py index 653785f8..b258a02d 100644 --- a/src/gstools/krige/methods.py +++ b/src/gstools/krige/methods.py @@ -77,8 +77,12 @@ class Simple(Krige): Default: False fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -171,8 +175,12 @@ class Ordinary(Krige): Default: False fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -275,8 +283,12 @@ class Universal(Krige): Default: False fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -376,8 +388,12 @@ class ExtDrift(Krige): Default: False fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -467,8 +483,12 @@ class Detrended(Krige): Default: `"pinv"` fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ diff --git a/src/gstools/tools/__init__.py b/src/gstools/tools/__init__.py index bd329576..1f68dbaf 100644 --- a/src/gstools/tools/__init__.py +++ b/src/gstools/tools/__init__.py @@ -58,10 +58,19 @@ .. autosummary:: EARTH_RADIUS + KM_SCALE + DEGREE_SCALE + RADIAN_SCALE ---- .. autodata:: EARTH_RADIUS + +.. autodata:: KM_SCALE + +.. autodata:: DEGREE_SCALE + +.. autodata:: RADIAN_SCALE """ from gstools.tools.export import ( @@ -103,6 +112,15 @@ EARTH_RADIUS = 6371.0 """float: earth radius for WGS84 ellipsoid in km""" +KM_SCALE = 6371.0 +"""float: earth radius for WGS84 ellipsoid in km""" + +DEGREE_SCALE = 57.29577951308232 +"""float: radius for unit sphere in degree""" + +RADIAN_SCALE = 1.0 +"""float: radius for unit sphere""" + __all__ = [ "vtk_export", @@ -135,4 +153,7 @@ "generate_grid", "generate_st_grid", "EARTH_RADIUS", + "KM_SCALE", + "DEGREE_SCALE", + "RADIAN_SCALE", ] diff --git a/src/gstools/tools/geometric.py b/src/gstools/tools/geometric.py index afdcacaf..7f2ea10e 100644 --- a/src/gstools/tools/geometric.py +++ b/src/gstools/tools/geometric.py @@ -27,6 +27,7 @@ latlon2pos pos2latlon chordal_to_great_circle + great_circle_to_chordal """ # pylint: disable=C0103 import numpy as np @@ -624,71 +625,98 @@ 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, temporal=False, time_scale=1.0 +): """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` + Sphere 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 + temporal : :class:`bool`, optional + Whether latlon includes an appended time axis. + Default: False + time_scale : :class:`float`, optional + Scaling factor (e.g. anisotropy) for the time axis. + Default: `1.0` Returns ------- :class:`numpy.ndarray` the 3D position array """ - latlon = np.asarray(latlon, dtype=dtype).reshape((2, -1)) - 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, + latlon = np.asarray(latlon, dtype=dtype).reshape( + (3 if temporal else 2, -1) ) + lat, lon = np.deg2rad(latlon[:2]) + 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 temporal: + return np.array(pos_tuple + (latlon[2] / time_scale,), 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, temporal=False, time_scale=1.0 +): """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` + Sphere 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 + temporal : :class:`bool`, optional + Whether latlon includes an appended time axis. + Default: False + time_scale : :class:`float`, optional + Scaling factor (e.g. anisotropy) for the time axis. + Default: `1.0` 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 temporal 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]) - return np.rad2deg((lat, lon), dtype=dtype) + latlon = np.rad2deg((lat, lon), dtype=dtype) + if temporal: + return np.array( + (latlon[0], latlon[1], pos[3] * time_scale), dtype=dtype + ) + return latlon -def chordal_to_great_circle(dist): +def chordal_to_great_circle(dist, radius=1.0): """ Calculate great circle distance corresponding to given chordal distance. Parameters ---------- dist : array_like - Chordal distance of two points on the unit-sphere. + Chordal distance of two points on the sphere. + radius : :class:`float`, optional + Sphere radius. Default: `1.0` Returns ------- @@ -697,6 +725,29 @@ def chordal_to_great_circle(dist): Notes ----- - If given values are not in [0, 1], they will be truncated. + If given values are not in [0, 2 * radius], they will be truncated. + """ + diameter = 2 * radius + return diameter * np.arcsin( + np.maximum(np.minimum(np.divide(dist, diameter), 1), 0) + ) + + +def great_circle_to_chordal(dist, radius=1.0): + """ + Calculate chordal distance corresponding to given great circle distance. + + Parameters + ---------- + dist : array_like + Great circle distance of two points on the sphere. + radius : :class:`float`, optional + Sphere radius. Default: `1.0` + + Returns + ------- + :class:`numpy.ndarray` + Chordal distance corresponding to given great circle distance. """ - return 2 * np.arcsin(np.maximum(np.minimum(np.divide(dist, 2), 1), 0)) + diameter = 2 * radius + return diameter * np.sin(np.divide(dist, diameter)) diff --git a/src/gstools/transform/field.py b/src/gstools/transform/field.py index 9ac33b6c..81824739 100644 --- a/src/gstools/transform/field.py +++ b/src/gstools/transform/field.py @@ -26,7 +26,7 @@ normal_to_arcsin normal_to_uquad """ -# pylint: disable=C0103, C0123, R0911 +# pylint: disable=C0103, C0123, R0911, R1735 import numpy as np from gstools.normalizer import ( diff --git a/src/gstools/variogram/binning.py b/src/gstools/variogram/binning.py index be490110..e8e42f38 100644 --- a/src/gstools/variogram/binning.py +++ b/src/gstools/variogram/binning.py @@ -10,6 +10,7 @@ """ import numpy as np +from gstools.tools import RADIAN_SCALE from gstools.tools.geometric import ( chordal_to_great_circle, format_struct_pos_dim, @@ -31,6 +32,7 @@ def standard_bins( mesh_type="unstructured", bin_no=None, max_dist=None, + geo_scale=RADIAN_SCALE, ): r""" Get standard binning. @@ -62,6 +64,13 @@ def standard_bins( Cut of length for the bins. If None is given, it will be set to one third of the box-diameter from the given points. Default: None + geo_scale : :class:`float`, optional + Geographic unit scaling in case of latlon coordinates to get a + meaningful bins unit. + By default, bins are assumed to be in radians with latlon=True. + Can be set to :any:`KM_SCALE` to have bins in km or + :any:`DEGREE_SCALE` to have bins in degrees. + Default: :any:`RADIAN_SCALE` Returns ------- @@ -80,7 +89,7 @@ def standard_bins( pos = generate_grid(format_struct_pos_dim(pos, dim)[0]) else: pos = np.asarray(pos, dtype=np.double).reshape(dim, -1) - pos = latlon2pos(pos) if latlon else pos + pos = latlon2pos(pos, radius=geo_scale) if latlon else pos pnt_cnt = len(pos[0]) box = [] for axis in pos: @@ -88,7 +97,7 @@ def standard_bins( box = np.asarray(box) diam = np.linalg.norm(box[:, 0] - box[:, 1]) # convert diameter to great-circle distance if using latlon - diam = chordal_to_great_circle(diam) if latlon else diam + diam = chordal_to_great_circle(diam, geo_scale) if latlon else diam bin_no = _sturges(pnt_cnt) if bin_no is None else int(bin_no) max_dist = diam / 3 if max_dist is None else float(max_dist) return np.linspace(0, max_dist, num=bin_no + 1, dtype=np.double) diff --git a/src/gstools/variogram/variogram.py b/src/gstools/variogram/variogram.py index b57e0354..69746658 100644 --- a/src/gstools/variogram/variogram.py +++ b/src/gstools/variogram/variogram.py @@ -14,6 +14,7 @@ from gstools import config from gstools.normalizer.tools import remove_trend_norm_mean +from gstools.tools import RADIAN_SCALE from gstools.tools.geometric import ( ang2dir, format_struct_pos_shape, @@ -92,6 +93,8 @@ def vario_estimate( normalizer=None, trend=None, fit_normalizer=False, + geo_scale=RADIAN_SCALE, + **std_bins, ): r""" Estimates the empirical variogram. @@ -222,10 +225,20 @@ def vario_estimate( fit_normalizer : :class:`bool`, optional Whether to fit the data-normalizer to the given (detrended) field. Default: False + geo_scale : :class:`float`, optional + Geographic unit scaling in case of latlon coordinates to get a + meaningful bins unit. + By default, bins are assumed to be in radians with latlon=True. + Can be set to :any:`KM_SCALE` to have bins in km or + :any:`DEGREE_SCALE` to have bins in degrees. + Default: :any:`RADIAN_SCALE` + **std_bins + Optional arguments that are forwarded to the :any:`standard_bins` routine + if no bins are given (bin_no, max_dist). Returns ------- - bin_center : (n), :class:`numpy.ndarray` + bin_centers : (n), :class:`numpy.ndarray` The bin centers. gamma : (n) or (d, n), :class:`numpy.ndarray` The estimated variogram values at bin centers. @@ -250,18 +263,18 @@ def vario_estimate( """ if bin_edges is not None: bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double, copy=False) - bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0 # allow multiple fields at same positions (ndmin=2: first axis -> field ID) # need to convert to ma.array, since list of ma.array is not recognised field = np.ma.array(field, ndmin=2, dtype=np.double, copy=True) masked = np.ma.is_masked(field) or np.any(mask) # catch special case if everything is masked if masked and np.all(mask): - bin_centres = np.empty(0) if bin_edges is None else bin_centres - estimates = np.zeros_like(bin_centres) + bin_centers = np.empty(0) if bin_edges is None else bin_centers + estimates = np.zeros_like(bin_centers) if return_counts: - return bin_centres, estimates, np.zeros_like(estimates, dtype=int) - return bin_centres, estimates + return bin_centers, estimates, np.zeros_like(estimates, dtype=int) + return bin_centers, estimates if not masked: field = field.filled() # check mesh shape @@ -331,10 +344,15 @@ def vario_estimate( ) field = field[:, sampled_idx] pos = pos[:, sampled_idx] - # create bining if not given + # create bins if bin_edges is None: - bin_edges = standard_bins(pos, dim, latlon) - bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + bin_edges = standard_bins( + pos, dim, latlon, geo_scale=geo_scale, **std_bins + ) + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + if latlon: + # internally we always use radians + bin_edges /= geo_scale # normalize field norm_field_out = remove_trend_norm_mean( *(pos, field, mean, normalizer, trend), @@ -371,7 +389,7 @@ def vario_estimate( if dir_no == 1: estimates, counts = estimates[0], counts[0] est_out = (estimates, counts) - return (bin_centres,) + est_out[: 2 if return_counts else 1] + norm_out + return (bin_centers,) + est_out[: 2 if return_counts else 1] + norm_out def vario_estimate_axis( diff --git a/tests/test_krige.py b/tests/test_krige.py index a37bf1e1..d702b0ee 100644 --- a/tests/test_krige.py +++ b/tests/test_krige.py @@ -132,7 +132,6 @@ def test_universal(self): ) def test_detrended(self): - for Model in self.cov_models: for dim in self.dims: model = Model( diff --git a/tests/test_latlon.py b/tests/test_latlon.py index 10d9f682..98088db8 100644 --- a/tests/test_latlon.py +++ b/tests/test_latlon.py @@ -24,7 +24,7 @@ def fix_dim(self): class TestLatLon(unittest.TestCase): def setUp(self): self.cmod = gs.Gaussian( - latlon=True, var=2, len_scale=777, rescale=gs.EARTH_RADIUS + latlon=True, var=2, len_scale=777, geo_scale=gs.KM_SCALE ) self.lat = self.lon = range(-80, 81) @@ -66,7 +66,10 @@ def test_conv(self): 6, self.cmod.anisometrize(self.cmod.isometrize((8, 6)))[1, 0] ) self.assertAlmostEqual( - 1, self.cmod.isometrize(self.cmod.anisometrize((1, 0, 0)))[0, 0] + gs.EARTH_RADIUS, + self.cmod.isometrize( + self.cmod.anisometrize((gs.EARTH_RADIUS, 0, 0)) + )[0, 0], ) def test_cov_model(self): @@ -91,15 +94,16 @@ def test_vario_est(self): srf = gs.SRF(self.cmod, seed=12345) field = srf.structured((self.lat, self.lon)) - bin_edges = [0.01 * i for i in range(30)] + bin_edges = np.linspace(0, 3 * 777, 30) bin_center, emp_vario = gs.vario_estimate( *((self.lat, self.lon), field, bin_edges), latlon=True, mesh_type="structured", sampling_size=2000, sampling_seed=12345, + geo_scale=gs.KM_SCALE, ) - mod = gs.Gaussian(latlon=True, rescale=gs.EARTH_RADIUS) + mod = gs.Gaussian(latlon=True, geo_scale=gs.KM_SCALE) mod.fit_variogram(bin_center, emp_vario, nugget=False) # allow 10 percent relative error self.assertLess(_rel_err(mod.var, self.cmod.var), 0.1) @@ -114,7 +118,7 @@ def test_krige(self): bin_edges, latlon=True, ) - mod = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS) + mod = gs.Spherical(latlon=True, geo_scale=gs.KM_SCALE) mod.fit_variogram(*emp_vario, nugget=False) kri = gs.krige.Ordinary( mod, @@ -134,7 +138,7 @@ def test_cond_srf(self): bin_edges, latlon=True, ) - mod = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS) + mod = gs.Spherical(latlon=True, geo_scale=gs.KM_SCALE) mod.fit_variogram(*emp_vario, nugget=False) krige = gs.krige.Ordinary( mod, (self.data[:, 0], self.data[:, 1]), self.data[:, 2] @@ -144,7 +148,7 @@ def test_cond_srf(self): for i, dat in enumerate(self.data[:, 2]): self.assertAlmostEqual(field[i], dat, 3) - def error_test(self): + def test_error(self): # try fitting directional variogram mod = gs.Gaussian(latlon=True) with self.assertRaises(ValueError): diff --git a/tests/test_temporal.py b/tests/test_temporal.py new file mode 100644 index 00000000..c179db1c --- /dev/null +++ b/tests/test_temporal.py @@ -0,0 +1,78 @@ +""" +This is the unittest for temporal related routines. +""" + +import unittest + +import numpy as np + +import gstools as gs + + +class TestTemporal(unittest.TestCase): + def setUp(self): + self.mod = gs.Gaussian( + latlon=True, + temporal=True, + len_scale=1000, + anis=0.5, + geo_scale=gs.KM_SCALE, + ) + + def test_latlon(self): + mod = gs.Gaussian( + latlon=True, temporal=True, angles=[1, 2, 3, 4, 5, 6] + ) + self.assertEqual(mod.dim, 4) + self.assertEqual(mod.field_dim, 3) + self.assertEqual(mod.spatial_dim, 2) + self.assertTrue(np.allclose(mod.angles, 0)) + + mod1 = gs.Gaussian(latlon=True, temporal=True, len_scale=[10, 5]) + mod2 = gs.Gaussian(latlon=True, temporal=True, len_scale=10, anis=0.5) + + self.assertTrue(np.allclose(mod1.anis, mod2.anis)) + self.assertAlmostEqual(mod1.len_scale, mod2.len_scale) + + def test_latlon2pos(self): + self.assertAlmostEqual( + 8, self.mod.anisometrize(self.mod.isometrize((8, 6, 9)))[0, 0] + ) + self.assertAlmostEqual( + 6, self.mod.anisometrize(self.mod.isometrize((8, 6, 9)))[1, 0] + ) + self.assertAlmostEqual( + 9, self.mod.anisometrize(self.mod.isometrize((8, 6, 9)))[2, 0] + ) + self.assertAlmostEqual( + gs.EARTH_RADIUS, + self.mod.isometrize( + self.mod.anisometrize((gs.EARTH_RADIUS, 0, 0, 10)) + )[0, 0], + ) + self.assertAlmostEqual( + 10, + self.mod.isometrize( + self.mod.anisometrize((gs.EARTH_RADIUS, 0, 0, 10)) + )[3, 0], + ) + + def test_rotation(self): + mod = gs.Gaussian( + spatial_dim=3, temporal=True, angles=[1, 2, 3, 4, 5, 6] + ) + self.assertTrue(np.allclose(mod.angles, [1, 2, 3, 0, 0, 0])) + self.assertEqual(mod.dim, 4) + + def test_krige(self): + # auto-fitting latlon-temporal model in kriging not possible + with self.assertRaises(ValueError): + kri = gs.Krige(self.mod, 3 * [[1, 2]], [1, 2], fit_variogram=True) + + def test_field(self): + srf = gs.SRF(self.mod) + self.assertTrue(srf.temporal) + + +if __name__ == "__main__": + unittest.main()