From 7b868629846333a5a0289b0972d4a99d2b3d652f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 6 Jun 2023 12:15:46 +0200 Subject: [PATCH 01/40] Covmodel: all kwargs after dim are now keyword only --- src/gstools/covmodel/base.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 686768fa..5a1a1f5d 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -123,6 +123,7 @@ class CovModel: def __init__( self, dim=3, + *, var=1.0, len_scale=1.0, nugget=0.0, From 985f62b8a62f37f84f8455deddde670448487910 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 6 Jun 2023 14:21:17 +0200 Subject: [PATCH 02/40] tests: minimal black fixes --- tests/test_krige.py | 1 - 1 file changed, 1 deletion(-) 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( From 83bcf137931758ab6e5c934297ea0dc134a3cbb3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 6 Jun 2023 14:21:29 +0200 Subject: [PATCH 03/40] covmodel: add time attribute --- src/gstools/covmodel/base.py | 33 ++++++++++++++++++++++++------ src/gstools/covmodel/fit.py | 1 + src/gstools/covmodel/tools.py | 9 +++++---- src/gstools/tools/geometric.py | 37 ++++++++++++++++++++++++---------- 4 files changed, 59 insertions(+), 21 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 5a1a1f5d..3a923c1c 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -103,6 +103,12 @@ class CovModel: disabled. `rescale` can be set to e.g. earth's radius, to have a meaningful `len_scale` parameter. Default: False + time : :class:`bool`, optional + Create a metric spatio-temporal covariance model. + Setting this to true will increase `dim` and `field_dim` by 1. + `spatial_dim` will be `field_dim - 1`. + The time-dimension is appended, meaning the pos tuple is (x,y,z,...,t). + Default: False var_raw : :class:`float` or :any:`None`, optional raw variance of the model which will be multiplied with :any:`CovModel.var_factor` to result in the actual variance. @@ -132,6 +138,7 @@ def __init__( integral_scale=None, rescale=None, latlon=False, + time=False, var_raw=None, hankel_kw=None, **opt_arg, @@ -157,11 +164,13 @@ def __init__( self._nugget_bounds = None self._anis_bounds = None self._opt_arg_bounds = {} - # Set latlon first + # Set latlon and time first self._latlon = bool(latlon) + self._time = bool(time) # SFT class will be created within dim.setter but needs hankel_kw self.hankel_kw = hankel_kw - self.dim = dim + # using time increases model dimension + self.dim = dim + int(self.time) # optional arguments for the variogram-model set_opt_args(self, opt_arg) @@ -177,7 +186,9 @@ def __init__( # set anisotropy and len_scale, disable anisotropy for latlon models self._len_scale, anis = set_len_anis(self.dim, len_scale, anis) if self.latlon: + # keep time anisotropy for metric spatio-temporal model self._anis = np.array((self.dim - 1) * [1], dtype=np.double) + self._anis[-1] = anis[-1] if self.time else 1.0 self._angles = np.array(self.dim * [0], dtype=np.double) else: self._anis = anis @@ -531,14 +542,14 @@ def isometrize(self, pos): """Make a position tuple ready for isotropic operations.""" pos = np.asarray(pos, dtype=np.double).reshape((self.field_dim, -1)) if self.latlon: - return latlon2pos(pos) + return latlon2pos(pos, time=self.time) return np.dot(matrix_isometrize(self.dim, self.angles, self.anis), pos) def anisometrize(self, pos): """Bring a position tuple into the anisotropic coordinate-system.""" pos = np.asarray(pos, dtype=np.double).reshape((self.dim, -1)) if self.latlon: - return pos2latlon(pos) + return pos2latlon(pos, time=self.time) return np.dot( matrix_anisometrize(self.dim, self.angles, self.anis), pos ) @@ -863,6 +874,11 @@ def arg_bounds(self): res.update(self.opt_arg_bounds) return res + @property + def time(self): + """:class:`bool`: Whether the model is a metric spatio-temporal one.""" + return self._time + # geographical coordinates related @property @@ -872,8 +888,13 @@ def latlon(self): @property def field_dim(self): - """:class:`int`: The field dimension of the model.""" - return 2 if self.latlon else self.dim + """:class:`int`: The (parametric) field dimension of the model (with time).""" + return 2 + int(self._time) if self.latlon else self.dim + + @property + def spatial_dim(self): + """:class:`int`: The spatial field dimension of the model (without time).""" + return 2 if self.latlon else self.dim - int(self._time) # standard parameters diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index 2ff5398b..963e6a33 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -522,6 +522,7 @@ def logistic_weights(p=0.1, mean=0.7): # pragma: no cover callable Weighting function. """ + # define the callable weights function def func(x_data): """Callable function for the weights.""" diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 98ed3b8a..37a5dae7 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -498,13 +498,13 @@ def set_dim(model, dim): AttributeWarning, ) dim = model.fix_dim() - if model.latlon and dim != 3: + if model.latlon and dim != (3 + int(model.time)): raise ValueError( f"{model.name}: using fixed dimension {model.fix_dim()}, " - "which is not compatible with a latlon model." + f"which is not compatible with a latlon model (with time={model.time})." ) - # force dim=3 for latlon models - dim = 3 if model.latlon else dim + # force dim=3 (or 4 when time=True) for latlon models + dim = (3 + int(model.time)) if model.latlon else dim # set the dimension if dim < 1: raise ValueError("Only dimensions of d >= 1 are supported.") @@ -551,6 +551,7 @@ def compare(this, that): equal &= np.all(np.isclose(this.angles, that.angles)) equal &= np.isclose(this.rescale, that.rescale) equal &= this.latlon == that.latlon + equal &= this.time == that.time for opt in this.opt_arg: equal &= np.isclose(getattr(this, opt), getattr(that, opt)) return equal diff --git a/src/gstools/tools/geometric.py b/src/gstools/tools/geometric.py index afdcacaf..c05d15ab 100644 --- a/src/gstools/tools/geometric.py +++ b/src/gstools/tools/geometric.py @@ -624,60 +624,75 @@ def ang2dir(angles, dtype=np.double, dim=None): return vec -def latlon2pos(latlon, radius=1.0, dtype=np.double): +def latlon2pos(latlon, radius=1.0, dtype=np.double, time=False): """Convert lat-lon geo coordinates to 3D position tuple. Parameters ---------- latlon : :class:`list` of :class:`numpy.ndarray` latitude and longitude given in degrees. + May includes an appended time axis if `time=True`. radius : :class:`float`, optional Earth radius. Default: `1.0` dtype : data-type, optional The desired data-type for the array. If not given, then the type will be determined as the minimum type required to hold the objects in the sequence. Default: None + time : bool, optional + Whether latlon includes an appended time axis. + Default: False Returns ------- :class:`numpy.ndarray` the 3D position array """ - latlon = np.asarray(latlon, dtype=dtype).reshape((2, -1)) + latlon = np.asarray(latlon, dtype=dtype).reshape((3 if time else 2, -1)) + if time: + timeax = latlon[2] + latlon = latlon[:2] lat, lon = np.deg2rad(latlon) - return np.array( - ( - radius * np.cos(lat) * np.cos(lon), - radius * np.cos(lat) * np.sin(lon), - radius * np.sin(lat) * np.ones_like(lon), - ), - dtype=dtype, + pos_tuple = ( + radius * np.cos(lat) * np.cos(lon), + radius * np.cos(lat) * np.sin(lon), + radius * np.sin(lat) * np.ones_like(lon), ) + if time: + return np.array(pos_tuple + (timeax,), dtype=dtype) + return np.array(pos_tuple, dtype=dtype) -def pos2latlon(pos, radius=1.0, dtype=np.double): +def pos2latlon(pos, radius=1.0, dtype=np.double, time=False): """Convert 3D position tuple from sphere to lat-lon geo coordinates. Parameters ---------- pos : :class:`list` of :class:`numpy.ndarray` The position tuple containing points on a unit-sphere. + May includes an appended time axis if `time=True`. radius : :class:`float`, optional Earth radius. Default: `1.0` dtype : data-type, optional The desired data-type for the array. If not given, then the type will be determined as the minimum type required to hold the objects in the sequence. Default: None + time : bool, optional + Whether latlon includes an appended time axis. + Default: False Returns ------- :class:`numpy.ndarray` the 3D position array """ - pos = np.asarray(pos, dtype=dtype).reshape((3, -1)) + pos = np.asarray(pos, dtype=dtype).reshape((4 if time else 3, -1)) # prevent numerical errors in arcsin lat = np.arcsin(np.maximum(np.minimum(pos[2] / radius, 1.0), -1.0)) lon = np.arctan2(pos[1], pos[0]) + if time: + timeax = pos[3] + lat, lon = np.rad2deg((lat, lon), dtype=dtype) + return np.array((lat, lon, timeax), dtype=dtype) return np.rad2deg((lat, lon), dtype=dtype) From 2ef2f6e7897b3da24a382b54f87e5c7acd3c211d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 6 Jun 2023 15:15:43 +0200 Subject: [PATCH 04/40] time: add time axis to plotter --- src/gstools/field/base.py | 5 +++++ src/gstools/field/plot.py | 44 ++++++++++++++++++++++++++++++--------- 2 files changed, 39 insertions(+), 10 deletions(-) diff --git a/src/gstools/field/base.py b/src/gstools/field/base.py index 903f3893..6098e219 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 time(self): + """:class:`bool`: Whether the field depends on time.""" + return False if self.model is None else self.model.time + @property def name(self): """:class:`str`: The name of the class.""" diff --git a/src/gstools/field/plot.py b/src/gstools/field/plot.py index 528cdcc3..f6242535 100644 --- a/src/gstools/field/plot.py +++ b/src/gstools/field/plot.py @@ -54,7 +54,14 @@ def plot_field( if fld.dim == 1: return plot_1d(fld.pos, fld[field], fig, ax, **kwargs) return plot_nd( - fld.pos, fld[field], fld.mesh_type, fig, ax, fld.latlon, **kwargs + fld.pos, + fld[field], + fld.mesh_type, + fig, + ax, + fld.latlon, + fld.time, + **kwargs, ) @@ -104,6 +111,7 @@ def plot_nd( fig=None, ax=None, latlon=False, + time=False, resolution=128, ax_names=None, aspect="quad", @@ -136,6 +144,11 @@ def plot_nd( ValueError will be raised, if a direction was specified. Bin edges need to be given in radians in this case. Default: False + time : :class:`bool`, optional + Indicate a metric spatio-temporal covariance model. + The time-dimension is assumed to be appended, + meaning the pos tuple is (x,y,z,...,t) or (lat, lon, t). + Default: False resolution : :class:`int`, optional Resolution of the imshow plot. The default is 128. ax_names : :class:`list` of :class:`str`, optional @@ -159,14 +172,20 @@ def plot_nd( """ dim = len(pos) assert dim > 1 - assert not latlon or dim == 2 + assert not latlon or dim == 2 + int(bool(time)) if dim == 2 and contour_plot: return _plot_2d( pos, field, mesh_type, fig, ax, latlon, ax_names, **kwargs ) - pos = pos[::-1] if latlon else pos - field = field.T if (latlon and mesh_type != "unstructured") else field - ax_names = _ax_names(dim, latlon, ax_names) + if latlon: + # swap lat-lon to lon-lat (x-y) + if time: + pos = (pos[1], pos[0], pos[2]) + else: + pos = (pos[1], pos[0]) + if mesh_type != "unstructured": + field = np.moveaxis(field, [0, 1], [1, 0]) + ax_names = _ax_names(dim, latlon, time, ax_names) # init planes planes = rotation_planes(dim) plane_names = [f" {ax_names[p[0]]} - {ax_names[p[1]]}" for p in planes] @@ -323,15 +342,20 @@ def plot_vec_field(fld, field="field", fig=None, ax=None): # pragma: no cover return ax -def _ax_names(dim, latlon=False, ax_names=None): +def _ax_names(dim, latlon=False, time=False, ax_names=None): + t_fac = int(bool(time)) if ax_names is not None: assert len(ax_names) >= dim return ax_names[:dim] - if dim == 2 and latlon: - return ["lon", "lat"] + if dim == 2 + t_fac and latlon: + return ["lon", "lat"] + t_fac * ["time"] if dim <= 3: - return ["$x$", "$y$", "$z$"][:dim] + (dim == 1) * ["field"] - return [f"$x_{{{i}}}$" for i in range(dim)] + return ( + ["$x$", "$y$", "$z$"][:dim] + + t_fac * ["time"] + + (dim == 1) * ["field"] + ) + return [f"$x_{{{i}}}$" for i in range(dim - t_fac)] + t_fac * ["time"] def _plot_2d( From cab5baf89c09de1b73c9f50eb0c3a1510bb61af6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 6 Jun 2023 15:17:47 +0200 Subject: [PATCH 05/40] Examples: update examples with new time attribute --- examples/09_spatio_temporal/01_precip_1d.py | 4 +-- examples/09_spatio_temporal/02_precip_2d.py | 4 +-- .../03_geografic_coordinates.py | 34 +++++++++++++++++++ 3 files changed, 36 insertions(+), 6 deletions(-) create mode 100644 examples/09_spatio_temporal/03_geografic_coordinates.py diff --git a/examples/09_spatio_temporal/01_precip_1d.py b/examples/09_spatio_temporal/01_precip_1d.py index 4671b2f7..2c431c7f 100644 --- a/examples/09_spatio_temporal/01_precip_1d.py +++ b/examples/09_spatio_temporal/01_precip_1d.py @@ -26,13 +26,11 @@ # half daily timesteps over three months t = np.arange(0.0, 90.0, 0.5) -# total spatio-temporal dimension -st_dim = 1 + 1 # space-time anisotropy ratio given in units d / km st_anis = 0.4 # an exponential variogram with a corr. lengths of 2d and 5km -model = gs.Exponential(dim=st_dim, var=1.0, len_scale=5.0, anis=st_anis) +model = gs.Exponential(dim=1, time=True, var=1.0, len_scale=5.0, anis=st_anis) # create a spatial random field instance srf = gs.SRF(model, seed=seed) diff --git a/examples/09_spatio_temporal/02_precip_2d.py b/examples/09_spatio_temporal/02_precip_2d.py index 049225d3..d3d781b3 100644 --- a/examples/09_spatio_temporal/02_precip_2d.py +++ b/examples/09_spatio_temporal/02_precip_2d.py @@ -27,13 +27,11 @@ # half daily timesteps over three months t = np.arange(0.0, 90.0, 0.5) -# total spatio-temporal dimension -st_dim = 2 + 1 # space-time anisotropy ratio given in units d / km st_anis = 0.4 # an exponential variogram with a corr. lengths of 5km, 5km, and 2d -model = gs.Exponential(dim=st_dim, var=1.0, len_scale=5.0, anis=st_anis) +model = gs.Exponential(dim=2, time=True, var=1.0, len_scale=5.0, anis=st_anis) # create a spatial random field instance srf = gs.SRF(model, seed=seed) diff --git a/examples/09_spatio_temporal/03_geografic_coordinates.py b/examples/09_spatio_temporal/03_geografic_coordinates.py new file mode 100644 index 00000000..ea65d128 --- /dev/null +++ b/examples/09_spatio_temporal/03_geografic_coordinates.py @@ -0,0 +1,34 @@ +""" +Working with spatio-temporal lat-lon fields +------------------------------------------- + +In this example, we demonstrate how to generate a spatio-temporal +random field on geographical coordinates. + +First we setup a model, with ``latlon=True`` and ``time=True``, +to get the associated spatio-temporal Yadrenko model. + +In addition, we will use the earth radius provided by :any:`EARTH_RADIUS`, +to have a meaningful length scale in km. + +To generate the field, we simply pass ``(lat, lon, time)`` as the position tuple +to the :any:`SRF` class. + +The anisotropy factor of `0.1` will result in a time length-scale of `77.7` days. +""" +import gstools as gs + +model = gs.Gaussian( + latlon=True, + time=True, + var=1, + len_scale=777, + anis=0.1, + rescale=gs.EARTH_RADIUS, +) + +lat = lon = range(-80, 81) +time = range(0, 110, 10) +srf = gs.SRF(model, seed=1234) +field = srf.structured((lat, lon, time)) +srf.plot() From e6867d65351f31f0d7489373ed2c51effcad1262 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 6 Jun 2023 15:26:56 +0200 Subject: [PATCH 06/40] pylint: ignore 'use-dict-literal', increase max limits --- pyproject.toml | 4 ++-- src/gstools/transform/field.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) 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/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 ( From a42daedd669e529386f888676bbff7d47b4984da Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 8 Jun 2023 16:46:40 +0200 Subject: [PATCH 07/40] CovModel: add radius property; correctly scale time axis for latlon; add DEGREE_SCALE --- src/gstools/__init__.py | 5 ++++- src/gstools/covmodel/base.py | 33 ++++++++++++++++++++++++++++----- src/gstools/covmodel/fit.py | 2 +- src/gstools/covmodel/plot.py | 6 +++--- src/gstools/tools/__init__.py | 6 ++++++ src/gstools/tools/geometric.py | 32 +++++++++++++++++++------------- 6 files changed, 61 insertions(+), 23 deletions(-) diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py index 60dd1e33..8e72fdfa 100644 --- a/src/gstools/__init__.py +++ b/src/gstools/__init__.py @@ -125,6 +125,7 @@ .. autosummary:: EARTH_RADIUS + DEGREE_SCALE """ # Hooray! from gstools import ( @@ -161,6 +162,7 @@ from gstools.field import SRF, CondSRF from gstools.krige import Krige from gstools.tools import ( + DEGREE_SCALE, EARTH_RADIUS, generate_grid, generate_st_grid, @@ -188,7 +190,7 @@ __all__ = ["__version__"] __all__ += ["covmodel", "field", "variogram", "krige", "random", "tools"] -__all__ += ["transform", "normalizer"] +__all__ += ["transform", "normalizer", "config"] __all__ += [ "CovModel", "Gaussian", @@ -226,6 +228,7 @@ "generate_grid", "generate_st_grid", "EARTH_RADIUS", + "DEGREE_SCALE", "vtk_export", "vtk_export_structured", "vtk_export_unstructured", diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 3a923c1c..a2740acf 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -103,6 +103,12 @@ class CovModel: disabled. `rescale` can be set to e.g. earth's radius, to have a meaningful `len_scale` parameter. Default: False + radius : :class:`float`, optional + Sphere radius in case of latlon coordinates to get a meaningful length + scale. By default, len_scale is assumed to be in radians with latlon=True. + Can be set to :any:`EARTH_RADIUS` to have len_scale in km or + :any:`DEGREE_SCALE` to have len_scale in degree. + Default: ``1.0`` time : :class:`bool`, optional Create a metric spatio-temporal covariance model. Setting this to true will increase `dim` and `field_dim` by 1. @@ -138,6 +144,7 @@ def __init__( integral_scale=None, rescale=None, latlon=False, + radius=1.0, time=False, var_raw=None, hankel_kw=None, @@ -167,6 +174,7 @@ def __init__( # Set latlon and time first self._latlon = bool(latlon) self._time = bool(time) + self._radius = abs(float(radius)) # SFT class will be created within dim.setter but needs hankel_kw self.hankel_kw = hankel_kw # using time increases model dimension @@ -254,15 +262,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(2 * np.sin(zeta / 2) * self.radius) 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(2 * np.sin(zeta / 2) * self.radius) 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(2 * np.sin(zeta / 2) * self.radius) def vario_spatial(self, pos): r"""Spatial variogram respecting anisotropy and rotation.""" @@ -542,14 +550,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, time=self.time) + return latlon2pos( + pos, + radius=self.radius, + time=self.time, + 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, time=self.time) + return pos2latlon( + pos, + radius=self.radius, + time=self.time, + time_scale=self.anis[-1], + ) return np.dot( matrix_anisometrize(self.dim, self.angles, self.anis), pos ) @@ -886,6 +904,11 @@ def latlon(self): """:class:`bool`: Whether the model depends on geographical coords.""" return self._latlon + @property + def radius(self): + """:class:`float`: Sphere radius for geographical coords.""" + return self._radius + @property def field_dim(self): """:class:`int`: The (parametric) field dimension of the model (with time).""" diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index 963e6a33..1ab04e9f 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -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 = 2 * np.sin(x_data / 2) * model.radius return x_data, y_data, is_dir_vario diff --git a/src/gstools/covmodel/plot.py b/src/gstools/covmodel/plot.py index efcc2630..cab3091f 100644 --- a/src/gstools/covmodel/plot.py +++ b/src/gstools/covmodel/plot.py @@ -150,7 +150,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.radius, 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 +165,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.radius, 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 +180,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.radius, 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/tools/__init__.py b/src/gstools/tools/__init__.py index bd329576..79319c49 100644 --- a/src/gstools/tools/__init__.py +++ b/src/gstools/tools/__init__.py @@ -58,10 +58,13 @@ .. autosummary:: EARTH_RADIUS + DEGREE_SCALE ---- .. autodata:: EARTH_RADIUS + +.. autodata:: DEGREE_SCALE """ from gstools.tools.export import ( @@ -103,6 +106,9 @@ EARTH_RADIUS = 6371.0 """float: earth radius for WGS84 ellipsoid in km""" +DEGREE_SCALE = 57.29577951308232 +"""float: radius for unit sphere in degree""" + __all__ = [ "vtk_export", diff --git a/src/gstools/tools/geometric.py b/src/gstools/tools/geometric.py index c05d15ab..4734dcbb 100644 --- a/src/gstools/tools/geometric.py +++ b/src/gstools/tools/geometric.py @@ -624,7 +624,9 @@ def ang2dir(angles, dtype=np.double, dim=None): return vec -def latlon2pos(latlon, radius=1.0, dtype=np.double, time=False): +def latlon2pos( + latlon, radius=1.0, dtype=np.double, time=False, time_scale=1.0 +): """Convert lat-lon geo coordinates to 3D position tuple. Parameters @@ -638,9 +640,12 @@ def latlon2pos(latlon, radius=1.0, dtype=np.double, time=False): The desired data-type for the array. If not given, then the type will be determined as the minimum type required to hold the objects in the sequence. Default: None - time : bool, optional + time : :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 ------- @@ -648,21 +653,18 @@ def latlon2pos(latlon, radius=1.0, dtype=np.double, time=False): the 3D position array """ latlon = np.asarray(latlon, dtype=dtype).reshape((3 if time else 2, -1)) - if time: - timeax = latlon[2] - latlon = latlon[:2] - lat, lon = np.deg2rad(latlon) + 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 time: - return np.array(pos_tuple + (timeax,), dtype=dtype) + 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, time=False): +def pos2latlon(pos, radius=1.0, dtype=np.double, time=False, time_scale=1.0): """Convert 3D position tuple from sphere to lat-lon geo coordinates. Parameters @@ -676,9 +678,12 @@ def pos2latlon(pos, radius=1.0, dtype=np.double, time=False): The desired data-type for the array. If not given, then the type will be determined as the minimum type required to hold the objects in the sequence. Default: None - time : bool, optional + time : :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 ------- @@ -689,11 +694,12 @@ def pos2latlon(pos, radius=1.0, dtype=np.double, time=False): # 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]) + latlon = np.rad2deg((lat, lon), dtype=dtype) if time: - timeax = pos[3] - lat, lon = np.rad2deg((lat, lon), dtype=dtype) - return np.array((lat, lon, timeax), dtype=dtype) - return np.rad2deg((lat, lon), dtype=dtype) + return np.array( + (latlon[0], latlon[1], pos[3] * time_scale), dtype=dtype + ) + return latlon def chordal_to_great_circle(dist): From a4d3d3dc52c0ef3d43d846988c95b017cd411f90 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 8 Jun 2023 16:51:23 +0200 Subject: [PATCH 08/40] examples: use radius instead of rescale for latlon models now --- examples/03_variogram/06_auto_bin_latlon.py | 2 +- examples/08_geo_coordinates/01_dwd_krige.py | 4 ++-- examples/08_geo_coordinates/README.rst | 4 ++-- ...fic_coordinates.py => 03_geographic_coordinates.py} | 10 ++++++---- 4 files changed, 11 insertions(+), 9 deletions(-) rename examples/09_spatio_temporal/{03_geografic_coordinates.py => 03_geographic_coordinates.py} (78%) diff --git a/examples/03_variogram/06_auto_bin_latlon.py b/examples/03_variogram/06_auto_bin_latlon.py index 70c0a09b..f0e9d8d2 100644 --- a/examples/03_variogram/06_auto_bin_latlon.py +++ b/examples/03_variogram/06_auto_bin_latlon.py @@ -45,7 +45,7 @@ # data-variance as additional information during the fit of the variogram. emp_v = gs.vario_estimate(pos, field, latlon=True) -sph = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS) +sph = gs.Spherical(latlon=True, radius=gs.EARTH_RADIUS) sph.fit_variogram(*emp_v, sill=np.var(field)) ax = sph.plot(x_max=2 * np.max(emp_v[0])) ax.scatter(*emp_v, label="Empirical variogram") diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py index b37e7fa0..4c665b9d 100755 --- a/examples/08_geo_coordinates/01_dwd_krige.py +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -86,7 +86,7 @@ def get_dwd_temperature(date="2020-06-09 12:00:00"): # Now we can use this estimated variogram to fit a model to it. # Here we will use a :any:`Spherical` model. We select the ``latlon`` option # to use the `Yadrenko` variant of the model to gain a valid model for lat-lon -# coordinates and we rescale it to the earth-radius. Otherwise the length +# coordinates and we set the radius to the earth-radius. Otherwise the length # scale would be given in radians representing the great-circle distance. # # We deselect the nugget from fitting and plot the result afterwards. @@ -97,7 +97,7 @@ def get_dwd_temperature(date="2020-06-09 12:00:00"): # still holds the ordinary routine that is not respecting the great-circle # distance. -model = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS) +model = gs.Spherical(latlon=True, radius=gs.EARTH_RADIUS) model.fit_variogram(bin_c, vario, nugget=False) ax = model.plot("vario_yadrenko", x_max=bins[-1]) ax.scatter(bin_c, vario) diff --git a/examples/08_geo_coordinates/README.rst b/examples/08_geo_coordinates/README.rst index 87b419df..895cd250 100644 --- a/examples/08_geo_coordinates/README.rst +++ b/examples/08_geo_coordinates/README.rst @@ -22,14 +22,14 @@ in your desired model (see :any:`CovModel`): By doing so, the model will use the associated `Yadrenko` model on a sphere (see `here `_). 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 ``radius`` 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, radius=gs.EARTH_RADIUS) Then ``len_scale`` can be interpreted as given in km. diff --git a/examples/09_spatio_temporal/03_geografic_coordinates.py b/examples/09_spatio_temporal/03_geographic_coordinates.py similarity index 78% rename from examples/09_spatio_temporal/03_geografic_coordinates.py rename to examples/09_spatio_temporal/03_geographic_coordinates.py index ea65d128..f1d2959d 100644 --- a/examples/09_spatio_temporal/03_geografic_coordinates.py +++ b/examples/09_spatio_temporal/03_geographic_coordinates.py @@ -14,8 +14,10 @@ To generate the field, we simply pass ``(lat, lon, time)`` as the position tuple to the :any:`SRF` class. -The anisotropy factor of `0.1` will result in a time length-scale of `77.7` days. +The anisotropy factor of `0.1` (days/km) will result in a time length-scale of `77.7` days. """ +import numpy as np + import gstools as gs model = gs.Gaussian( @@ -24,11 +26,11 @@ var=1, len_scale=777, anis=0.1, - rescale=gs.EARTH_RADIUS, + radius=gs.EARTH_RADIUS, ) -lat = lon = range(-80, 81) -time = range(0, 110, 10) +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() From 456f34d20f65d0ca838d6d134ee76327f581ceaf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 8 Jun 2023 17:56:26 +0200 Subject: [PATCH 09/40] CovModel: update __repr__ --- src/gstools/covmodel/tools.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 37a5dae7..5ac3bbd2 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -569,21 +569,29 @@ def model_repr(model): # pragma: no cover m = model p = model._prec opt_str = "" + t_str = ", time=True" if m.time 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.time else f", anis={m.anis[-1]:.{p}}" + ) + r_str = "" if np.isclose(m.radius, 1) else f", radius={m.radius:.{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})" ) From bd536f87008845b2b388c739a61ef263dc3f9e7e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 09:19:27 +0200 Subject: [PATCH 10/40] CovModel: rename 'radius' to 'geo_scale' and add more scales --- examples/03_variogram/06_auto_bin_latlon.py | 2 +- .../08_geo_coordinates/00_field_generation.py | 8 +++--- examples/08_geo_coordinates/01_dwd_krige.py | 4 +-- examples/08_geo_coordinates/README.rst | 4 +-- .../03_geographic_coordinates.py | 12 ++++----- src/gstools/__init__.py | 4 +++ src/gstools/covmodel/base.py | 26 +++++++++---------- src/gstools/covmodel/fit.py | 2 +- src/gstools/covmodel/plot.py | 6 ++--- src/gstools/covmodel/tools.py | 6 ++++- src/gstools/tools/__init__.py | 9 +++++++ src/gstools/tools/geometric.py | 4 +-- 12 files changed, 53 insertions(+), 34 deletions(-) diff --git a/examples/03_variogram/06_auto_bin_latlon.py b/examples/03_variogram/06_auto_bin_latlon.py index f0e9d8d2..83216b2a 100644 --- a/examples/03_variogram/06_auto_bin_latlon.py +++ b/examples/03_variogram/06_auto_bin_latlon.py @@ -45,7 +45,7 @@ # data-variance as additional information during the fit of the variogram. emp_v = gs.vario_estimate(pos, field, latlon=True) -sph = gs.Spherical(latlon=True, radius=gs.EARTH_RADIUS) +sph = gs.Spherical(latlon=True, geo_scale=gs.EARTH_RADIUS) sph.fit_variogram(*emp_v, sill=np.var(field)) ax = sph.plot(x_max=2 * np.max(emp_v[0])) ax.scatter(*emp_v, label="Empirical variogram") diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py index b5685c7d..e2530f36 100755 --- a/examples/08_geo_coordinates/00_field_generation.py +++ b/examples/08_geo_coordinates/00_field_generation.py @@ -8,15 +8,17 @@ First we setup a model, with ``latlon=True``, to get the associated Yadrenko model. -In addition, we will use the earth radius provided by :any:`EARTH_RADIUS`, -to have a meaningful length scale in km. +In addition, we will use the earth radius provided by :any:`EARTH_RADIUS` +as ``geo_scale`` to have a meaningful length scale in km. To generate the field, we simply pass ``(lat, lon)`` as the position tuple to the :any:`SRF` class. """ import gstools as gs -model = gs.Gaussian(latlon=True, var=1, len_scale=777, rescale=gs.EARTH_RADIUS) +model = gs.Gaussian( + latlon=True, var=1, len_scale=777, geo_scale=gs.EARTH_RADIUS +) lat = lon = range(-80, 81) srf = gs.SRF(model, seed=1234) diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py index 4c665b9d..39de1b55 100755 --- a/examples/08_geo_coordinates/01_dwd_krige.py +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -86,7 +86,7 @@ def get_dwd_temperature(date="2020-06-09 12:00:00"): # Now we can use this estimated variogram to fit a model to it. # Here we will use a :any:`Spherical` model. We select the ``latlon`` option # to use the `Yadrenko` variant of the model to gain a valid model for lat-lon -# coordinates and we set the radius to the earth-radius. Otherwise the length +# coordinates and we set the ``geo_scale`` to the earth-radius. Otherwise the length # scale would be given in radians representing the great-circle distance. # # We deselect the nugget from fitting and plot the result afterwards. @@ -97,7 +97,7 @@ def get_dwd_temperature(date="2020-06-09 12:00:00"): # still holds the ordinary routine that is not respecting the great-circle # distance. -model = gs.Spherical(latlon=True, radius=gs.EARTH_RADIUS) +model = gs.Spherical(latlon=True, geo_scale=gs.EARTH_RADIUS) model.fit_variogram(bin_c, vario, nugget=False) ax = model.plot("vario_yadrenko", x_max=bins[-1]) ax.scatter(bin_c, vario) diff --git a/examples/08_geo_coordinates/README.rst b/examples/08_geo_coordinates/README.rst index 895cd250..76083ce8 100644 --- a/examples/08_geo_coordinates/README.rst +++ b/examples/08_geo_coordinates/README.rst @@ -22,14 +22,14 @@ in your desired model (see :any:`CovModel`): By doing so, the model will use the associated `Yadrenko` model on a sphere (see `here `_). 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 ``radius`` +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, radius=gs.EARTH_RADIUS) + model = gs.Gaussian(latlon=True, var=2, len_scale=500, geo_scale=gs.EARTH_RADIUS) Then ``len_scale`` can be interpreted as given in km. diff --git a/examples/09_spatio_temporal/03_geographic_coordinates.py b/examples/09_spatio_temporal/03_geographic_coordinates.py index f1d2959d..39cc31ff 100644 --- a/examples/09_spatio_temporal/03_geographic_coordinates.py +++ b/examples/09_spatio_temporal/03_geographic_coordinates.py @@ -8,25 +8,25 @@ First we setup a model, with ``latlon=True`` and ``time=True``, to get the associated spatio-temporal Yadrenko model. -In addition, we will use the earth radius provided by :any:`EARTH_RADIUS`, -to have a meaningful length scale in km. +In addition, we will use the earth radius provided by :any:`EARTH_RADIUS` + as ``geo_scale`` to have a meaningful length scale in km. To generate the field, we simply pass ``(lat, lon, time)`` as the position tuple to the :any:`SRF` class. -The anisotropy factor of `0.1` (days/km) will result in a time length-scale of `77.7` days. +The anisotropy factor of `0.1` (days/km) will result in a time length-scale of `100` days. """ import numpy as np import gstools as gs -model = gs.Gaussian( +model = gs.Matern( latlon=True, time=True, var=1, - len_scale=777, + len_scale=1000, anis=0.1, - radius=gs.EARTH_RADIUS, + geo_scale=gs.EARTH_RADIUS, ) lat = lon = np.linspace(-80, 81, 50) diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py index 8e72fdfa..2404a4ff 100644 --- a/src/gstools/__init__.py +++ b/src/gstools/__init__.py @@ -164,6 +164,8 @@ from gstools.tools import ( DEGREE_SCALE, EARTH_RADIUS, + KM_SCALE, + RADIAN_SCALE, generate_grid, generate_st_grid, rotated_main_axes, @@ -228,7 +230,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 a2740acf..4bbfd78f 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -100,11 +100,11 @@ 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 - radius : :class:`float`, optional - Sphere radius in case of latlon coordinates to get a meaningful length + geo_scale : :class:`float`, optional + Geographic scaling in case of latlon coordinates to get a meaningful length scale. By default, len_scale is assumed to be in radians with latlon=True. Can be set to :any:`EARTH_RADIUS` to have len_scale in km or :any:`DEGREE_SCALE` to have len_scale in degree. @@ -144,7 +144,7 @@ def __init__( integral_scale=None, rescale=None, latlon=False, - radius=1.0, + geo_scale=1.0, time=False, var_raw=None, hankel_kw=None, @@ -174,7 +174,7 @@ def __init__( # Set latlon and time first self._latlon = bool(latlon) self._time = bool(time) - self._radius = abs(float(radius)) + 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 @@ -262,15 +262,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) * self.radius) + return self.variogram(2 * np.sin(zeta / 2) * 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) * self.radius) + return self.covariance(2 * np.sin(zeta / 2) * 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) * self.radius) + return self.correlation(2 * np.sin(zeta / 2) * self.geo_scale) def vario_spatial(self, pos): r"""Spatial variogram respecting anisotropy and rotation.""" @@ -552,7 +552,7 @@ def isometrize(self, pos): if self.latlon: return latlon2pos( pos, - radius=self.radius, + radius=self.geo_scale, time=self.time, time_scale=self.anis[-1], ) @@ -564,7 +564,7 @@ def anisometrize(self, pos): if self.latlon: return pos2latlon( pos, - radius=self.radius, + radius=self.geo_scale, time=self.time, time_scale=self.anis[-1], ) @@ -905,9 +905,9 @@ def latlon(self): return self._latlon @property - def radius(self): - """:class:`float`: Sphere radius for geographical coords.""" - return self._radius + def geo_scale(self): + """:class:`float`: Geographic scaling for geographical coords.""" + return self._geo_scale @property def field_dim(self): diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index 1ab04e9f..00b8b2e5 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -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) * model.radius + x_data = 2 * np.sin(x_data / 2) * model.geo_scale return x_data, y_data, is_dir_vario diff --git a/src/gstools/covmodel/plot.py b/src/gstools/covmodel/plot.py index cab3091f..414a95d8 100644 --- a/src/gstools/covmodel/plot.py +++ b/src/gstools/covmodel/plot.py @@ -150,7 +150,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_scale / model.radius, 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 +165,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_scale / model.radius, 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 +180,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_scale / model.radius, 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 5ac3bbd2..a1bef143 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -578,7 +578,11 @@ def model_repr(model): # pragma: no cover ani_str = ( "" if m.is_isotropic or not m.time else f", anis={m.anis[-1]:.{p}}" ) - r_str = "" if np.isclose(m.radius, 1) else f", radius={m.radius:.{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}{t_str}, var={m.var:.{p}}, " f"len_scale={m.len_scale:.{p}}, nugget={m.nugget:.{p}}" diff --git a/src/gstools/tools/__init__.py b/src/gstools/tools/__init__.py index 79319c49..a5977839 100644 --- a/src/gstools/tools/__init__.py +++ b/src/gstools/tools/__init__.py @@ -106,9 +106,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", @@ -141,4 +147,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 4734dcbb..dcb85cb3 100644 --- a/src/gstools/tools/geometric.py +++ b/src/gstools/tools/geometric.py @@ -635,7 +635,7 @@ def latlon2pos( 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 @@ -673,7 +673,7 @@ def pos2latlon(pos, radius=1.0, dtype=np.double, time=False, time_scale=1.0): 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 From cc6d75b44f940bace44b770546fc9c2ba525c7d1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 09:30:10 +0200 Subject: [PATCH 11/40] Better geo_scale documentation --- examples/08_geo_coordinates/00_field_generation.py | 4 +--- src/gstools/__init__.py | 2 ++ src/gstools/covmodel/base.py | 5 +++-- src/gstools/tools/__init__.py | 6 ++++++ 4 files changed, 12 insertions(+), 5 deletions(-) diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py index e2530f36..9feef660 100755 --- a/examples/08_geo_coordinates/00_field_generation.py +++ b/examples/08_geo_coordinates/00_field_generation.py @@ -16,9 +16,7 @@ """ import gstools as gs -model = gs.Gaussian( - latlon=True, var=1, len_scale=777, geo_scale=gs.EARTH_RADIUS -) +model = gs.Gaussian(latlon=True, len_scale=777, geo_scale=gs.EARTH_RADIUS) lat = lon = range(-80, 81) srf = gs.SRF(model, seed=1234) diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py index 2404a4ff..6d62d558 100644 --- a/src/gstools/__init__.py +++ b/src/gstools/__init__.py @@ -125,7 +125,9 @@ .. autosummary:: EARTH_RADIUS + KM_SCALE DEGREE_SCALE + RADIAN_SCALE """ # Hooray! from gstools import ( diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 4bbfd78f..2e965841 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -31,6 +31,7 @@ set_opt_args, spectral_rad_pdf, ) +from gstools.tools import RADIAN_SCALE from gstools.tools.geometric import ( latlon2pos, matrix_anisometrize, @@ -108,7 +109,7 @@ class CovModel: scale. By default, len_scale is assumed to be in radians with latlon=True. Can be set to :any:`EARTH_RADIUS` to have len_scale in km or :any:`DEGREE_SCALE` to have len_scale in degree. - Default: ``1.0`` + Default: :any:`RADIAN_SCALE` time : :class:`bool`, optional Create a metric spatio-temporal covariance model. Setting this to true will increase `dim` and `field_dim` by 1. @@ -144,7 +145,7 @@ def __init__( integral_scale=None, rescale=None, latlon=False, - geo_scale=1.0, + geo_scale=RADIAN_SCALE, time=False, var_raw=None, hankel_kw=None, diff --git a/src/gstools/tools/__init__.py b/src/gstools/tools/__init__.py index a5977839..1f68dbaf 100644 --- a/src/gstools/tools/__init__.py +++ b/src/gstools/tools/__init__.py @@ -58,13 +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 ( From 4746c8d93dc91d23ed703fa7f7810ec4fcb30de4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 10:50:51 +0200 Subject: [PATCH 12/40] examples: fix typo --- examples/09_spatio_temporal/03_geographic_coordinates.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/09_spatio_temporal/03_geographic_coordinates.py b/examples/09_spatio_temporal/03_geographic_coordinates.py index 39cc31ff..3ec8d28c 100644 --- a/examples/09_spatio_temporal/03_geographic_coordinates.py +++ b/examples/09_spatio_temporal/03_geographic_coordinates.py @@ -9,7 +9,7 @@ to get the associated spatio-temporal Yadrenko model. In addition, we will use the earth radius provided by :any:`EARTH_RADIUS` - as ``geo_scale`` to have a meaningful length scale in km. +as ``geo_scale`` to have a meaningful length scale in km. To generate the field, we simply pass ``(lat, lon, time)`` as the position tuple to the :any:`SRF` class. From af6f5220f51c90ed22ab6fec342f8b94f07b230f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 10:55:35 +0200 Subject: [PATCH 13/40] vario: rename 'bin_centres' to 'bin_center' following doc-string --- src/gstools/variogram/variogram.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/gstools/variogram/variogram.py b/src/gstools/variogram/variogram.py index b57e0354..f259bf00 100644 --- a/src/gstools/variogram/variogram.py +++ b/src/gstools/variogram/variogram.py @@ -250,18 +250,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_center = (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_center = np.empty(0) if bin_edges is None else bin_center + estimates = np.zeros_like(bin_center) if return_counts: - return bin_centres, estimates, np.zeros_like(estimates, dtype=int) - return bin_centres, estimates + return bin_center, estimates, np.zeros_like(estimates, dtype=int) + return bin_center, estimates if not masked: field = field.filled() # check mesh shape @@ -334,7 +334,7 @@ def vario_estimate( # create bining if not given if bin_edges is None: bin_edges = standard_bins(pos, dim, latlon) - bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + bin_center = (bin_edges[:-1] + bin_edges[1:]) / 2.0 # normalize field norm_field_out = remove_trend_norm_mean( *(pos, field, mean, normalizer, trend), @@ -371,7 +371,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_center,) + est_out[: 2 if return_counts else 1] + norm_out def vario_estimate_axis( From 81e96167c84ebf348becd2ecfe26a61484b7b600 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 12:10:44 +0200 Subject: [PATCH 14/40] tools: add great_circle_to_chordal; add radius to chordal_to_great_circle --- src/gstools/tools/geometric.py | 34 ++++++++++++++++++++++++++++++---- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/src/gstools/tools/geometric.py b/src/gstools/tools/geometric.py index dcb85cb3..fcbd7c68 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 @@ -702,14 +703,16 @@ def pos2latlon(pos, radius=1.0, dtype=np.double, time=False, time_scale=1.0): 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 ------- @@ -718,6 +721,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)) From 97c30b485ee688aa6a80c4f51cd165146df4ed8a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 12:11:59 +0200 Subject: [PATCH 15/40] vario: add geo_scale to variogram estimation routines --- src/gstools/covmodel/base.py | 7 ++++--- src/gstools/covmodel/fit.py | 4 ++-- src/gstools/covmodel/plot.py | 6 +++--- src/gstools/krige/base.py | 5 ++++- src/gstools/variogram/binning.py | 12 ++++++++++-- src/gstools/variogram/variogram.py | 13 ++++++++++++- 6 files changed, 35 insertions(+), 12 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 2e965841..8502be99 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -33,6 +33,7 @@ ) from gstools.tools import RADIAN_SCALE from gstools.tools.geometric import ( + great_circle_to_chordal, latlon2pos, matrix_anisometrize, matrix_isometrize, @@ -263,15 +264,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) * self.geo_scale) + 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) * self.geo_scale) + 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) * self.geo_scale) + return self.correlation(great_circle_to_chordal(zeta, self.geo_scale)) def vario_spatial(self, pos): r"""Spatial variogram respecting anisotropy and rotation.""" diff --git a/src/gstools/covmodel/fit.py b/src/gstools/covmodel/fit.py index 00b8b2e5..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) * model.geo_scale + x_data = great_circle_to_chordal(x_data, model.geo_scale) return x_data, y_data, is_dir_vario diff --git a/src/gstools/covmodel/plot.py b/src/gstools/covmodel/plot.py index 414a95d8..fa526612 100644 --- a/src/gstools/covmodel/plot.py +++ b/src/gstools/covmodel/plot.py @@ -150,7 +150,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_scale / model.geo_scale, np.pi) + x_max = 3 * model.len_scale 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 +165,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_scale / model.geo_scale, np.pi) + x_max = 3 * model.len_scale 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 +180,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_scale / model.geo_scale, np.pi) + x_max = 3 * model.len_scale 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/krige/base.py b/src/gstools/krige/base.py index 77c6832e..4d1ee824 100755 --- a/src/gstools/krige/base.py +++ b/src/gstools/krige/base.py @@ -527,7 +527,10 @@ def set_condition( 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/variogram/binning.py b/src/gstools/variogram/binning.py index be490110..891b39e5 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,12 @@ 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 scaling in case of latlon coordinates to get meaningful bins. + By default, bins are assumed to be given in radians with latlon=True. + Can be set to :any:`EARTH_RADIUS` to have units in km or + :any:`DEGREE_SCALE` to have units in degree. + Default: :any:`RADIAN_SCALE` Returns ------- @@ -80,7 +88,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 +96,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 f259bf00..e179dc28 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,7 @@ def vario_estimate( normalizer=None, trend=None, fit_normalizer=False, + geo_scale=RADIAN_SCALE, ): r""" Estimates the empirical variogram. @@ -222,6 +224,12 @@ 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 scaling in case of latlon coordinates to get meaningful bins. + By default, bins are assumed to be given in radians with latlon=True. + Can be set to :any:`EARTH_RADIUS` to have units in km or + :any:`DEGREE_SCALE` to have units in degree. + Default: :any:`RADIAN_SCALE` Returns ------- @@ -333,8 +341,11 @@ def vario_estimate( pos = pos[:, sampled_idx] # create bining if not given if bin_edges is None: - bin_edges = standard_bins(pos, dim, latlon) + bin_edges = standard_bins(pos, dim, latlon, geo_scale=geo_scale) bin_center = (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), From b96f844862e3f2f33139e67c2e554f119e8bb6c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 12:20:28 +0200 Subject: [PATCH 16/40] vario: forward kwargs to standard_bins routine --- src/gstools/variogram/variogram.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/gstools/variogram/variogram.py b/src/gstools/variogram/variogram.py index e179dc28..8c5785c6 100644 --- a/src/gstools/variogram/variogram.py +++ b/src/gstools/variogram/variogram.py @@ -94,6 +94,7 @@ def vario_estimate( trend=None, fit_normalizer=False, geo_scale=RADIAN_SCALE, + **std_bins, ): r""" Estimates the empirical variogram. @@ -230,6 +231,9 @@ def vario_estimate( Can be set to :any:`EARTH_RADIUS` to have units in km or :any:`DEGREE_SCALE` to have units in degree. 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 ------- @@ -341,7 +345,9 @@ def vario_estimate( pos = pos[:, sampled_idx] # create bining if not given if bin_edges is None: - bin_edges = standard_bins(pos, dim, latlon, geo_scale=geo_scale) + bin_edges = standard_bins( + pos, dim, latlon, geo_scale=geo_scale, **std_bins + ) bin_center = (bin_edges[:-1] + bin_edges[1:]) / 2.0 if latlon: # internally we always use radians From 5486a1c6be8f277c48d38f2c97fee2d03546b13d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 12:20:49 +0200 Subject: [PATCH 17/40] examples: update examples for geo_scale --- .../08_geo_coordinates/00_field_generation.py | 9 ++++++--- examples/08_geo_coordinates/01_dwd_krige.py | 15 ++++++++------- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py index 9feef660..83700d61 100755 --- a/examples/08_geo_coordinates/00_field_generation.py +++ b/examples/08_geo_coordinates/00_field_generation.py @@ -14,6 +14,8 @@ To generate the field, we simply pass ``(lat, lon)`` as the position tuple to the :any:`SRF` class. """ +import numpy as np + import gstools as gs model = gs.Gaussian(latlon=True, len_scale=777, geo_scale=gs.EARTH_RADIUS) @@ -32,7 +34,7 @@ # # As we will see, everthing went well... phew! -bin_edges = [0.01 * i for i in range(30)] +bin_edges = np.linspace(0, 777 * 3, 30) bin_center, emp_vario = gs.vario_estimate( (lat, lon), field, @@ -41,11 +43,12 @@ mesh_type="structured", sampling_size=2000, sampling_seed=12345, + geo_scale=gs.EARTH_RADIUS, ) -ax = model.plot("vario_yadrenko", x_max=0.3) +ax = model.plot("vario_yadrenko", x_max=max(bin_center)) model.fit_variogram(bin_center, emp_vario, nugget=False) -model.plot("vario_yadrenko", ax=ax, label="fitted", x_max=0.3) +model.plot("vario_yadrenko", ax=ax, label="fitted", x_max=max(bin_center)) ax.scatter(bin_center, emp_vario, color="k") print(model) diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py index 39de1b55..f032d0e3 100755 --- a/examples/08_geo_coordinates/01_dwd_krige.py +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -76,11 +76,11 @@ def get_dwd_temperature(date="2020-06-09 12:00:00"): ############################################################################### # First we will estimate the variogram of our temperature data. -# As the maximal bin distance we choose 8 degrees, which corresponds to a -# chordal length of about 900 km. +# As the maximal bin distance we choose 900 km. -bins = gs.standard_bins((lat, lon), max_dist=np.deg2rad(8), latlon=True) -bin_c, vario = gs.vario_estimate((lat, lon), temp, bins, latlon=True) +bin_center, vario = gs.vario_estimate( + (lat, lon), temp, latlon=True, geo_scale=gs.EARTH_RADIUS, max_dist=900 +) ############################################################################### # Now we can use this estimated variogram to fit a model to it. @@ -98,9 +98,9 @@ def get_dwd_temperature(date="2020-06-09 12:00:00"): # distance. model = gs.Spherical(latlon=True, geo_scale=gs.EARTH_RADIUS) -model.fit_variogram(bin_c, vario, nugget=False) -ax = model.plot("vario_yadrenko", x_max=bins[-1]) -ax.scatter(bin_c, vario) +model.fit_variogram(bin_center, vario, nugget=False) +ax = model.plot("vario_yadrenko", x_max=max(bin_center)) +ax.scatter(bin_center, vario) print(model) ############################################################################### @@ -171,3 +171,4 @@ def north_south_drift(lat, lon): ############################################################################### # Interpretion of the results is now up to you! ;-) +plt.show() From 1e214655b05e564f82c4ce73b4979000f0d397a3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 12:26:20 +0200 Subject: [PATCH 18/40] debug fix --- examples/08_geo_coordinates/01_dwd_krige.py | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py index f032d0e3..98c570b4 100755 --- a/examples/08_geo_coordinates/01_dwd_krige.py +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -171,4 +171,3 @@ def north_south_drift(lat, lon): ############################################################################### # Interpretion of the results is now up to you! ;-) -plt.show() From 14ccf3d48d1d7369920adc13b05d14414524f3a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 12:26:39 +0200 Subject: [PATCH 19/40] update latlon auto-bin example with geo_scale --- examples/03_variogram/06_auto_bin_latlon.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/03_variogram/06_auto_bin_latlon.py b/examples/03_variogram/06_auto_bin_latlon.py index 83216b2a..7f3b97fb 100644 --- a/examples/03_variogram/06_auto_bin_latlon.py +++ b/examples/03_variogram/06_auto_bin_latlon.py @@ -44,10 +44,10 @@ # Since the overall range of these meteo-stations is too low, we can use the # data-variance as additional information during the fit of the variogram. -emp_v = gs.vario_estimate(pos, field, latlon=True) +emp_v = gs.vario_estimate(pos, field, latlon=True, geo_scale=gs.EARTH_RADIUS) sph = gs.Spherical(latlon=True, geo_scale=gs.EARTH_RADIUS) sph.fit_variogram(*emp_v, sill=np.var(field)) -ax = sph.plot(x_max=2 * np.max(emp_v[0])) +ax = sph.plot("vario_yadrenko", x_max=2 * np.max(emp_v[0])) ax.scatter(*emp_v, label="Empirical variogram") ax.legend() print(sph) From 8223e2bea418362e3f2873a5eca38f02a13399e2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 12:49:34 +0200 Subject: [PATCH 20/40] examples: update readme for geo_scale --- examples/08_geo_coordinates/README.rst | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/examples/08_geo_coordinates/README.rst b/examples/08_geo_coordinates/README.rst index 76083ce8..b664512f 100644 --- a/examples/08_geo_coordinates/README.rst +++ b/examples/08_geo_coordinates/README.rst @@ -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:: From d8a1309b920384ef77ef6fa48e6b33de2b03f3cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 16:35:33 +0200 Subject: [PATCH 21/40] krige: auto fitting not possible for spatio-temporal latlon models; update docs --- src/gstools/krige/base.py | 19 +++++++++++++---- src/gstools/krige/methods.py | 40 +++++++++++++++++++++++++++--------- 2 files changed, 45 insertions(+), 14 deletions(-) diff --git a/src/gstools/krige/base.py b/src/gstools/krige/base.py index 4d1ee824..774d2b58 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,6 +530,9 @@ 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.time: + 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) 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 """ From 868d4b113dd323cb317bac8e9ab0a9a4cb778731 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 17:09:56 +0200 Subject: [PATCH 22/40] CovModel: spatial vario/cov/cor now also use xyz with latlon models --- src/gstools/covmodel/base.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 8502be99..fe0fc3e9 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -580,7 +580,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 From 12012c694a73731606451286b0c50543b963f46a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Mon, 12 Jun 2023 17:10:13 +0200 Subject: [PATCH 23/40] plot: minor fixes for latlon --- src/gstools/covmodel/plot.py | 16 ++++++++-------- src/gstools/field/plot.py | 4 ++-- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/gstools/covmodel/plot.py b/src/gstools/covmodel/plot.py index fa526612..168dc1b4 100644 --- a/src/gstools/covmodel/plot.py +++ b/src/gstools/covmodel/plot.py @@ -52,12 +52,12 @@ # plotting routines ####################################################### -def _plot_spatial(dim, pos, field, fig, ax, latlon, **kwargs): +def _plot_spatial(dim, pos, field, fig, ax, time, **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_nd(pos, field, "structured", fig, ax, time=time, **kwargs) def plot_vario_spatial( @@ -70,7 +70,7 @@ 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.time, **kwargs) def plot_cov_spatial( @@ -83,7 +83,7 @@ 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.time, **kwargs) def plot_cor_spatial( @@ -96,7 +96,7 @@ 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.time, **kwargs) def plot_variogram( @@ -150,7 +150,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 = 3 * model.len_scale + 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 +165,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 = 3 * model.len_scale + 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 +180,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 = 3 * model.len_scale + 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/field/plot.py b/src/gstools/field/plot.py index f6242535..38d744ff 100644 --- a/src/gstools/field/plot.py +++ b/src/gstools/field/plot.py @@ -349,9 +349,9 @@ def _ax_names(dim, latlon=False, time=False, ax_names=None): return ax_names[:dim] if dim == 2 + t_fac and latlon: return ["lon", "lat"] + t_fac * ["time"] - if dim <= 3: + if dim - t_fac <= 3: return ( - ["$x$", "$y$", "$z$"][:dim] + ["$x$", "$y$", "$z$"][: dim - t_fac] + t_fac * ["time"] + (dim == 1) * ["field"] ) From db3e019d7c9c1e950bcad7289137fd268f09667e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 13 Jun 2023 17:47:12 +0200 Subject: [PATCH 24/40] CovModel: rename 'time' to 'temporal' --- examples/09_spatio_temporal/01_precip_1d.py | 2 +- examples/09_spatio_temporal/02_precip_2d.py | 2 +- .../03_geographic_coordinates.py | 4 +-- examples/09_spatio_temporal/README.rst | 32 +++++++++++++------ src/gstools/covmodel/base.py | 24 +++++++------- src/gstools/covmodel/plot.py | 18 ++++++++--- src/gstools/covmodel/tools.py | 16 ++++++---- src/gstools/field/base.py | 4 +-- src/gstools/field/plot.py | 16 +++++----- src/gstools/tools/geometric.py | 20 +++++++----- 10 files changed, 83 insertions(+), 55 deletions(-) diff --git a/examples/09_spatio_temporal/01_precip_1d.py b/examples/09_spatio_temporal/01_precip_1d.py index 2c431c7f..2da77783 100644 --- a/examples/09_spatio_temporal/01_precip_1d.py +++ b/examples/09_spatio_temporal/01_precip_1d.py @@ -30,7 +30,7 @@ st_anis = 0.4 # an exponential variogram with a corr. lengths of 2d and 5km -model = gs.Exponential(dim=1, time=True, var=1.0, len_scale=5.0, anis=st_anis) +model = gs.Exponential(dim=1, temporal=True, 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 d3d781b3..3406803b 100644 --- a/examples/09_spatio_temporal/02_precip_2d.py +++ b/examples/09_spatio_temporal/02_precip_2d.py @@ -31,7 +31,7 @@ st_anis = 0.4 # an exponential variogram with a corr. lengths of 5km, 5km, and 2d -model = gs.Exponential(dim=2, time=True, var=1.0, len_scale=5.0, anis=st_anis) +model = gs.Exponential(dim=2, temporal=True, 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 index 3ec8d28c..63ac29b0 100644 --- a/examples/09_spatio_temporal/03_geographic_coordinates.py +++ b/examples/09_spatio_temporal/03_geographic_coordinates.py @@ -5,7 +5,7 @@ In this example, we demonstrate how to generate a spatio-temporal random field on geographical coordinates. -First we setup a model, with ``latlon=True`` and ``time=True``, +First we setup a model, with ``latlon=True`` and ``temporal=True``, to get the associated spatio-temporal Yadrenko model. In addition, we will use the earth radius provided by :any:`EARTH_RADIUS` @@ -22,7 +22,7 @@ model = gs.Matern( latlon=True, - time=True, + temporal=True, var=1, len_scale=1000, anis=0.1, diff --git a/examples/09_spatio_temporal/README.rst b/examples/09_spatio_temporal/README.rst index 07aa5faa..08cc2fcd 100644 --- a/examples/09_spatio_temporal/README.rst +++ b/examples/09_spatio_temporal/README.rst @@ -5,28 +5,42 @@ 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 and setting a +spatio-temporal anisotropy ratio 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(dim=dim, temporal=True, 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. +There are 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 one 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(dim=dim, temporal=True, 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/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index fe0fc3e9..fad19671 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -111,7 +111,7 @@ class CovModel: Can be set to :any:`EARTH_RADIUS` to have len_scale in km or :any:`DEGREE_SCALE` to have len_scale in degree. Default: :any:`RADIAN_SCALE` - time : :class:`bool`, optional + 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`. @@ -147,7 +147,7 @@ def __init__( rescale=None, latlon=False, geo_scale=RADIAN_SCALE, - time=False, + temporal=False, var_raw=None, hankel_kw=None, **opt_arg, @@ -173,14 +173,14 @@ def __init__( self._nugget_bounds = None self._anis_bounds = None self._opt_arg_bounds = {} - # Set latlon and time first + # Set latlon and temporal first self._latlon = bool(latlon) - self._time = bool(time) + 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 - self.dim = dim + int(self.time) + self.dim = dim + int(self.temporal) # optional arguments for the variogram-model set_opt_args(self, opt_arg) @@ -198,7 +198,7 @@ def __init__( if self.latlon: # keep time anisotropy for metric spatio-temporal model self._anis = np.array((self.dim - 1) * [1], dtype=np.double) - self._anis[-1] = anis[-1] if self.time else 1.0 + self._anis[-1] = anis[-1] if self.temporal else 1.0 self._angles = np.array(self.dim * [0], dtype=np.double) else: self._anis = anis @@ -555,7 +555,7 @@ def isometrize(self, pos): return latlon2pos( pos, radius=self.geo_scale, - time=self.time, + temporal=self.temporal, time_scale=self.anis[-1], ) return np.dot(matrix_isometrize(self.dim, self.angles, self.anis), pos) @@ -567,7 +567,7 @@ def anisometrize(self, pos): return pos2latlon( pos, radius=self.geo_scale, - time=self.time, + temporal=self.temporal, time_scale=self.anis[-1], ) return np.dot( @@ -897,9 +897,9 @@ def arg_bounds(self): return res @property - def time(self): + def temporal(self): """:class:`bool`: Whether the model is a metric spatio-temporal one.""" - return self._time + return self._temporal # geographical coordinates related @@ -916,12 +916,12 @@ def geo_scale(self): @property def field_dim(self): """:class:`int`: The (parametric) field dimension of the model (with time).""" - return 2 + int(self._time) if self.latlon else self.dim + 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._time) + return 2 if self.latlon else self.dim - int(self.temporal) # standard parameters diff --git a/src/gstools/covmodel/plot.py b/src/gstools/covmodel/plot.py index 168dc1b4..43f94df6 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, time, **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, time=time, **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.time, **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.time, **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.time, **kwargs) + return _plot_spatial( + model.dim, pos, fld, fig, ax, model.temporal, **kwargs + ) def plot_variogram( diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index a1bef143..2c466dc5 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -498,13 +498,13 @@ def set_dim(model, dim): AttributeWarning, ) dim = model.fix_dim() - if model.latlon and dim != (3 + int(model.time)): + if model.latlon and dim != (3 + int(model.temporal)): raise ValueError( f"{model.name}: using fixed dimension {model.fix_dim()}, " - f"which is not compatible with a latlon model (with time={model.time})." + f"which is not compatible with a latlon model (with temporal={model.temporal})." ) - # force dim=3 (or 4 when time=True) for latlon models - dim = (3 + int(model.time)) 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.") @@ -551,7 +551,7 @@ def compare(this, that): equal &= np.all(np.isclose(this.angles, that.angles)) equal &= np.isclose(this.rescale, that.rescale) equal &= this.latlon == that.latlon - equal &= this.time == that.time + equal &= this.temporal == that.temporal for opt in this.opt_arg: equal &= np.isclose(getattr(this, opt), getattr(that, opt)) return equal @@ -569,14 +569,16 @@ def model_repr(model): # pragma: no cover m = model p = model._prec opt_str = "" - t_str = ", time=True" if m.time else "" + 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}}" if m.latlon: ani_str = ( - "" if m.is_isotropic or not m.time else f", anis={m.anis[-1]:.{p}}" + "" + if m.is_isotropic or not m.temporal + else f", anis={m.anis[-1]:.{p}}" ) r_str = ( "" diff --git a/src/gstools/field/base.py b/src/gstools/field/base.py index 6098e219..bb514141 100755 --- a/src/gstools/field/base.py +++ b/src/gstools/field/base.py @@ -679,9 +679,9 @@ def latlon(self): return False if self.model is None else self.model.latlon @property - def time(self): + def temporal(self): """:class:`bool`: Whether the field depends on time.""" - return False if self.model is None else self.model.time + return False if self.model is None else self.model.temporal @property def name(self): diff --git a/src/gstools/field/plot.py b/src/gstools/field/plot.py index 38d744ff..37af18bd 100644 --- a/src/gstools/field/plot.py +++ b/src/gstools/field/plot.py @@ -60,7 +60,7 @@ def plot_field( fig, ax, fld.latlon, - fld.time, + fld.temporal, **kwargs, ) @@ -111,7 +111,7 @@ def plot_nd( fig=None, ax=None, latlon=False, - time=False, + temporal=False, resolution=128, ax_names=None, aspect="quad", @@ -144,7 +144,7 @@ def plot_nd( ValueError will be raised, if a direction was specified. Bin edges need to be given in radians in this case. Default: False - time : :class:`bool`, optional + 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). @@ -172,20 +172,20 @@ def plot_nd( """ dim = len(pos) assert dim > 1 - assert not latlon or dim == 2 + int(bool(time)) + 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 ) if latlon: # swap lat-lon to lon-lat (x-y) - if time: + 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, time, ax_names) + 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] @@ -342,8 +342,8 @@ def plot_vec_field(fld, field="field", fig=None, ax=None): # pragma: no cover return ax -def _ax_names(dim, latlon=False, time=False, ax_names=None): - t_fac = int(bool(time)) +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] diff --git a/src/gstools/tools/geometric.py b/src/gstools/tools/geometric.py index fcbd7c68..7f2ea10e 100644 --- a/src/gstools/tools/geometric.py +++ b/src/gstools/tools/geometric.py @@ -626,7 +626,7 @@ def ang2dir(angles, dtype=np.double, dim=None): def latlon2pos( - latlon, radius=1.0, dtype=np.double, time=False, time_scale=1.0 + latlon, radius=1.0, dtype=np.double, temporal=False, time_scale=1.0 ): """Convert lat-lon geo coordinates to 3D position tuple. @@ -641,7 +641,7 @@ def latlon2pos( The desired data-type for the array. If not given, then the type will be determined as the minimum type required to hold the objects in the sequence. Default: None - time : :class:`bool`, optional + temporal : :class:`bool`, optional Whether latlon includes an appended time axis. Default: False time_scale : :class:`float`, optional @@ -653,19 +653,23 @@ def latlon2pos( :class:`numpy.ndarray` the 3D position array """ - latlon = np.asarray(latlon, dtype=dtype).reshape((3 if time else 2, -1)) + 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 time: + 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, time=False, time_scale=1.0): +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 @@ -679,7 +683,7 @@ def pos2latlon(pos, radius=1.0, dtype=np.double, time=False, time_scale=1.0): The desired data-type for the array. If not given, then the type will be determined as the minimum type required to hold the objects in the sequence. Default: None - time : :class:`bool`, optional + temporal : :class:`bool`, optional Whether latlon includes an appended time axis. Default: False time_scale : :class:`float`, optional @@ -691,12 +695,12 @@ def pos2latlon(pos, radius=1.0, dtype=np.double, time=False, time_scale=1.0): :class:`numpy.ndarray` the 3D position array """ - pos = np.asarray(pos, dtype=dtype).reshape((4 if time else 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]) latlon = np.rad2deg((lat, lon), dtype=dtype) - if time: + if temporal: return np.array( (latlon[0], latlon[1], pos[3] * time_scale), dtype=dtype ) From 2f97808e96f14d253fb319b493da168b67c01adf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 13 Jun 2023 18:10:05 +0200 Subject: [PATCH 25/40] geo_scale: better docs; always use KM_SCALE in examples --- README.md | 2 +- docs/source/index.rst | 2 +- examples/03_variogram/06_auto_bin_latlon.py | 4 ++-- examples/08_geo_coordinates/00_field_generation.py | 8 +++++--- examples/08_geo_coordinates/01_dwd_krige.py | 4 ++-- examples/08_geo_coordinates/README.rst | 2 +- .../09_spatio_temporal/03_geographic_coordinates.py | 11 ++++++----- src/gstools/covmodel/base.py | 9 +++++---- src/gstools/variogram/binning.py | 9 +++++---- src/gstools/variogram/variogram.py | 9 +++++---- tests/test_latlon.py | 8 ++++---- 11 files changed, 37 insertions(+), 31 deletions(-) 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 Date: Tue, 13 Jun 2023 18:24:07 +0200 Subject: [PATCH 26/40] Plot: minor fixes for st plots --- src/gstools/covmodel/plot.py | 2 +- src/gstools/field/plot.py | 26 +++++++++++++++++++++----- 2 files changed, 22 insertions(+), 6 deletions(-) diff --git a/src/gstools/covmodel/plot.py b/src/gstools/covmodel/plot.py index 43f94df6..334ccbe2 100644 --- a/src/gstools/covmodel/plot.py +++ b/src/gstools/covmodel/plot.py @@ -56,7 +56,7 @@ 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_1d(pos, field, fig, ax, temporal, **kwargs) return plot_nd( pos, field, "structured", fig, ax, temporal=temporal, **kwargs ) diff --git a/src/gstools/field/plot.py b/src/gstools/field/plot.py index 37af18bd..d06a22a1 100644 --- a/src/gstools/field/plot.py +++ b/src/gstools/field/plot.py @@ -52,7 +52,7 @@ 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], @@ -65,7 +65,9 @@ def plot_field( ) -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. @@ -76,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` @@ -95,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]) @@ -175,7 +182,15 @@ def plot_nd( 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, ) if latlon: # swap lat-lon to lon-lat (x-y) @@ -365,6 +380,7 @@ def _plot_2d( fig=None, ax=None, latlon=False, + temporal=False, ax_names=None, levels=64, antialias=True, @@ -372,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) From b62550cc58526438de7e80a0ef0754c05b31a40c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 13 Jun 2023 18:25:32 +0200 Subject: [PATCH 27/40] Krige: bugfix for renamed attribute temporal --- src/gstools/krige/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gstools/krige/base.py b/src/gstools/krige/base.py index 774d2b58..336f4cc0 100755 --- a/src/gstools/krige/base.py +++ b/src/gstools/krige/base.py @@ -530,7 +530,7 @@ 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.time: + 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) From 1e47ab773058defb34dd34622e887742e9866cd9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 13 Jun 2023 19:14:26 +0200 Subject: [PATCH 28/40] CovModel: prevent rotation between spatial and temporal dims --- src/gstools/covmodel/base.py | 35 ++++++++++------------- src/gstools/covmodel/tools.py | 54 +++++++++++++++++++++++++++++++++-- 2 files changed, 66 insertions(+), 23 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index f5d90d14..7ed5960d 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -28,6 +28,7 @@ set_arg_bounds, set_dim, set_len_anis, + set_model_angles, set_opt_args, spectral_rad_pdf, ) @@ -39,7 +40,6 @@ matrix_isometrize, pos2latlon, rotated_main_axes, - set_angles, ) __all__ = ["CovModel"] @@ -194,16 +194,15 @@ def __init__( # set parameters self.rescale = rescale self._nugget = float(nugget) + # set anisotropy and len_scale, disable anisotropy for latlon models - self._len_scale, anis = set_len_anis(self.dim, len_scale, anis) - if self.latlon: - # keep time anisotropy for metric spatio-temporal model - self._anis = np.array((self.dim - 1) * [1], dtype=np.double) - self._anis[-1] = anis[-1] if self.temporal else 1.0 - self._angles = np.array(self.dim * [0], dtype=np.double) - else: - self._anis = anis - self._angles = set_angles(self.dim, angles) + self._len_scale, self._anis = set_len_anis( + self.dim, len_scale, anis, self.latlon, self.temporal + ) + self._angles = set_model_angles( + self.dim, angles, self.latlon, self.temporal + ) + # set var at last, because of the var_factor (to be right initialized) if var_raw is None: self._var = None @@ -1004,12 +1003,9 @@ def anis(self): @anis.setter def anis(self, anis): - if self.latlon: - self._anis = np.array((self.dim - 1) * [1], dtype=np.double) - else: - self._len_scale, self._anis = set_len_anis( - self.dim, self.len_scale, anis - ) + self._len_scale, self._anis = set_len_anis( + self.dim, self.len_scale, anis, self.latlon, self.temporal + ) self.check_arg_bounds() @property @@ -1019,10 +1015,9 @@ def angles(self): @angles.setter def angles(self, angles): - if self.latlon: - self._angles = np.array(self.dim * [0], dtype=np.double) - else: - self._angles = set_angles(self.dim, angles) + self._angles = set_model_angles( + self.dim, angles, self.latlon, self.temporal + ) self.check_arg_bounds() @property diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 2c466dc5..17c76531 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, temporal=False): """Set the length scale and anisotropy factors for the given dimension. Parameters @@ -194,6 +195,13 @@ def set_len_anis(dim, len_scale, anis): the length scale of the SRF in x direction or in x- (y-, ...) direction anis : :class:`float` or :class:`list` the anisotropy of length scales along the transversal axes + latlon : :class:`bool`, optional + Whether the model is describing 2D fields on earths surface described + by latitude and longitude. + Default: False + temporal : :class:`bool`, optional + Whether a time-dimension is appended. + Default: False Returns ------- @@ -235,9 +243,47 @@ def set_len_anis(dim, len_scale, anis): raise ValueError( "anisotropy-ratios needs to be > 0, " + "got: " + str(out_anis) ) + # no anisotropy for latlon (only when temporal, but then only in time-dimension) + if latlon: + out_anis[: dim - (2 if temporal else 1)] = 1.0 return out_len_scale, out_anis +def set_model_angles(dim, angles, latlon=False, temporal=False): + """Set the model angles for the given dimension. + + Parameters + ---------- + dim : :class:`int` + spatial dimension + angles : :class:`float` or :class:`list` + the angles of the SRF + latlon : :class:`bool`, optional + Whether the model is describing 2D fields on earths surface described + by latitude and longitude. + Default: False + temporal : :class:`bool`, optional + Whether a time-dimension is appended. + Default: False + + Returns + ------- + angles : :class:`float` + the angles fitting to the dimension + + Notes + ----- + If too few angles are given, they are filled up with `0`. + """ + if latlon: + return np.array(no_of_angles(dim) * [0], dtype=np.double) + out_angles = set_angles(dim, angles) + if temporal: + # no rotation between spatial dimensions and temporal dimension + out_angles[no_of_angles(dim - 1) :] = 0.0 + return out_angles + + def check_bounds(bounds): """ Check if given bounds are valid. @@ -522,7 +568,9 @@ def set_dim(model, dim): model.dim, model._len_scale, model._anis ) if model._angles is not None: - model._angles = set_angles(model.dim, model._angles) + model._angles = set_model_angles( + model.dim, model._angles, model.latlon, model.temporal + ) model.check_arg_bounds() From e821a6a423e79c207d477cb8c834fa8180667bd5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 13 Jun 2023 19:24:00 +0200 Subject: [PATCH 29/40] CovModel: saver setting of anis --- src/gstools/covmodel/base.py | 8 +++++--- src/gstools/covmodel/tools.py | 13 +++++-------- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 7ed5960d..b6a533f7 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -197,7 +197,7 @@ def __init__( # 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.temporal + self.dim, len_scale, anis, self.latlon ) self._angles = set_model_angles( self.dim, angles, self.latlon, self.temporal @@ -974,7 +974,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: @@ -1004,7 +1006,7 @@ def anis(self): @anis.setter def anis(self, anis): self._len_scale, self._anis = set_len_anis( - self.dim, self.len_scale, anis, self.latlon, self.temporal + self.dim, self.len_scale, anis, self.latlon ) self.check_arg_bounds() diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 17c76531..a029512f 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -184,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, latlon=False, temporal=False): +def set_len_anis(dim, len_scale, anis, latlon=False): """Set the length scale and anisotropy factors for the given dimension. Parameters @@ -197,10 +197,7 @@ def set_len_anis(dim, len_scale, anis, latlon=False, temporal=False): the anisotropy of length scales along the transversal axes latlon : :class:`bool`, optional Whether the model is describing 2D fields on earths surface described - by latitude and longitude. - Default: False - temporal : :class:`bool`, optional - Whether a time-dimension is appended. + by latitude and longitude. In this case there is no spatial anisotropy. Default: False Returns @@ -241,11 +238,11 @@ def set_len_anis(dim, len_scale, anis, latlon=False, temporal=False): 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 anisotropy for latlon (only when temporal, but then only in time-dimension) + # no spatial anisotropy for latlon if latlon: - out_anis[: dim - (2 if temporal else 1)] = 1.0 + out_anis[:2] = 1.0 return out_len_scale, out_anis From 677fd1a1111b075c2514c9b08b1c4ffdc41a43a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Tue, 13 Jun 2023 19:27:27 +0200 Subject: [PATCH 30/40] minor f-string fix --- src/gstools/covmodel/tools.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index a029512f..8fcca7a9 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -238,7 +238,7 @@ def set_len_anis(dim, len_scale, anis, latlon=False): for ani in out_anis: if not ani > 0.0: raise ValueError( - f"anisotropy-ratios needs to be > 0, " + "got: {out_anis}" + f"anisotropy-ratios needs to be > 0, got: {out_anis}" ) # no spatial anisotropy for latlon if latlon: From f4bc0d1f6fe595b16d9155523bc5cd783101d4a8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 14 Jun 2023 12:01:11 +0200 Subject: [PATCH 31/40] test temporal related stuff --- tests/test_latlon.py | 2 +- tests/test_temporal.py | 43 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 1 deletion(-) create mode 100644 tests/test_temporal.py diff --git a/tests/test_latlon.py b/tests/test_latlon.py index 05cb9176..bd347281 100644 --- a/tests/test_latlon.py +++ b/tests/test_latlon.py @@ -144,7 +144,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..800743e5 --- /dev/null +++ b/tests/test_temporal.py @@ -0,0 +1,43 @@ +""" +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): + ... + + 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_rotation(self): + mod = gs.Gaussian(dim=3, temporal=True, angles=[1, 2, 3, 4, 5, 6]) + self.assertTrue(np.allclose(mod.angles, [1, 2, 3, 0, 0, 0])) + + def test_krige(self): + mod = gs.Gaussian(latlon=True, temporal=True) + # auto-fitting latlon-temporal model in kriging not possible + with self.assertRaises(ValueError): + kri = gs.Krige(mod, 3 * [[1, 2, 3]], [1, 2, 3], fit_variogram=True) + + +if __name__ == "__main__": + unittest.main() From fad4afe85b89ec464070831147ab2e8b846cd05e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 14 Jun 2023 13:46:32 +0200 Subject: [PATCH 32/40] more temporal tests --- tests/test_latlon.py | 16 ++++++++++------ tests/test_temporal.py | 38 +++++++++++++++++++++++++++++++++++--- 2 files changed, 45 insertions(+), 9 deletions(-) diff --git a/tests/test_latlon.py b/tests/test_latlon.py index bd347281..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.KM_SCALE + 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.KM_SCALE) + 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.KM_SCALE) + 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.KM_SCALE) + 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] diff --git a/tests/test_temporal.py b/tests/test_temporal.py index 800743e5..98893693 100644 --- a/tests/test_temporal.py +++ b/tests/test_temporal.py @@ -11,7 +11,13 @@ 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( @@ -28,15 +34,41 @@ def test_latlon(self): 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(dim=3, temporal=True, angles=[1, 2, 3, 4, 5, 6]) self.assertTrue(np.allclose(mod.angles, [1, 2, 3, 0, 0, 0])) def test_krige(self): - mod = gs.Gaussian(latlon=True, temporal=True) # auto-fitting latlon-temporal model in kriging not possible with self.assertRaises(ValueError): - kri = gs.Krige(mod, 3 * [[1, 2, 3]], [1, 2, 3], fit_variogram=True) + 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__": From 44e174c327596ff1fc7adf8946cfa324e453ae78 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Wed, 14 Jun 2023 16:36:14 +0200 Subject: [PATCH 33/40] CovModel: add 'spatial_dim' argument --- examples/09_spatio_temporal/01_precip_1d.py | 4 +++- examples/09_spatio_temporal/02_precip_2d.py | 4 +++- examples/09_spatio_temporal/README.rst | 17 +++++++++++------ src/gstools/covmodel/base.py | 17 ++++++++++++++--- tests/test_temporal.py | 5 ++++- 5 files changed, 35 insertions(+), 12 deletions(-) diff --git a/examples/09_spatio_temporal/01_precip_1d.py b/examples/09_spatio_temporal/01_precip_1d.py index 2da77783..4b4c6b8a 100644 --- a/examples/09_spatio_temporal/01_precip_1d.py +++ b/examples/09_spatio_temporal/01_precip_1d.py @@ -30,7 +30,9 @@ st_anis = 0.4 # an exponential variogram with a corr. lengths of 2d and 5km -model = gs.Exponential(dim=1, temporal=True, var=1, len_scale=5, 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 3406803b..81c78964 100644 --- a/examples/09_spatio_temporal/02_precip_2d.py +++ b/examples/09_spatio_temporal/02_precip_2d.py @@ -31,7 +31,9 @@ st_anis = 0.4 # an exponential variogram with a corr. lengths of 5km, 5km, and 2d -model = gs.Exponential(dim=2, temporal=True, var=1, len_scale=5, 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/README.rst b/examples/09_spatio_temporal/README.rst index 08cc2fcd..3cb06b9e 100644 --- a/examples/09_spatio_temporal/README.rst +++ b/examples/09_spatio_temporal/README.rst @@ -6,23 +6,28 @@ like rainfall, air temperature or crop yield. GSTools provides the metric spatio-temporal model for all covariance models by setting ``temporal=True``, which enhances the spatial model dimension with -a time dimension to result in the spatio-temporal dimension and setting a -spatio-temporal anisotropy ratio like this: +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_anis = 0.4 - st_model = gs.Exponential(dim=dim, temporal=True, anis=st_anis) + st_model = gs.Exponential(temporal=True, spatial_dim=dim, anis=st_anis) Since it is given in the name "spatio-temporal", time is always treated as last dimension. -There are three different dimension attributes giving information about (i) the +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 one with geographic coordinates, the model dimension is 4, +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: @@ -40,7 +45,7 @@ non-temporal models, without altering the behavior in the time dimension: 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=dim, temporal=True, 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/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index b6a533f7..077edaa1 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -54,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`` @@ -118,6 +121,11 @@ class CovModel: `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. @@ -149,6 +157,7 @@ def __init__( latlon=False, geo_scale=RADIAN_SCALE, temporal=False, + spatial_dim=None, var_raw=None, hankel_kw=None, **opt_arg, @@ -180,8 +189,10 @@ def __init__( 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 - self.dim = dim + int(self.temporal) + # 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) diff --git a/tests/test_temporal.py b/tests/test_temporal.py index 98893693..c179db1c 100644 --- a/tests/test_temporal.py +++ b/tests/test_temporal.py @@ -58,8 +58,11 @@ def test_latlon2pos(self): ) def test_rotation(self): - mod = gs.Gaussian(dim=3, temporal=True, angles=[1, 2, 3, 4, 5, 6]) + 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 From 0ece19c9a09bb85a76eb0c066f706acb4131b5ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 15 Jun 2023 10:40:51 +0200 Subject: [PATCH 34/40] minor f-string fixes --- src/gstools/covmodel/tools.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 8fcca7a9..dddeb441 100644 --- a/src/gstools/covmodel/tools.py +++ b/src/gstools/covmodel/tools.py @@ -421,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): @@ -537,7 +535,7 @@ 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() From 02806be571087e492d92c83c07a0397cb272290e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 15 Jun 2023 10:41:07 +0200 Subject: [PATCH 35/40] update changelog --- CHANGELOG.md | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 61aa289f..7b956b3a 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 +- CovModels expect all arguments by keyword now (except `dim`) [#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 From d44c819570b10168ee5be45a6e986196bce499ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 15 Jun 2023 10:46:12 +0200 Subject: [PATCH 36/40] CovModel: be less strict about key-word-only args (don't want to bother people) --- CHANGELOG.md | 2 +- src/gstools/covmodel/base.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7b956b3a..070e0a71 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,7 +24,7 @@ All notable changes to **GSTools** will be documented in this file. - `vario_estimate` now forwards additional kwargs to `standard_bins` (`bin_no`, `max_dist`) [#308](https://github.com/GeoStat-Framework/GSTools/pull/308) ### Changes -- CovModels expect all arguments by keyword now (except `dim`) [#308](https://github.com/GeoStat-Framework/GSTools/pull/308) +- CovModels 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 diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 077edaa1..be530075 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -146,12 +146,12 @@ class CovModel: 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, From 8efe29dd2e3bc8ab124e1559885d924b50fb910d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 15 Jun 2023 17:29:45 +0200 Subject: [PATCH 37/40] variogram: rename bin_center to bin_centers --- src/gstools/variogram/variogram.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/gstools/variogram/variogram.py b/src/gstools/variogram/variogram.py index afb03e70..b1edddf4 100644 --- a/src/gstools/variogram/variogram.py +++ b/src/gstools/variogram/variogram.py @@ -238,7 +238,7 @@ def vario_estimate( 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. @@ -261,20 +261,17 @@ def vario_estimate( "Geostatistics for environmental scientists.", John Wiley & Sons. (2007) """ - if bin_edges is not None: - bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double, copy=False) - bin_center = (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_center = np.empty(0) if bin_edges is None else bin_center - estimates = np.zeros_like(bin_center) + bin_centers = np.empty(0) if bin_edges is None else bin_centers + estimates = np.zeros_like(bin_centers) if return_counts: - return bin_center, estimates, np.zeros_like(estimates, dtype=int) - return bin_center, estimates + return bin_centers, estimates, np.zeros_like(estimates, dtype=int) + return bin_centers, estimates if not masked: field = field.filled() # check mesh shape @@ -344,12 +341,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, geo_scale=geo_scale, **std_bins ) - bin_center = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + else: + bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double, copy=False) + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + if latlon: # internally we always use radians bin_edges /= geo_scale @@ -389,7 +389,7 @@ def vario_estimate( if dir_no == 1: estimates, counts = estimates[0], counts[0] est_out = (estimates, counts) - return (bin_center,) + 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( From 6356ce0bd9df6ba28b7660a883a42c26f4c2bb95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 15 Jun 2023 17:30:44 +0200 Subject: [PATCH 38/40] changelog: minor fix --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 070e0a71..477fc00f 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -24,7 +24,7 @@ All notable changes to **GSTools** will be documented in this file. - `vario_estimate` now forwards additional kwargs to `standard_bins` (`bin_no`, `max_dist`) [#308](https://github.com/GeoStat-Framework/GSTools/pull/308) ### Changes -- CovModels expect special arguments by keyword now [#308](https://github.com/GeoStat-Framework/GSTools/pull/308) +- `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 From e72bfefc1953fcbe85c1db90bb4c0ef98967302c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 15 Jun 2023 17:38:41 +0200 Subject: [PATCH 39/40] vario: revert moving code-block --- src/gstools/variogram/variogram.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/gstools/variogram/variogram.py b/src/gstools/variogram/variogram.py index b1edddf4..69746658 100644 --- a/src/gstools/variogram/variogram.py +++ b/src/gstools/variogram/variogram.py @@ -261,6 +261,9 @@ def vario_estimate( "Geostatistics for environmental scientists.", John Wiley & Sons. (2007) """ + if bin_edges is not None: + bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double, copy=False) + 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) @@ -346,10 +349,7 @@ def vario_estimate( bin_edges = standard_bins( pos, dim, latlon, geo_scale=geo_scale, **std_bins ) - else: - bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double, copy=False) - bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0 - + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.0 if latlon: # internally we always use radians bin_edges /= geo_scale From 0881384cb3cc22dcd768a6a8a4a0ee42bee8c41a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Thu, 15 Jun 2023 17:41:25 +0200 Subject: [PATCH 40/40] changelog: minor markdown fixes --- CHANGELOG.md | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 477fc00f..733d75e3 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,16 +5,16 @@ 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) +- 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 + - 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) +- 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 @@ -53,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) @@ -209,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)): @@ -236,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** @@ -260,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 @@ -276,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)) @@ -294,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 @@ -307,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) @@ -325,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 @@ -333,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 @@ -381,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