diff --git a/examples/03_variogram/06_auto_bin_latlon.py b/examples/03_variogram/06_auto_bin_latlon.py index 70c0a09b..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) -sph = gs.Spherical(latlon=True, rescale=gs.EARTH_RADIUS) +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) diff --git a/examples/08_geo_coordinates/00_field_generation.py b/examples/08_geo_coordinates/00_field_generation.py index b5685c7d..83700d61 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 numpy as np + import gstools as gs -model = gs.Gaussian(latlon=True, var=1, len_scale=777, rescale=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) @@ -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 b37e7fa0..98c570b4 100755 --- a/examples/08_geo_coordinates/01_dwd_krige.py +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -76,17 +76,17 @@ 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. # 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 ``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,10 +97,10 @@ 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.fit_variogram(bin_c, vario, nugget=False) -ax = model.plot("vario_yadrenko", x_max=bins[-1]) -ax.scatter(bin_c, vario) +model = gs.Spherical(latlon=True, geo_scale=gs.EARTH_RADIUS) +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) ############################################################################### diff --git a/examples/08_geo_coordinates/README.rst b/examples/08_geo_coordinates/README.rst index 87b419df..b664512f 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 ``geo_scale`` argument: .. code-block:: python import gstools as gs - model = gs.Gaussian(latlon=True, var=2, len_scale=500, rescale=gs.EARTH_RADIUS) + model = gs.Gaussian(latlon=True, var=2, len_scale=500, geo_scale=gs.EARTH_RADIUS) Then ``len_scale`` can be interpreted as given in km. @@ -37,20 +37,21 @@ A `Yadrenko` model :math:`C` is derived from a valid isotropic covariance model in 3D :math:`C_{3D}` by the following relation: .. math:: - C(\zeta)=C_{3D}\left(2 \cdot \sin\left(\frac{\zeta}{2}\right)\right) + C(\zeta)=C_{3D}\left(2r \cdot \sin\left(\frac{\zeta}{2r}\right)\right) Where :math:`\zeta` is the -`great-circle distance `_. +`great-circle distance `_ +and :math:`r` is the ``geo_scale``. .. note:: ``lat`` and ``lon`` are given in degree, whereas the great-circle distance - :math:`zeta` is given in radians. + :math:`zeta` is given in units of the ``geo_scale``. -Note, that :math:`2 \cdot \sin(\frac{\zeta}{2})` is the +Note, that :math:`2r \cdot \sin(\frac{\zeta}{2r})` is the `chordal distance `_ -of two points on a sphere, which means we simply think of the earth surface -as a sphere, that is cut out of the surrounding three dimensional space, +of two points on a sphere with radius :math:`r`, which means we simply think of the +earth surface as a sphere, that is cut out of the surrounding three dimensional space, when using the `Yadrenko` model. .. note:: diff --git a/examples/09_spatio_temporal/01_precip_1d.py b/examples/09_spatio_temporal/01_precip_1d.py index 4671b2f7..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_geographic_coordinates.py b/examples/09_spatio_temporal/03_geographic_coordinates.py new file mode 100644 index 00000000..3ec8d28c --- /dev/null +++ b/examples/09_spatio_temporal/03_geographic_coordinates.py @@ -0,0 +1,36 @@ +""" +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` +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 `100` days. +""" +import numpy as np + +import gstools as gs + +model = gs.Matern( + latlon=True, + time=True, + var=1, + len_scale=1000, + anis=0.1, + geo_scale=gs.EARTH_RADIUS, +) + +lat = lon = np.linspace(-80, 81, 50) +time = np.linspace(0, 777, 50) +srf = gs.SRF(model, seed=1234) +field = srf.structured((lat, lon, time)) +srf.plot() diff --git a/pyproject.toml b/pyproject.toml index 41090458..3a6c7ec0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -147,9 +147,9 @@ target-version = [ max-args = 20 max-locals = 50 max-branches = 30 - max-statements = 80 + max-statements = 85 max-attributes = 25 - max-public-methods = 75 + max-public-methods = 80 [tool.cibuildwheel] # Switch to using build diff --git a/src/gstools/__init__.py b/src/gstools/__init__.py index 60dd1e33..6d62d558 100644 --- a/src/gstools/__init__.py +++ b/src/gstools/__init__.py @@ -125,6 +125,9 @@ .. autosummary:: EARTH_RADIUS + KM_SCALE + DEGREE_SCALE + RADIAN_SCALE """ # Hooray! from gstools import ( @@ -161,7 +164,10 @@ from gstools.field import SRF, CondSRF from gstools.krige import Krige from gstools.tools import ( + DEGREE_SCALE, EARTH_RADIUS, + KM_SCALE, + RADIAN_SCALE, generate_grid, generate_st_grid, rotated_main_axes, @@ -188,7 +194,7 @@ __all__ = ["__version__"] __all__ += ["covmodel", "field", "variogram", "krige", "random", "tools"] -__all__ += ["transform", "normalizer"] +__all__ += ["transform", "normalizer", "config"] __all__ += [ "CovModel", "Gaussian", @@ -226,6 +232,9 @@ "generate_grid", "generate_st_grid", "EARTH_RADIUS", + "KM_SCALE", + "DEGREE_SCALE", + "RADIAN_SCALE", "vtk_export", "vtk_export_structured", "vtk_export_unstructured", diff --git a/src/gstools/covmodel/base.py b/src/gstools/covmodel/base.py index 686768fa..fe0fc3e9 100644 --- a/src/gstools/covmodel/base.py +++ b/src/gstools/covmodel/base.py @@ -31,7 +31,9 @@ set_opt_args, spectral_rad_pdf, ) +from gstools.tools import RADIAN_SCALE from gstools.tools.geometric import ( + great_circle_to_chordal, latlon2pos, matrix_anisometrize, matrix_isometrize, @@ -100,9 +102,21 @@ class CovModel: :math:`2\sin(\alpha/2)`, where :math:`\alpha` is the great-circle distance, which is equal to the spatial distance of two points in 3D. As a consequence, `dim` will be set to `3` and anisotropy will be - disabled. `rescale` can be set to e.g. earth's radius, + disabled. `geo_scale` can be set to e.g. earth's radius, to have a meaningful `len_scale` parameter. Default: False + geo_scale : :class:`float`, optional + Geographic 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. + 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. + `spatial_dim` will be `field_dim - 1`. + The time-dimension is appended, meaning the pos tuple is (x,y,z,...,t). + Default: False var_raw : :class:`float` or :any:`None`, optional raw variance of the model which will be multiplied with :any:`CovModel.var_factor` to result in the actual variance. @@ -123,6 +137,7 @@ class CovModel: def __init__( self, dim=3, + *, var=1.0, len_scale=1.0, nugget=0.0, @@ -131,6 +146,8 @@ def __init__( integral_scale=None, rescale=None, latlon=False, + geo_scale=RADIAN_SCALE, + time=False, var_raw=None, hankel_kw=None, **opt_arg, @@ -156,11 +173,14 @@ 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) + self._geo_scale = abs(float(geo_scale)) # SFT class will be created within dim.setter but needs hankel_kw self.hankel_kw = hankel_kw - self.dim = dim + # using time increases model dimension + self.dim = dim + int(self.time) # optional arguments for the variogram-model set_opt_args(self, opt_arg) @@ -176,7 +196,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 @@ -242,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)) + return self.variogram(great_circle_to_chordal(zeta, self.geo_scale)) def cov_yadrenko(self, zeta): r"""Yadrenko covariance for great-circle distance from latlon-pos.""" - return self.covariance(2 * np.sin(zeta / 2)) + return self.covariance(great_circle_to_chordal(zeta, self.geo_scale)) def cor_yadrenko(self, zeta): r"""Yadrenko correlation for great-circle distance from latlon-pos.""" - return self.correlation(2 * np.sin(zeta / 2)) + return self.correlation(great_circle_to_chordal(zeta, self.geo_scale)) def vario_spatial(self, pos): r"""Spatial variogram respecting anisotropy and rotation.""" @@ -530,14 +552,24 @@ def isometrize(self, pos): """Make a position tuple ready for isotropic operations.""" pos = np.asarray(pos, dtype=np.double).reshape((self.field_dim, -1)) if self.latlon: - return latlon2pos(pos) + return latlon2pos( + pos, + radius=self.geo_scale, + 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) + return pos2latlon( + pos, + radius=self.geo_scale, + time=self.time, + time_scale=self.anis[-1], + ) return np.dot( matrix_anisometrize(self.dim, self.angles, self.anis), pos ) @@ -548,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 @@ -862,6 +896,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 @@ -869,10 +908,20 @@ def latlon(self): """:class:`bool`: Whether the model depends on geographical coords.""" return self._latlon + @property + def geo_scale(self): + """:class:`float`: Geographic scaling for geographical coords.""" + return self._geo_scale + @property def field_dim(self): - """:class:`int`: The field dimension of the model.""" - return 2 if self.latlon else self.dim + """:class:`int`: The (parametric) field dimension of the model (with time).""" + return 2 + int(self._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..dc2d5a3a 100755 --- a/src/gstools/covmodel/fit.py +++ b/src/gstools/covmodel/fit.py @@ -13,7 +13,7 @@ from scipy.optimize import curve_fit from gstools.covmodel.tools import check_arg_in_bounds, default_arg_from_bounds -from gstools.tools.geometric import set_anis +from gstools.tools.geometric import great_circle_to_chordal, set_anis __all__ = ["fit_variogram"] @@ -341,7 +341,7 @@ def _check_vario(model, x_data, y_data): ) if model.latlon: # convert to yadrenko model - x_data = 2 * np.sin(x_data / 2) + x_data = great_circle_to_chordal(x_data, model.geo_scale) return x_data, y_data, is_dir_vario @@ -522,6 +522,7 @@ def logistic_weights(p=0.1, mean=0.7): # pragma: no cover callable Weighting function. """ + # define the callable weights function def func(x_data): """Callable function for the weights.""" diff --git a/src/gstools/covmodel/plot.py b/src/gstools/covmodel/plot.py index efcc2630..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 = min(3 * model.len_rescaled, np.pi) + x_max = min(3 * model.len_scale, model.geo_scale * np.pi) x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} Yadrenko variogram") ax.plot(x_s, model.vario_yadrenko(x_s), **kwargs) @@ -165,7 +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.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_rescaled, np.pi) + x_max = min(3 * model.len_scale, model.geo_scale * np.pi) x_s = np.linspace(x_min, x_max) kwargs.setdefault("label", f"{model.name} Yadrenko correlation") ax.plot(x_s, model.cor_yadrenko(x_s), **kwargs) diff --git a/src/gstools/covmodel/tools.py b/src/gstools/covmodel/tools.py index 98ed3b8a..a1bef143 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 @@ -568,21 +569,33 @@ 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.geo_scale, 1) + else f", geo_scale={m.geo_scale:.{p}}" + ) repr_str = ( - f"{m.name}(latlon={m.latlon}, var={m.var:.{p}}, " - f"len_scale={m.len_scale:.{p}}, nugget={m.nugget:.{p}}{opt_str})" + f"{m.name}(latlon={m.latlon}{t_str}, var={m.var:.{p}}, " + f"len_scale={m.len_scale:.{p}}, nugget={m.nugget:.{p}}" + f"{ani_str}{r_str}{opt_str})" ) else: + # only print anis and angles if model is anisotropic or rotated + ani_str = "" if m.is_isotropic else f", anis={list_format(m.anis, p)}" + ang_str = ( + f", angles={list_format(m.angles, p)}" if m.do_rotation else "" + ) repr_str = ( - f"{m.name}(dim={m.dim}, var={m.var:.{p}}, " + f"{m.name}(dim={m.spatial_dim}{t_str}, var={m.var:.{p}}, " f"len_scale={m.len_scale:.{p}}, nugget={m.nugget:.{p}}" f"{ani_str}{ang_str}{opt_str})" ) diff --git a/src/gstools/field/base.py b/src/gstools/field/base.py index 903f3893..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..38d744ff 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 <= 3: - return ["$x$", "$y$", "$z$"][:dim] + (dim == 1) * ["field"] - return [f"$x_{{{i}}}$" for i in range(dim)] + if dim == 2 + t_fac and latlon: + return ["lon", "lat"] + t_fac * ["time"] + if dim - t_fac <= 3: + return ( + ["$x$", "$y$", "$z$"][: dim - t_fac] + + t_fac * ["time"] + + (dim == 1) * ["field"] + ) + return [f"$x_{{{i}}}$" for i in range(dim - t_fac)] + t_fac * ["time"] def _plot_2d( diff --git a/src/gstools/krige/base.py b/src/gstools/krige/base.py index 77c6832e..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,12 +530,18 @@ def set_condition( self.normalizer.fit(self.cond_val - self.cond_trend) if fit_variogram: # fitting model to empirical variogram of data # normalize field + if self.model.latlon and self.model.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) if self.model.is_isotropic: emp_vario = vario_estimate( - self.cond_pos, field, latlon=self.model.latlon + self.cond_pos, + field, + latlon=self.model.latlon, + geo_scale=self.model.geo_scale, ) else: axes = rotated_main_axes(self.model.dim, self.model.angles) diff --git a/src/gstools/krige/methods.py b/src/gstools/krige/methods.py index 653785f8..b258a02d 100644 --- a/src/gstools/krige/methods.py +++ b/src/gstools/krige/methods.py @@ -77,8 +77,12 @@ class Simple(Krige): Default: False fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -171,8 +175,12 @@ class Ordinary(Krige): Default: False fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -275,8 +283,12 @@ class Universal(Krige): Default: False fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -376,8 +388,12 @@ class ExtDrift(Krige): Default: False fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ @@ -467,8 +483,12 @@ class Detrended(Krige): Default: `"pinv"` fit_variogram : :class:`bool`, optional Whether to fit the given variogram model to the data. - This is done by using isotropy settings of the given model, - assuming the sill to be the data variance and with the + Directional variogram fitting is triggered by setting + any anisotropy factor of the model to anything unequal 1 + but the main axes of correlation are taken from the model + rotation angles. If the model is a spatio-temporal latlon + model, this will raise an error. + This assumes the sill to be the data variance and with standard bins provided by the :any:`standard_bins` routine. Default: False """ diff --git a/src/gstools/tools/__init__.py b/src/gstools/tools/__init__.py index bd329576..1f68dbaf 100644 --- a/src/gstools/tools/__init__.py +++ b/src/gstools/tools/__init__.py @@ -58,10 +58,19 @@ .. autosummary:: EARTH_RADIUS + KM_SCALE + DEGREE_SCALE + RADIAN_SCALE ---- .. autodata:: EARTH_RADIUS + +.. autodata:: KM_SCALE + +.. autodata:: DEGREE_SCALE + +.. autodata:: RADIAN_SCALE """ from gstools.tools.export import ( @@ -103,6 +112,15 @@ EARTH_RADIUS = 6371.0 """float: earth radius for WGS84 ellipsoid in km""" +KM_SCALE = 6371.0 +"""float: earth radius for WGS84 ellipsoid in km""" + +DEGREE_SCALE = 57.29577951308232 +"""float: radius for unit sphere in degree""" + +RADIAN_SCALE = 1.0 +"""float: radius for unit sphere""" + __all__ = [ "vtk_export", @@ -135,4 +153,7 @@ "generate_grid", "generate_st_grid", "EARTH_RADIUS", + "KM_SCALE", + "DEGREE_SCALE", + "RADIAN_SCALE", ] diff --git a/src/gstools/tools/geometric.py b/src/gstools/tools/geometric.py index afdcacaf..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 @@ -624,71 +625,94 @@ 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, time_scale=1.0 +): """Convert lat-lon geo coordinates to 3D position tuple. Parameters ---------- latlon : :class:`list` of :class:`numpy.ndarray` latitude and longitude given in degrees. + May includes an appended time axis if `time=True`. radius : :class:`float`, optional - Earth radius. Default: `1.0` + Sphere radius. Default: `1.0` dtype : data-type, optional The desired data-type for the array. If not given, then the type will be determined as the minimum type required to hold the objects in the sequence. Default: None + 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 ------- :class:`numpy.ndarray` the 3D position array """ - latlon = np.asarray(latlon, dtype=dtype).reshape((2, -1)) - lat, lon = np.deg2rad(latlon) - return np.array( - ( - radius * np.cos(lat) * np.cos(lon), - radius * np.cos(lat) * np.sin(lon), - radius * np.sin(lat) * np.ones_like(lon), - ), - dtype=dtype, + latlon = np.asarray(latlon, dtype=dtype).reshape((3 if time 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: + return np.array(pos_tuple + (latlon[2] / time_scale,), dtype=dtype) + return np.array(pos_tuple, dtype=dtype) -def pos2latlon(pos, radius=1.0, dtype=np.double): +def pos2latlon(pos, radius=1.0, dtype=np.double, time=False, time_scale=1.0): """Convert 3D position tuple from sphere to lat-lon geo coordinates. Parameters ---------- pos : :class:`list` of :class:`numpy.ndarray` The position tuple containing points on a unit-sphere. + May includes an appended time axis if `time=True`. radius : :class:`float`, optional - Earth radius. Default: `1.0` + Sphere radius. Default: `1.0` dtype : data-type, optional The desired data-type for the array. If not given, then the type will be determined as the minimum type required to hold the objects in the sequence. Default: None + 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 ------- :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]) - return np.rad2deg((lat, lon), dtype=dtype) + latlon = np.rad2deg((lat, lon), dtype=dtype) + if time: + return np.array( + (latlon[0], latlon[1], pos[3] * time_scale), dtype=dtype + ) + return latlon -def chordal_to_great_circle(dist): +def chordal_to_great_circle(dist, radius=1.0): """ Calculate great circle distance corresponding to given chordal distance. Parameters ---------- dist : array_like - Chordal distance of two points on the unit-sphere. + Chordal distance of two points on the sphere. + radius : :class:`float`, optional + Sphere radius. Default: `1.0` Returns ------- @@ -697,6 +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)) diff --git a/src/gstools/transform/field.py b/src/gstools/transform/field.py index 9ac33b6c..81824739 100644 --- a/src/gstools/transform/field.py +++ b/src/gstools/transform/field.py @@ -26,7 +26,7 @@ normal_to_arcsin normal_to_uquad """ -# pylint: disable=C0103, C0123, R0911 +# pylint: disable=C0103, C0123, R0911, R1735 import numpy as np from gstools.normalizer import ( diff --git a/src/gstools/variogram/binning.py b/src/gstools/variogram/binning.py index be490110..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 b57e0354..8c5785c6 100644 --- a/src/gstools/variogram/variogram.py +++ b/src/gstools/variogram/variogram.py @@ -14,6 +14,7 @@ from gstools import config from gstools.normalizer.tools import remove_trend_norm_mean +from gstools.tools import RADIAN_SCALE from gstools.tools.geometric import ( ang2dir, format_struct_pos_shape, @@ -92,6 +93,8 @@ def vario_estimate( normalizer=None, trend=None, fit_normalizer=False, + geo_scale=RADIAN_SCALE, + **std_bins, ): r""" Estimates the empirical variogram. @@ -222,6 +225,15 @@ 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` + **std_bins + Optional arguments that are forwarded to the :any:`standard_bins` routine + if no bins are given (bin_no, max_dist). Returns ------- @@ -250,18 +262,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 @@ -333,8 +345,13 @@ 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_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + 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 + bin_edges /= geo_scale # normalize field norm_field_out = remove_trend_norm_mean( *(pos, field, mean, normalizer, trend), @@ -371,7 +388,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( 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(