diff --git a/CHANGELOG.md b/CHANGELOG.md index 6e6d92ace..be0e69b14 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,35 @@ All notable changes to **GSTools** will be documented in this file. + +## [1.3.3] - Pure Pink - ? + +### Enhancements +See: [#197](https://github.com/GeoStat-Framework/GSTools/issues/197) +- `gstools.transform`: + - add keywords `field`, `store` and `process` keyword to all transformations to control storage and respect `normalizer` + - added `apply_function` transformation + - added `apply` as wrapper for all transformations + - added `transform` method to all `Field` (sub)classes as interface to `transform.apply` + - added checks for normal fields to work smoothly with recently added `normalizer` submodule +- `Field`: + - allow naming fields when generating and control storage with `store` keyword + - all subclasses now have the `post_process` keyword (apply mean, normalizer, trend) + - added subscription to access fields by name (`Field["field"]`) + - added `set_pos` method to set position tuple + - allow reusing present `pos` tuple + - added `pos`, `mesh_type`, `field_names`, `field_shape`, `all_fields` properties +- `CondSRF`: + - memory optimization by forwarding `pos` from underlying `krige` instance + - only recalculate kriging field if `pos` tuple changed (optimized ensemble generation) +- performance improvement by using `np.asarray` instead of `np.array` where possible +- updated examples to use new features +- added incomplete lower gamma function `inc_gamma_low` (for TPLGaussian spectral density) + +### Bugfixes +- `inc_gamma` was defined wrong for integer `s < 0` + + ## [1.3.2] - Pure Pink - 2021-07 ### Bugfixes @@ -280,7 +309,7 @@ All notable changes to **GSTools** will be documented in this file. First release of GSTools. -[Unreleased]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.2...HEAD +[1.3.3]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.2...HEAD [1.3.2]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.1...v1.3.2 [1.3.1]: https://github.com/GeoStat-Framework/gstools/compare/v1.3.0...v1.3.1 [1.3.0]: https://github.com/GeoStat-Framework/gstools/compare/v1.2.1...v1.3.0 diff --git a/docs/source/random.rst b/docs/source/random.rst index 4bce61b0e..c43314996 100644 --- a/docs/source/random.rst +++ b/docs/source/random.rst @@ -2,8 +2,6 @@ gstools.random ============== .. automodule:: gstools.random - :members: - :undoc-members: .. raw:: latex diff --git a/docs/source/tools.rst b/docs/source/tools.rst index f6c86fe4a..e3e91643a 100644 --- a/docs/source/tools.rst +++ b/docs/source/tools.rst @@ -2,8 +2,6 @@ gstools.tools ============= .. automodule:: gstools.tools - :members: - :undoc-members: .. raw:: latex diff --git a/docs/source/transform.rst b/docs/source/transform.rst index 1d41d590a..a456cb121 100644 --- a/docs/source/transform.rst +++ b/docs/source/transform.rst @@ -2,8 +2,6 @@ gstools.transform ================= .. automodule:: gstools.transform - :members: - :undoc-members: .. raw:: latex diff --git a/examples/01_random_field/01_srf_ensemble.py b/examples/01_random_field/01_srf_ensemble.py index 6174d4ce2..2563d9601 100644 --- a/examples/01_random_field/01_srf_ensemble.py +++ b/examples/01_random_field/01_srf_ensemble.py @@ -4,6 +4,8 @@ Creating an ensemble of random fields would also be a great idea. Let's reuse most of the previous code. + +We will set the position tuple `pos` before generation to reuse it afterwards. """ import numpy as np @@ -14,26 +16,26 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model) +srf.set_pos([x, y], "structured") ############################################################################### # This time, we did not provide a seed to :any:`SRF`, as the seeds will used # during the actual computation of the fields. We will create four ensemble -# members, for better visualisation and save them in a list and in a first +# members, for better visualisation, save them in to srf class and in a first # step, we will be using the loop counter as the seeds. - ens_no = 4 -field = [] for i in range(ens_no): - field.append(srf.structured([x, y], seed=i)) + srf(seed=i, store=f"field{i}") ############################################################################### -# Now let's have a look at the results: +# Now let's have a look at the results. We can access the fields by name or +# index: fig, ax = pt.subplots(2, 2, sharex=True, sharey=True) ax = ax.flatten() for i in range(ens_no): - ax[i].imshow(field[i].T, origin="lower") + ax[i].imshow(srf[i].T, origin="lower") pt.show() ############################################################################### @@ -44,9 +46,8 @@ # provides a seed generator :any:`MasterRNG`. The loop, in which the fields are # generated would then look like - from gstools.random import MasterRNG seed = MasterRNG(20170519) for i in range(ens_no): - field.append(srf.structured([x, y], seed=seed())) + srf(seed=seed(), store=f"better_field{i}") diff --git a/examples/06_conditioned_fields/00_condition_ensemble.py b/examples/06_conditioned_fields/00_condition_ensemble.py index a622fb463..551e867b2 100644 --- a/examples/06_conditioned_fields/00_condition_ensemble.py +++ b/examples/06_conditioned_fields/00_condition_ensemble.py @@ -26,15 +26,20 @@ model = gs.Gaussian(dim=1, var=0.5, len_scale=1.5) krige = gs.krige.Ordinary(model, cond_pos, cond_val) cond_srf = gs.CondSRF(krige) +cond_srf.set_pos(gridx) ############################################################################### +# To generate the ensemble we will use a seed-generator. +# We can specify individual names for each field by the keyword `store`: -fields = [] +seed = gs.random.MasterRNG(20170519) for i in range(100): - fields.append(cond_srf(gridx, seed=i)) + cond_srf(seed=seed(), store=f"f{i}") label = "Conditioned ensemble" if i == 0 else None - plt.plot(gridx, fields[i], color="k", alpha=0.1, label=label) -plt.plot(gridx, cond_srf.krige(gridx, only_mean=True), label="estimated mean") + plt.plot(gridx, cond_srf[f"f{i}"], color="k", alpha=0.1, label=label) + +fields = [cond_srf[f"f{i}"] for i in range(100)] +plt.plot(gridx, cond_srf.krige(only_mean=True), label="estimated mean") plt.plot(gridx, np.mean(fields, axis=0), linestyle=":", label="Ensemble mean") plt.plot(gridx, cond_srf.krige.field, linestyle="dashed", label="kriged field") plt.scatter(cond_pos, cond_val, color="k", zorder=10, label="Conditions") diff --git a/examples/06_conditioned_fields/01_2D_condition_ensemble.py b/examples/06_conditioned_fields/01_2D_condition_ensemble.py index 35d00bfd9..b77ee0f9a 100644 --- a/examples/06_conditioned_fields/01_2D_condition_ensemble.py +++ b/examples/06_conditioned_fields/01_2D_condition_ensemble.py @@ -20,14 +20,19 @@ model = gs.Gaussian(dim=2, var=0.5, len_scale=5, anis=0.5, angles=-0.5) krige = gs.Krige(model, cond_pos=cond_pos, cond_val=cond_val) cond_srf = gs.CondSRF(krige) +cond_srf.set_pos([x, y], "structured") ############################################################################### -# We create a list containing the generated conditioned fields. +# To generate the ensemble we will use a seed-generator. +# By specifying ``store=[f"fld{i}", False, False]``, only the conditioned field +# is stored with the specified name. The raw random field and the raw kriging +# field is not stored. This way, we can access each conditioned field by index +# ``cond_srf[i]``: +seed = gs.random.MasterRNG(20170519) ens_no = 4 -field = [] for i in range(ens_no): - field.append(cond_srf.structured([x, y], seed=i)) + cond_srf(seed=seed(), store=[f"fld{i}", False, False]) ############################################################################### # Now let's have a look at the pairwise differences between the generated @@ -35,19 +40,20 @@ fig, ax = plt.subplots(ens_no + 1, ens_no + 1, figsize=(8, 8)) # plotting kwargs for scatter and image -sc_kwargs = dict(c=cond_val, edgecolors="k", vmin=0, vmax=np.max(field)) -im_kwargs = dict(extent=2 * [0, 5], origin="lower", vmin=0, vmax=np.max(field)) +vmax = np.max(cond_srf.all_fields) +sc_kw = dict(c=cond_val, edgecolors="k", vmin=0, vmax=vmax) +im_kw = dict(extent=2 * [0, 5], origin="lower", vmin=0, vmax=vmax) for i in range(ens_no): # conditioned fields and conditions - ax[i + 1, 0].imshow(field[i].T, **im_kwargs) - ax[i + 1, 0].scatter(*cond_pos, **sc_kwargs) - ax[i + 1, 0].set_ylabel(f"Field {i+1}", fontsize=10) - ax[0, i + 1].imshow(field[i].T, **im_kwargs) - ax[0, i + 1].scatter(*cond_pos, **sc_kwargs) - ax[0, i + 1].set_title(f"Field {i+1}", fontsize=10) + ax[i + 1, 0].imshow(cond_srf[i].T, **im_kw) + ax[i + 1, 0].scatter(*cond_pos, **sc_kw) + ax[i + 1, 0].set_ylabel(f"Field {i}", fontsize=10) + ax[0, i + 1].imshow(cond_srf[i].T, **im_kw) + ax[0, i + 1].scatter(*cond_pos, **sc_kw) + ax[0, i + 1].set_title(f"Field {i}", fontsize=10) # absolute differences for j in range(ens_no): - ax[i + 1, j + 1].imshow(np.abs(field[i] - field[j]).T, **im_kwargs) + ax[i + 1, j + 1].imshow(np.abs(cond_srf[i] - cond_srf[j]).T, **im_kw) # beautify plots ax[0, 0].axis("off") @@ -56,3 +62,9 @@ a.set_xticks([]), a.set_yticks([]) fig.subplots_adjust(wspace=0, hspace=0) fig.show() + +############################################################################### +# To check if the generated fields are correct, we can have a look at their +# names: + +print(cond_srf.field_names) diff --git a/examples/07_transformations/00_log_normal.py b/examples/07_transformations/00_log_normal.py index a44f50fb2..3d2f45325 100755 --- a/examples/07_transformations/00_log_normal.py +++ b/examples/07_transformations/00_log_normal.py @@ -3,6 +3,8 @@ ----------------- Here we transform a field to a log-normal distribution: + +See :any:`transform.normal_to_lognormal` """ import gstools as gs @@ -11,5 +13,5 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519) srf.structured([x, y]) -gs.transform.normal_to_lognormal(srf) +srf.transform("normal_to_lognormal") # also "lognormal" works srf.plot() diff --git a/examples/07_transformations/01_binary.py b/examples/07_transformations/01_binary.py index 20b419191..a403abb7f 100755 --- a/examples/07_transformations/01_binary.py +++ b/examples/07_transformations/01_binary.py @@ -5,6 +5,8 @@ Here we transform a field to a binary field with only two values. The dividing value is the mean by default and the upper and lower values are derived to preserve the variance. + +See :any:`transform.binary` """ import gstools as gs @@ -13,5 +15,5 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519) srf.structured([x, y]) -gs.transform.binary(srf) +srf.transform("binary") srf.plot() diff --git a/examples/07_transformations/02_discrete.py b/examples/07_transformations/02_discrete.py index 325018984..a490fb5b5 100755 --- a/examples/07_transformations/02_discrete.py +++ b/examples/07_transformations/02_discrete.py @@ -6,32 +6,38 @@ If we do not give thresholds, the pairwise means of the given values are taken as thresholds. If thresholds are given, arbitrary values can be applied to the field. + +See :any:`transform.discrete` """ import numpy as np import gstools as gs -# structured field with a size of 100x100 and a grid-size of 0.5x0.5 +# Structured field with a size of 100x100 and a grid-size of 0.5x0.5 x = y = np.arange(200) * 0.5 model = gs.Gaussian(dim=2, var=1, len_scale=5) srf = gs.SRF(model, seed=20170519) - -# create 5 equidistanly spaced values, thresholds are the arithmetic means srf.structured([x, y]) -discrete_values = np.linspace(np.min(srf.field), np.max(srf.field), 5) -gs.transform.discrete(srf, discrete_values) -srf.plot() -# calculate thresholds for equal shares +############################################################################### +# Create 5 equidistanly spaced values, thresholds are the arithmetic means + +values1 = np.linspace(np.min(srf.field), np.max(srf.field), 5) +srf.transform("discrete", store="f1", values=values1) +srf.plot("f1") + +############################################################################### +# Calculate thresholds for equal shares # but apply different values to the separated classes -discrete_values2 = [0, -1, 2, -3, 4] -srf.structured([x, y]) -gs.transform.discrete(srf, discrete_values2, thresholds="equal") -srf.plot() -# user defined thresholds +values2 = [0, -1, 2, -3, 4] +srf.transform("discrete", store="f2", values=values2, thresholds="equal") +srf.plot("f2") + +############################################################################### +# Create user defined thresholds +# and apply different values to the separated classes + +values3 = [0, 1, 10] thresholds = [-1, 1] -# apply different values to the separated classes -discrete_values3 = [0, 1, 10] -srf.structured([x, y]) -gs.transform.discrete(srf, discrete_values3, thresholds=thresholds) -srf.plot() +srf.transform("discrete", store="f3", values=values3, thresholds=thresholds) +srf.plot("f3") diff --git a/examples/07_transformations/03_zinn_harvey.py b/examples/07_transformations/03_zinn_harvey.py index d3d28a1a7..01e65790c 100755 --- a/examples/07_transformations/03_zinn_harvey.py +++ b/examples/07_transformations/03_zinn_harvey.py @@ -6,6 +6,8 @@ `Zinn & Harvey (2003) `__. With this transformation, one could overcome the restriction that in ordinary Gaussian random fields the mean values are the ones being the most connected. + +See :any:`transform.zinnharvey` """ import gstools as gs @@ -14,5 +16,5 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519) srf.structured([x, y]) -gs.transform.zinnharvey(srf, conn="high") +srf.transform("zinnharvey", conn="high") srf.plot() diff --git a/examples/07_transformations/04_bimodal.py b/examples/07_transformations/04_bimodal.py index e0d4b2ded..8234a486a 100755 --- a/examples/07_transformations/04_bimodal.py +++ b/examples/07_transformations/04_bimodal.py @@ -1,6 +1,6 @@ """ -bimodal fields ---------------- +Bimodal fields +-------------- We provide two transformations to obtain bimodal distributions: @@ -8,6 +8,8 @@ * `uquad `__. Both transformations will preserve the mean and variance of the given field by default. + +See: :any:`transform.normal_to_arcsin` and :any:`transform.normal_to_uquad` """ import gstools as gs @@ -16,5 +18,5 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, seed=20170519) field = srf.structured([x, y]) -gs.transform.normal_to_arcsin(srf) +srf.transform("normal_to_arcsin") # also "arcsin" works srf.plot() diff --git a/examples/07_transformations/05_combinations.py b/examples/07_transformations/05_combinations.py index 2edb1a180..fb27c5e69 100755 --- a/examples/07_transformations/05_combinations.py +++ b/examples/07_transformations/05_combinations.py @@ -9,7 +9,13 @@ Then we apply the Zinn & Harvey transformation to connect the low values. Afterwards the field is transformed to a binary field and last but not least, we transform it to log-values. + +We can select the desired field by its name and we can define an output name +to store the field. + +If you don't specify `field` and `store` everything happens inplace. """ +# sphinx_gallery_thumbnail_number = 1 import gstools as gs # structured field with a size of 100x100 and a grid-size of 1x1 @@ -17,13 +23,18 @@ model = gs.Gaussian(dim=2, var=1, len_scale=10) srf = gs.SRF(model, mean=-9, seed=20170519) srf.structured([x, y]) -gs.transform.normal_force_moments(srf) -gs.transform.zinnharvey(srf, conn="low") -gs.transform.binary(srf) -gs.transform.normal_to_lognormal(srf) -srf.plot() +srf.transform("force_moments", field="field", store="f_forced") +srf.transform("zinnharvey", field="f_forced", store="f_zinnharvey", conn="low") +srf.transform("binary", field="f_zinnharvey", store="f_binary") +srf.transform("lognormal", field="f_binary", store="f_result") +srf.plot(field="f_result") ############################################################################### # The resulting field could be interpreted as a transmissivity field, where # the values of low permeability are the ones being the most connected # and only two kinds of soil exist. +# +# All stored fields can be accessed and plotted by name: + +print("Max binary value:", srf.f_binary.max()) +srf.plot(field="f_zinnharvey") diff --git a/examples/07_transformations/README.rst b/examples/07_transformations/README.rst index 2be183ee0..d93c99307 100644 --- a/examples/07_transformations/README.rst +++ b/examples/07_transformations/README.rst @@ -20,18 +20,31 @@ common transformations: normal_to_uniform normal_to_arcsin normal_to_uquad + apply_function All the transformations take a field class, that holds a generated field, -as input and will manipulate this field inplace. +as input and will manipulate this field inplace or store it with a given name. -Simply import the transform submodule and apply a transformation to the srf class: +Simply apply a transformation to a field class: .. code-block:: python - from gstools import transform as tf + import gstools as gs ... - tf.normal_to_lognormal(srf) + srf = gs.SRF(model) + srf(...) + gs.transform.normal_to_lognormal(srf) + +Or use the provided wrapper: + +.. code-block:: python + + import gstools as gs + ... + srf = gs.SRF(model) + srf(...) + srf.transform("lognormal") Examples -------- diff --git a/examples/08_geo_coordinates/01_dwd_krige.py b/examples/08_geo_coordinates/01_dwd_krige.py index bac4b31cb..2ddc55ebd 100755 --- a/examples/08_geo_coordinates/01_dwd_krige.py +++ b/examples/08_geo_coordinates/01_dwd_krige.py @@ -130,8 +130,9 @@ def north_south_drift(lat, lon): g_lat = np.arange(47, 56.1, 0.1) g_lon = np.arange(5, 16.1, 0.1) -field, k_var = uk((g_lat, g_lon), mesh_type="structured") -mean = uk((g_lat, g_lon), mesh_type="structured", only_mean=True) +uk.set_pos((g_lat, g_lon), mesh_type="structured") +uk(return_var=False, store="temp_field") +uk(only_mean=True, store="mean_field") ############################################################################### # And that's it. Now let's have a look at the generated field and the input @@ -140,8 +141,8 @@ def north_south_drift(lat, lon): levels = np.linspace(5, 23, 64) fig, ax = plt.subplots(1, 3, figsize=[10, 5], sharey=True) sca = ax[0].scatter(lon, lat, c=temp, vmin=5, vmax=23, cmap="coolwarm") -co1 = ax[1].contourf(g_lon, g_lat, field, levels, cmap="coolwarm") -co2 = ax[2].contourf(g_lon, g_lat, mean, levels, cmap="coolwarm") +co1 = ax[1].contourf(g_lon, g_lat, uk["temp_field"], levels, cmap="coolwarm") +co2 = ax[2].contourf(g_lon, g_lat, uk["mean_field"], levels, cmap="coolwarm") [ax[i].plot(border[:, 0], border[:, 1], color="k") for i in range(3)] [ax[i].set_xlim([5, 16]) for i in range(3)] @@ -160,8 +161,8 @@ def north_south_drift(lat, lon): # a look at a cross-section at a longitude of 10 degree: fig, ax = plt.subplots() -ax.plot(g_lat, field[:, 50], label="Interpolated temperature") -ax.plot(g_lat, mean[:, 50], label="North-South mean drift") +ax.plot(g_lat, uk["temp_field"][:, 50], label="Interpolated temperature") +ax.plot(g_lat, uk["mean_field"][:, 50], label="North-South mean drift") ax.set_xlabel("Lat in deg") ax.set_ylabel("T in [°C]") ax.set_title("North-South cross-section at 10°") diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index f8c86589c..76de7d197 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -265,7 +265,7 @@ def cor_spatial(self, pos): def vario_nugget(self, r): """Isotropic variogram of the model respecting the nugget at r=0.""" - r = np.array(np.abs(r), dtype=np.double) + r = np.asarray(np.abs(r), dtype=np.double) r_gz = np.logical_not(np.isclose(r, 0)) res = np.empty_like(r, dtype=np.double) res[r_gz] = self.variogram(r[r_gz]) @@ -274,7 +274,7 @@ def vario_nugget(self, r): def cov_nugget(self, r): """Isotropic covariance of the model respecting the nugget at r=0.""" - r = np.array(np.abs(r), dtype=np.double) + r = np.asarray(np.abs(r), dtype=np.double) r_gz = np.logical_not(np.isclose(r, 0)) res = np.empty_like(r, dtype=np.double) res[r_gz] = self.covariance(r[r_gz]) @@ -501,7 +501,7 @@ def spectral_density(self, k): k : :class:`float` Radius of the phase: :math:`k=\left\Vert\mathbf{k}\right\Vert` """ - k = np.array(np.abs(k), dtype=np.double) + k = np.asarray(np.abs(k), dtype=np.double) return self._sft.transform(self.correlation, k, ret_err=False) def spectral_rad_pdf(self, r): @@ -525,14 +525,14 @@ def _has_ppf(self): def isometrize(self, pos): """Make a position tuple ready for isotropic operations.""" - pos = np.array(pos, dtype=np.double).reshape((self.field_dim, -1)) + pos = np.asarray(pos, dtype=np.double).reshape((self.field_dim, -1)) if self.latlon: return latlon2pos(pos) 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.array(pos, dtype=np.double).reshape((self.dim, -1)) + pos = np.asarray(pos, dtype=np.double).reshape((self.dim, -1)) if self.latlon: return pos2latlon(pos) return np.dot( diff --git a/gstools/covmodel/fit.py b/gstools/covmodel/fit.py index c67c9839b..76dfb1375 100755 --- a/gstools/covmodel/fit.py +++ b/gstools/covmodel/fit.py @@ -321,8 +321,8 @@ def _pre_init_guess(model, init_guess, mean_x=1.0, mean_y=1.0): def _check_vario(model, x_data, y_data): # prepare variogram data - x_data = np.array(x_data).reshape(-1) - y_data = np.array(y_data).reshape(-1) + x_data = np.asarray(x_data).reshape(-1) + y_data = np.asarray(y_data).reshape(-1) # if multiple variograms are given, they will be interpreted # as directional variograms along the main rotated axes of the model is_dir_vario = False @@ -352,12 +352,10 @@ def _set_weights(model, weights, x_data, curve_fit_kwargs, is_dir_vario): weights = 1.0 / weights(x_data) elif isinstance(weights, str) and weights == "inv": weights = 1.0 + x_data - elif is_dir_vario: - if weights.size * model.dim == x_data.size: - weights = np.tile(weights, model.dim) - weights = 1.0 / np.array(weights).reshape(-1) else: - weights = 1.0 / np.array(weights).reshape(-1) + if is_dir_vario and weights.size * model.dim == x_data.size: + weights = np.tile(weights, model.dim) + weights = 1.0 / np.asarray(weights).reshape(-1) curve_fit_kwargs["sigma"] = weights curve_fit_kwargs["absolute_sigma"] = True diff --git a/gstools/covmodel/models.py b/gstools/covmodel/models.py index c05a92dcb..9c770ea59 100644 --- a/gstools/covmodel/models.py +++ b/gstools/covmodel/models.py @@ -73,14 +73,14 @@ def default_rescale(self): return np.sqrt(np.pi) / 2.0 def spectral_density(self, k): # noqa: D102 - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) return (self.len_rescaled / 2.0 / np.sqrt(np.pi)) ** self.dim * np.exp( -((k * self.len_rescaled / 2.0) ** 2) ) def spectral_rad_cdf(self, r): """Gaussian radial spectral cdf.""" - r = np.array(r, dtype=np.double) + r = np.asarray(r, dtype=np.double) if self.dim == 1: return sps.erf(r * self.len_rescaled / 2.0) if self.dim == 2: @@ -100,7 +100,7 @@ def spectral_rad_ppf(self, u): ----- Not defined for 3D. """ - u = np.array(u, dtype=np.double) + u = np.asarray(u, dtype=np.double) if self.dim == 1: return 2.0 / self.len_rescaled * sps.erfinv(u) if self.dim == 2: @@ -143,7 +143,7 @@ def cor(self, h): return np.exp(-h) def spectral_density(self, k): # noqa: D102 - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) return ( self.len_rescaled ** self.dim * sps.gamma((self.dim + 1) / 2.0) @@ -153,7 +153,7 @@ def spectral_density(self, k): # noqa: D102 def spectral_rad_cdf(self, r): """Exponential radial spectral cdf.""" - r = np.array(r, dtype=np.double) + r = np.asarray(r, dtype=np.double) if self.dim == 1: return np.arctan(r * self.len_rescaled) * 2.0 / np.pi if self.dim == 2: @@ -178,7 +178,7 @@ def spectral_rad_ppf(self, u): ----- Not defined for 3D. """ - u = np.array(u, dtype=np.double) + u = np.asarray(u, dtype=np.double) if self.dim == 1: return np.tan(np.pi / 2 * u) / self.len_rescaled if self.dim == 2: @@ -341,7 +341,7 @@ def default_opt_arg_bounds(self): def cor(self, h): """Matérn normalized correlation function.""" - h = np.array(np.abs(h), dtype=np.double) + h = np.asarray(np.abs(h), dtype=np.double) # for nu > 20 we just use the gaussian model if self.nu > 20.0: return np.exp(-((h / 2.0) ** 2)) @@ -360,7 +360,7 @@ def cor(self, h): return res def spectral_density(self, k): # noqa: D102 - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) # for nu > 20 we just use an approximation of the gaussian model if self.nu > 20.0: return ( @@ -570,7 +570,7 @@ class Circular(CovModel): def cor(self, h): """Circular normalized correlation function.""" - h = np.array(np.abs(h), dtype=np.double) + h = np.asarray(np.abs(h), dtype=np.double) res = np.zeros_like(h) # arccos is instable around h=1 h_l1 = h < 1.0 @@ -661,7 +661,7 @@ class HyperSpherical(CovModel): def cor(self, h): """Hyper-Spherical normalized correlation function.""" - h = np.array(h, dtype=np.double) + h = np.asarray(h, dtype=np.double) res = np.zeros_like(h) h_l1 = h < 1 nu = (self.dim - 1.0) / 2.0 @@ -670,7 +670,7 @@ def cor(self, h): return res def spectral_density(self, k): # noqa: D102 - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) res = np.empty_like(k) kl = k * self.len_rescaled kl_gz = np.logical_not(np.isclose(k, 0)) @@ -751,7 +751,7 @@ def default_opt_arg_bounds(self): def cor(self, h): """Super-Spherical normalized correlation function.""" - h = np.array(h, dtype=np.double) + h = np.asarray(h, dtype=np.double) res = np.zeros_like(h) h_l1 = h < 1 fac = 1.0 / sps.hyp2f1(0.5, -self.nu, 1.5, 1.0) @@ -846,7 +846,7 @@ def check_opt_arg(self): def cor(self, h): """J-Bessel correlation.""" - h = np.array(h, dtype=np.double) + h = np.asarray(h, dtype=np.double) h_gz = np.logical_not(np.isclose(h, 0)) hh = h[h_gz] res = np.ones_like(h) @@ -855,7 +855,7 @@ def cor(self, h): return res def spectral_density(self, k): # noqa: D102 - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) k_ll = k < 1.0 / self.len_rescaled kk = k[k_ll] res = np.zeros_like(k) diff --git a/gstools/covmodel/tools.py b/gstools/covmodel/tools.py index b3abf1ab2..56e65ec0e 100644 --- a/gstools/covmodel/tools.py +++ b/gstools/covmodel/tools.py @@ -71,12 +71,12 @@ def correlation(self, r): def correlation_from_cor(self, r): """Correlation function of the model.""" - r = np.array(np.abs(r), dtype=np.double) + r = np.asarray(np.abs(r), dtype=np.double) return self.cor(r / self.len_rescaled) def cor_from_correlation(self, h): """Correlation taking a non-dimensional range.""" - h = np.array(np.abs(h), dtype=np.double) + h = np.asarray(np.abs(h), dtype=np.double) return self.correlation(h * self.len_rescaled) abstract = True @@ -270,7 +270,7 @@ def check_arg_in_bounds(model, arg, val=None): raise ValueError("check bounds: unknown argument: {}".format(arg)) bnd = list(model.arg_bounds[arg]) val = getattr(model, arg) if val is None else val - val = np.array(val) + val = np.asarray(val) error_case = 0 if len(bnd) == 2: bnd.append("cc") # use closed intervals by default @@ -332,7 +332,7 @@ def spectral_rad_pdf(model, r): PDF values. """ - r = np.array(np.abs(r), dtype=np.double) + r = np.asarray(np.abs(r), dtype=np.double) if model.dim > 1: r_gz = np.logical_not(np.isclose(r, 0)) # to prevent numerical errors, we just calculate where r>0 diff --git a/gstools/field/base.py b/gstools/field/base.py index 7c83f370b..10c22126b 100755 --- a/gstools/field/base.py +++ b/gstools/field/base.py @@ -11,6 +11,8 @@ """ # pylint: disable=C0103, C0415 from functools import partial +from collections.abc import Iterable +from copy import copy import numpy as np from gstools.covmodel.base import CovModel from gstools.tools.geometric import format_struct_pos_dim, generate_grid @@ -18,8 +20,10 @@ generate_on_mesh, to_vtk_helper, fmt_mean_norm_trend, + _names, ) from gstools.normalizer.tools import apply_mean_norm_trend, _check_normalizer +from gstools.transform.field import apply __all__ = ["Field"] @@ -27,10 +31,23 @@ """:class:`list` of :class:`str`: valid field value types.""" +def _pos_equal(pos1, pos2): + if pos1 is None or pos2 is None: + return False + if len(pos1) != len(pos2): + return False + for p1, p2 in zip(pos1, pos2): + if len(p1) != len(p2): + return False + if not np.allclose(p1, p2): + return False + return True + + def _set_mean_trend(value, dim): if callable(value) or value is None: return value - value = np.array(value, dtype=np.double).ravel() + value = np.asarray(value, dtype=np.double).ravel() if value.size > 1 and value.size != dim: # vector mean raise ValueError(f"Mean/Trend: Wrong size ({value})") return value if value.size > 1 else value.item() @@ -60,6 +77,9 @@ class Field: Dimension of the field if no model is given. """ + default_field_names = ["field"] + """:class:`list`: Default field names.""" + def __init__( self, model=None, @@ -70,10 +90,10 @@ def __init__( dim=None, ): # initialize attributes - self.pos = None - self.mesh_type = None - self.field = None - # initialize private attributes + self._mesh_type = "unstructured" # default + self._pos = None + self._field_shape = None + self._field_names = [] self._model = None self._value_type = None self._mean = None @@ -87,14 +107,58 @@ def __init__( self.normalizer = normalizer self.trend = trend + def __len__(self): + return len(self.field_names) + + def __contains__(self, item): + return item in self.field_names + + def __getitem__(self, key): + if key in self.field_names: + return getattr(self, key) + if isinstance(key, int): + return self[self.field_names[key]] + if isinstance(key, slice): + return [self[f] for f in self.field_names[key]] + if isinstance(key, Iterable) and not isinstance(key, str): + return [self[f] for f in key] + raise KeyError(f"{self.name}: requested field '{key}' not present") + + def __delitem__(self, key): + names = [] + if key in self.field_names: + names = [key] + elif isinstance(key, int): + names = [self.field_names[key]] + elif isinstance(key, slice): + names = self.field_names[key] + elif isinstance(key, Iterable) and not isinstance(key, str): + for k in key: + k = self.field_names[k] if isinstance(key, int) else k + names.append(k) + else: + raise KeyError(f"{self.name}: requested field '{key}' not present") + for name in names: + if name not in self.field_names: + raise KeyError( + f"{self.name}: requested field '{name}' not present" + ) + delattr(self, name) + del self._field_names[self._field_names.index(name)] + def __call__( - self, pos, field=None, mesh_type="unstructured", post_process=True + self, + pos=None, + field=None, + mesh_type="unstructured", + post_process=True, + store=True, ): """Generate the field. Parameters ---------- - pos : :class:`list` + pos : :class:`list`, optional the position tuple, containing main direction and transversal directions field : :class:`numpy.ndarray` or :any:`None`, optional @@ -104,24 +168,33 @@ def __call__( post_process : :class:`bool`, optional Whether to apply mean, normalizer and trend to the field. Default: `True` + store : :class:`str` or :class:`bool`, optional + Whether to store field (True/False) with default name + or with specified name. + The default is :any:`True` for default name "field". Returns ------- field : :class:`numpy.ndarray` the field values. """ + name, save = self.get_store_config(store) pos, shape = self.pre_pos(pos, mesh_type) if field is None: field = np.zeros(shape, dtype=np.double) else: - field = np.array(field, dtype=np.double).reshape(shape) - return self.post_field(field, process=post_process) + field = np.asarray(field, dtype=np.double).reshape(shape) + return self.post_field(field, name, post_process, save) def structured(self, *args, **kwargs): """Generate a field on a structured mesh. See :any:`__call__` """ + if self.pos is None: + self.mesh_type = "structured" + if not (args or "pos" in kwargs) and self.mesh_type == "unstructured": + raise ValueError("Field.structured: can't reuse present 'pos'") call = partial(self.__call__, mesh_type="structured") return call(*args, **kwargs) @@ -130,6 +203,10 @@ def unstructured(self, *args, **kwargs): See :any:`__call__` """ + if self.pos is None: + self.mesh_type = "unstructured" + if not (args or "pos" in kwargs) and self.mesh_type != "unstructured": + raise ValueError("Field.unstructured: can't reuse present 'pos'") call = partial(self.__call__, mesh_type="unstructured") return call(*args, **kwargs) @@ -173,7 +250,7 @@ def mesh( """ return generate_on_mesh(self, mesh, points, direction, name, **kwargs) - def pre_pos(self, pos, mesh_type="unstructured"): + def pre_pos(self, pos=None, mesh_type="unstructured", info=False): """ Preprocessing positions and mesh_type. @@ -185,34 +262,38 @@ def pre_pos(self, pos, mesh_type="unstructured"): mesh_type : :class:`str`, optional 'structured' / 'unstructured' Default: `"unstructured"` + info : :class:`bool`, optional + Whether to return information Returns ------- iso_pos : (d, n), :class:`numpy.ndarray` - the isometrized position tuple + Isometrized position tuple. shape : :class:`tuple` Shape of the resulting field. + info : :class:`dict`, optional + Information about settings. + + Warnings + -------- + When setting a new position tuple that differs from the present one, + all stored fields will be deleted. """ - # save mesh-type - self.mesh_type = mesh_type - # save pos tuple - if mesh_type != "unstructured": - pos, shape = format_struct_pos_dim(pos, self.dim) - self.pos = pos - pos = generate_grid(pos) + info_ret = {"deleted": False} + if pos is None: + if self.pos is None: + raise ValueError("Field: no position tuple 'pos' present") else: - pos = np.array(pos, dtype=np.double).reshape(self.dim, -1) - self.pos = pos - shape = np.shape(pos[0]) - # prepend dimension if we have a vector field - if self.value_type == "vector": - shape = (self.dim,) + shape - if self.latlon: - raise ValueError("Field: Vector fields not allowed for latlon") - # return isometrized pos tuple and resulting field shape + info_ret = self.set_pos(pos, mesh_type, info=True) + if self.mesh_type != "unstructured": + pos = generate_grid(self.pos) + else: + pos = self.pos + # return isometrized pos tuple, field shape and possible info + info_ret = (info_ret,) if self.model is None: - return pos, shape - return self.model.isometrize(pos), shape + return (pos, self.field_shape) + info * info_ret + return (self.model.isometrize(pos), self.field_shape) + info * info_ret def post_field(self, field, name="field", process=True, save=True): """ @@ -237,9 +318,10 @@ def post_field(self, field, name="field", process=True, save=True): field : :class:`numpy.ndarray` Processed field values. """ + if self.field_shape is None: + raise ValueError("post_field: no 'field_shape' present.") + field = np.asarray(field, dtype=np.double).reshape(self.field_shape) if process: - if self.pos is None: - raise ValueError("post_field: no 'pos' tuple set for field.") field = apply_mean_norm_trend( pos=self.pos, field=field, @@ -252,9 +334,59 @@ def post_field(self, field, name="field", process=True, save=True): stacked=False, ) if save: + name = str(name) + if not name.isidentifier() or ( + name not in self.field_names and name in dir(self) + ): + raise ValueError( + f"Field: given field name '{name}' is not valid" + ) + # allow resetting present fields + if name not in self._field_names: + self._field_names.append(name) setattr(self, name, field) return field + def delete_fields(self, select=None): + """Delete selected fields.""" + del self[self.field_names if select is None else select] + + def transform( + self, method, field="field", store=True, process=False, **kwargs + ): + """ + Apply field transformation. + + Parameters + ---------- + method : :class:`str` + Method to use. + See :any:`gstools.transform` for available transformations. + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + **kwargs + Keyword arguments forwarded to selected method. + + Raises + ------ + ValueError + When method is unknown. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + return apply( + self, method, field=field, store=store, process=process, **kwargs + ) + def to_pyvista( self, field_select="field", fieldname="field" ): # pragma: no cover @@ -343,6 +475,139 @@ def plot( return r + def set_pos(self, pos, mesh_type="unstructured", info=False): + """ + Set positions and mesh_type. + + Parameters + ---------- + pos : :any:`iterable` + the position tuple, containing main direction and transversal + directions + mesh_type : :class:`str`, optional + 'structured' / 'unstructured' + Default: `"unstructured"` + info : :class:`bool`, optional + Whether to return information + + Returns + ------- + info : :class:`dict`, optional + Information about settings. + + Warnings + -------- + When setting a new position tuple that differs from the present one, + all stored fields will be deleted. + """ + info_ret = {"deleted": False} + old_type = copy(self.mesh_type) + old_pos = copy(self.pos) + # save pos and mesh-type + self.mesh_type = mesh_type + self.pos = pos + # remove present fields if new pos is different from current + if old_type != self.mesh_type or not _pos_equal(old_pos, self.pos): + self.delete_fields() + info_ret["deleted"] = True + del old_pos + return info_ret if info else None + + def get_store_config(self, store, default=None, fld_cnt=None): + """ + Get storage configuration from given selection. + + Parameters + ---------- + store : :class:`str` or :class:`bool` or :class:`list`, optional + Whether to store fields (True/False) with default names + or with specified names. + The default is :any:`True` for default names. + default : :class:`str` or :class:`list`, optional + Default field names. The default is "field". + fld_cnt : :any:`None` or :class:`int`, optional + Number of fields when using lists. The default is None. + + Returns + ------- + name : :class:`str` or :class:`list` + Name(s) of field. + save : :class:`bool` or :class:`list` + Whether to save field(s). + """ + if default is None: + if fld_cnt is None: + default = self.default_field_names[0] + else: + default = self.default_field_names + # single field + if fld_cnt is None: + save = isinstance(store, str) or bool(store) + name = store if isinstance(store, str) else default + return name, save + # multiple fields + default = _names(default, fld_cnt) + save = [True] * fld_cnt + if isinstance(store, str): + store = [store] + if isinstance(store, Iterable): + store = list(store)[:fld_cnt] + store += [True] * (fld_cnt - len(store)) + name = [None] * fld_cnt + for i, val in enumerate(store): + save[i] = isinstance(val, str) or bool(val) + name[i] = val if isinstance(val, str) else default[i] + else: + save = [bool(store)] * fld_cnt + name = copy(default) + return name, save + + @property + def pos(self): + """:class:`tuple`: The position tuple of the field.""" + return self._pos + + @pos.setter + def pos(self, pos): + if self.mesh_type == "unstructured": + self._pos = np.asarray(pos, dtype=np.double).reshape(self.dim, -1) + self._field_shape = np.shape(self._pos[0]) + else: + self._pos, self._field_shape = format_struct_pos_dim(pos, self.dim) + # prepend dimension if we have a vector field + if self.value_type == "vector": + self._field_shape = (self.dim,) + self._field_shape + if self.latlon: + raise ValueError("Field: Vector fields not allowed for latlon") + + @property + def all_fields(self): + """:class:`list`: All fields as stacked list.""" + return self[self.field_names] + + @property + def field_names(self): + """:class:`list`: Names of present fields.""" + return self._field_names + + @field_names.deleter + def field_names(self): + self.delete_fields() + + @property + def field_shape(self): + """:class:`tuple`: The shape of the field.""" + return self._field_shape + + @property + def mesh_type(self): + """:class:`str`: The mesh type of the field.""" + return self._mesh_type + + @mesh_type.setter + def mesh_type(self, mesh_type): + self._mesh_type = str(mesh_type) + @property def model(self): """:any:`CovModel`: The covariance model of the field.""" diff --git a/gstools/field/cond_srf.py b/gstools/field/cond_srf.py index 609045c83..5d424c6ec 100644 --- a/gstools/field/cond_srf.py +++ b/gstools/field/cond_srf.py @@ -43,34 +43,56 @@ class CondSRF(Field): Have a look at the provided generators for further information. """ + default_field_names = ["field", "raw_field", "raw_krige"] + """:class:`list`: Default field names.""" + def __init__(self, krige, generator="RandMeth", **generator_kwargs): if not isinstance(krige, Krige): raise ValueError("CondSRF: krige should be an instance of Krige.") self._krige = krige # initialize attributes - self.pos = None - self.mesh_type = None - self.field = None - self.raw_field = None + self._field_names = [] # initialize private attributes self._generator = None # initialize attributes self.set_generator(generator, **generator_kwargs) - def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs): + def __call__( + self, + pos=None, + seed=np.nan, + mesh_type="unstructured", + post_process=True, + store=True, + krige_store=True, + **kwargs, + ): """Generate the conditioned spatial random field. The field is saved as `self.field` and is also returned. Parameters ---------- - pos : :class:`list` + pos : :class:`list`, optional the position tuple, containing main direction and transversal directions seed : :class:`int`, optional seed for RNG for reseting. Default: keep seed from generator mesh_type : :class:`str` 'structured' / 'unstructured' + post_process : :class:`bool`, optional + Whether to apply mean, normalizer and trend to the field. + Default: `True` + store : :class:`str` or :class:`bool` or :class:`list`, optional + Whether to store fields (True/False) with default names + or with specified names. + The default is :any:`True` for default names + ["field", "raw_field", "raw_krige"]. + krige_store : :class:`str` or :class:`bool` or :class:`list`, optional + Whether to store kriging fields (True/False) with default name + or with specified names. + The default is :any:`True` for default names + ["field", "krige_var"]. **kwargs keyword arguments that are forwarded to the kriging routine in use. @@ -79,23 +101,50 @@ def __call__(self, pos, seed=np.nan, mesh_type="unstructured", **kwargs): field : :class:`numpy.ndarray` the conditioned SRF """ + name, save = self.get_store_config(store=store, fld_cnt=3) + krige_name, krige_save = self.krige.get_store_config( + store=krige_store, fld_cnt=2 + ) kwargs["mesh_type"] = mesh_type kwargs["only_mean"] = False # overwrite if given kwargs["return_var"] = True # overwrite if given kwargs["post_process"] = False # overwrite if given + kwargs["store"] = [False, krige_name[1] if krige_save[1] else False] # update the model/seed in the generator if any changes were made self.generator.update(self.model, seed) # get isometrized positions and the resulting field-shape - iso_pos, shape = self.pre_pos(pos, mesh_type) + iso_pos, shape, info = self.pre_pos(pos, mesh_type, info=True) # generate the field - self.raw_field = np.reshape( - self.generator(iso_pos, add_nugget=False), shape - ) - field, krige_var = self.krige(pos, **kwargs) + rawfield = np.reshape(self.generator(iso_pos, add_nugget=False), shape) + # call krige on already set pos (reuse already calculated fields) + if ( + not info["deleted"] + and name[2] in self.field_names + and krige_name[1] in self.krige.field_names + ): + reuse = True + rawkrige, krige_var = self[name[2]], self.krige[krige_name[1]] + else: + reuse = False + rawkrige, krige_var = self.krige(**kwargs) var_scale, nugget = self.get_scaling(krige_var, shape) - # need to use a copy to not alter "field" by reference - self.krige.post_field(self.krige.field.copy()) - return self.post_field(field + var_scale * self.raw_field + nugget) + # store krige field (need a copy to not alter field by reference) + if not reuse or krige_name[0] not in self.krige.field_names: + self.krige.post_field( + rawkrige.copy(), krige_name[0], post_process, krige_save[0] + ) + # store raw krige field + if not reuse: + self.post_field(rawkrige, name[2], False, save[2]) + # store raw random field + self.post_field(rawfield, name[1], False, save[1]) + # store cond random field + return self.post_field( + field=rawkrige + var_scale * rawfield + nugget, + name=name[0], + process=post_process, + save=save[0], + ) def get_scaling(self, krige_var, shape): """ @@ -143,6 +192,59 @@ def set_generator(self, generator, **generator_kwargs): else: raise ValueError(f"gstools.CondSRF: Unknown generator {generator}") + def set_pos(self, pos, mesh_type="unstructured", info=False): + """ + Set positions and mesh_type. + + Parameters + ---------- + pos : :any:`iterable` + the position tuple, containing main direction and transversal + directions + mesh_type : :class:`str`, optional + 'structured' / 'unstructured' + Default: `"unstructured"` + info : :class:`bool`, optional + Whether to return information + + Returns + ------- + info : :class:`dict`, optional + Information about settings. + + Warnings + -------- + When setting a new position tuple that differs from the present one, + all stored fields will be deleted. + """ + info_ret = super().set_pos(pos, mesh_type, info=True) + if info_ret["deleted"]: + self.krige.delete_fields() + return info_ret if info else None + + @property + def pos(self): + """:class:`tuple`: The position tuple of the field.""" + return self.krige.pos + + @pos.setter + def pos(self, pos): + self.krige.pos = pos + + @property + def field_shape(self): + """:class:`tuple`: The shape of the field.""" + return self.krige.field_shape + + @property + def mesh_type(self): + """:class:`str`: The mesh type of the field.""" + return self.krige.mesh_type + + @mesh_type.setter + def mesh_type(self, mesh_type): + self.krige.mesh_type = mesh_type + @property def krige(self): """:any:`Krige`: The underlying kriging class.""" @@ -200,6 +302,6 @@ def value_type(self, value_type): def __repr__(self): """Return String representation.""" - return "CondSRF(krige={0}, generator={1})".format( - self.krige, self.generator.name + return "{0}(krige={1}, generator={2})".format( + self.name, self.krige, self.generator.name ) diff --git a/gstools/field/generator.py b/gstools/field/generator.py index 9b3eeca67..9c2761bf6 100644 --- a/gstools/field/generator.py +++ b/gstools/field/generator.py @@ -125,7 +125,7 @@ def __call__(self, pos, add_nugget=True): :class:`numpy.ndarray` the random modes """ - pos = np.array(pos, dtype=np.double) + pos = np.asarray(pos, dtype=np.double) summed_modes = summate(self._cov_sample, self._z_1, self._z_2, pos) nugget = self.get_nugget(summed_modes.shape) if add_nugget else 0.0 return np.sqrt(self.model.var / self._mode_no) * summed_modes + nugget @@ -415,7 +415,7 @@ def __call__(self, pos): :class:`numpy.ndarray` the random modes """ - pos = np.array(pos, dtype=np.double) + pos = np.asarray(pos, dtype=np.double) summed_modes = summate_incompr( self._cov_sample, self._z_1, self._z_2, pos ) diff --git a/gstools/field/plot.py b/gstools/field/plot.py index ca3a422b4..c0b589b21 100644 --- a/gstools/field/plot.py +++ b/gstools/field/plot.py @@ -51,12 +51,10 @@ def plot_field( **kwargs Forwarded to the plotting routine. """ - plt_fld = getattr(fld, field) - assert not (fld.pos is None or plt_fld is None) if fld.dim == 1: - return plot_1d(fld.pos, plt_fld, fig, ax, **kwargs) + return plot_1d(fld.pos, fld[field], fig, ax, **kwargs) return plot_nd( - fld.pos, plt_fld, fld.mesh_type, fig, ax, fld.latlon, **kwargs + fld.pos, fld[field], fld.mesh_type, fig, ax, fld.latlon, **kwargs ) @@ -302,9 +300,7 @@ def plot_vec_field(fld, field="field", fig=None, ax=None): # pragma: no cover "Only structured vector fields are supported " "for plotting. Please create one on a structured grid." ) - plt_fld = getattr(fld, field) - assert not (fld.pos is None or plt_fld is None) - + plt_fld = fld[field] norm = np.sqrt(plt_fld[0, :].T ** 2 + plt_fld[1, :].T ** 2) fig, ax = get_fig_ax(fig, ax) diff --git a/gstools/field/srf.py b/gstools/field/srf.py index de0951241..86d4ef6db 100644 --- a/gstools/field/srf.py +++ b/gstools/field/srf.py @@ -100,7 +100,13 @@ def __init__( self.set_generator(generator, **generator_kwargs) def __call__( - self, pos, seed=np.nan, point_volumes=0.0, mesh_type="unstructured" + self, + pos=None, + seed=np.nan, + point_volumes=0.0, + mesh_type="unstructured", + post_process=True, + store=True, ): """Generate the spatial random field. @@ -108,7 +114,7 @@ def __call__( Parameters ---------- - pos : :class:`list` + pos : :class:`list`, optional the position tuple, containing main direction and transversal directions seed : :class:`int`, optional @@ -121,12 +127,20 @@ def __call__( is changed. Default: ``0`` mesh_type : :class:`str` 'structured' / 'unstructured' + post_process : :class:`bool`, optional + Whether to apply mean, normalizer and trend to the field. + Default: `True` + store : :class:`str` or :class:`bool`, optional + Whether to store field (True/False) with default name + or with specified name. + The default is :any:`True` for default name "field". Returns ------- field : :class:`numpy.ndarray` the SRF """ + name, save = self.get_store_config(store) # update the model/seed in the generator if any changes were made self.generator.update(self.model, seed) # get isometrized positions and the resulting field-shape @@ -139,7 +153,7 @@ def __call__( if np.size(scaled_var) > 1: scaled_var = np.reshape(scaled_var, shape) field *= np.sqrt(scaled_var / self.model.sill) - return self.post_field(field) + return self.post_field(field, name, post_process, save) def upscaling_func(self, *args, **kwargs): """Upscaling method applied to the field variance.""" diff --git a/gstools/field/tools.py b/gstools/field/tools.py index fba20f94f..f46a0d305 100644 --- a/gstools/field/tools.py +++ b/gstools/field/tools.py @@ -70,11 +70,8 @@ def to_vtk_helper( fieldname : :class:`str`, optional Name of the field in the VTK file. Default: "field" """ + field = f_cls[field_select] if field_select in f_cls.field_names else None if f_cls.value_type == "vector": - if hasattr(f_cls, field_select): - field = getattr(f_cls, field_select) - else: - field = None if not (f_cls.pos is None or field is None or f_cls.mesh_type is None): suf = ["_X", "_Y", "_Z"] fields = {} @@ -85,10 +82,6 @@ def to_vtk_helper( return vtk_export(filename, f_cls.pos, fields, f_cls.mesh_type) raise ValueError(f"Field.to_vtk: '{field_select}' not available.") if f_cls.value_type == "scalar": - if hasattr(f_cls, field_select): - field = getattr(f_cls, field_select) - else: - field = None if not (f_cls.pos is None or field is None or f_cls.mesh_type is None): if filename is None: return to_vtk(f_cls.pos, {fieldname: field}, f_cls.mesh_type) diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 80cc21b49..0c87d0f1a 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -125,6 +125,9 @@ class Krige(Field): Springer, Berlin, Heidelberg (2003) """ + default_field_names = ["field", "krige_var", "mean_field"] + """:class:`list`: Default field names.""" + def __init__( self, model, @@ -144,8 +147,6 @@ def __init__( fit_variogram=False, ): super().__init__(model, mean=mean, normalizer=normalizer, trend=trend) - self.mean_field = None - self.krige_var = None self._unbiased = bool(unbiased) self._exact = bool(exact) self._pseudo_inv = bool(pseudo_inv) @@ -172,13 +173,14 @@ def __init__( def __call__( self, - pos, + pos=None, mesh_type="unstructured", ext_drift=None, chunk_size=None, only_mean=False, return_var=True, post_process=True, + store=True, ): """ Generate the kriging field. @@ -188,7 +190,7 @@ def __call__( Parameters ---------- - pos : :class:`list` + pos : :class:`list`, optional the position tuple, containing main direction and transversal directions (x, [y, z]) mesh_type : :class:`str`, optional @@ -208,6 +210,11 @@ def __call__( post_process : :class:`bool`, optional Whether to apply mean, normalizer and trend to the field. Default: `True` + store : :class:`str` or :class:`bool` or :class:`list`, optional + Whether to store kriging fields (True/False) with default name + or with specified names. + The default is :any:`True` for default names + ["field", "krige_var"] or "mean_field" if `only_mean=True`. Returns ------- @@ -218,6 +225,10 @@ def __call__( (if return_var is True and only_mean is False) """ return_var &= not only_mean # don't return variance when calc. mean + fld_cnt = 2 if return_var else 1 + default = self.default_field_names[2] if only_mean else None + name, save = self.get_store_config(store, default, fld_cnt) + iso_pos, shape = self.pre_pos(pos, mesh_type) pnt_cnt = len(iso_pos[0]) @@ -248,18 +259,14 @@ def __call__( self._summate(field, krige_var, c_slice, k_vec, return_var) # reshape field if we got a structured mesh field = np.reshape(field, shape) - if only_mean: # care about 'kriging the mean' - return self.post_field(field, "mean_field", process=post_process) # save field to class - field = self.post_field(field, "field", process=post_process) + field = self.post_field(field, name[0], post_process, save[0]) if return_var: # care about the estimated error variance krige_var = np.reshape( np.maximum(self.model.sill - krige_var, 0), shape ) - krige_var = self.post_field(krige_var, "krige_var", process=False) + krige_var = self.post_field(krige_var, name[1], False, save[1]) return field, krige_var - # if we only calculate the field, overwrite the error variance - self.krige_var = None return field def _summate(self, field, krige_var, c_slice, k_vec, return_var): @@ -364,7 +371,9 @@ def _pre_ext_drift(self, pnt_cnt, ext_drift=None, set_cond=False): the drift values at the given positions """ if ext_drift is not None: - ext_drift = np.array(ext_drift, dtype=np.double, ndmin=2) + ext_drift = np.array( + ext_drift, dtype=np.double, ndmin=2, copy=False + ) if ext_drift.size == 0: # treat empty array as no ext_drift return np.array([]) if set_cond: @@ -377,7 +386,7 @@ def _pre_ext_drift(self, pnt_cnt, ext_drift=None, set_cond=False): raise ValueError("Krige: wrong number of external drifts.") if np.prod(ext_shape) != np.prod(shape): raise ValueError("Krige: wrong number of ext. drift values.") - return np.array(ext_drift, dtype=np.double).reshape(shape) + return np.asarray(ext_drift, dtype=np.double).reshape(shape) if not set_cond and self._cond_ext_drift.size > 0: raise ValueError("Krige: wrong number of ext. drift values.") return np.array([]) @@ -603,7 +612,7 @@ def cond_err(self, value): "krige.cond_err: measurement errors can't be given, " "when interpolator should be exact." ) - value = np.array(value, dtype=np.double).reshape(-1) + value = np.asarray(value, dtype=np.double).reshape(-1) if value.size == 1: self._cond_err = value.item() else: @@ -689,11 +698,6 @@ def ext_drift_no(self): """:class:`int`: Number of external drift values per point.""" return self.cond_ext_drift.shape[0] - @property - def name(self): - """:class:`str`: The name of the kriging class.""" - return self.__class__.__name__ - def __repr__(self): """Return String representation.""" return "{0}(model={1}, cond_no={2}{3})".format( diff --git a/gstools/krige/tools.py b/gstools/krige/tools.py index 5d220f034..d84c44855 100644 --- a/gstools/krige/tools.py +++ b/gstools/krige/tools.py @@ -44,8 +44,8 @@ def set_condition(cond_pos, cond_val, dim): the error checked cond_val """ # convert the input for right shapes and dimension checks - cond_val = np.array(cond_val, dtype=np.double).reshape(-1) - cond_pos = np.array(cond_pos, dtype=np.double).reshape(dim, -1) + cond_val = np.asarray(cond_val, dtype=np.double).reshape(-1) + cond_pos = np.asarray(cond_pos, dtype=np.double).reshape(dim, -1) if len(cond_pos[0]) != len(cond_val): raise ValueError( "Please check your 'cond_pos' and 'cond_val' parameters. " diff --git a/gstools/normalizer/base.py b/gstools/normalizer/base.py index da6e0a9cb..a9fcf9b0c 100644 --- a/gstools/normalizer/base.py +++ b/gstools/normalizer/base.py @@ -67,10 +67,10 @@ def _kernel_loglikelihood(self, data): return res + np.sum(np.log(np.maximum(1e-16, self._derivative(data)))) def _check_input(self, data, data_range=None, return_output_template=True): - is_data = np.array(np.logical_not(np.isnan(data))) + is_data = np.logical_not(np.isnan(data)) if return_output_template: out = np.full_like(data, np.nan, dtype=np.double) - data = np.array(data, dtype=np.double)[is_data] + data = np.asarray(data, dtype=np.double)[is_data] if data_range is not None and np.min(np.abs(data_range)) < np.inf: dat_in = np.logical_and(data > data_range[0], data < data_range[1]) if not np.all(dat_in): diff --git a/gstools/normalizer/tools.py b/gstools/normalizer/tools.py index f71b758f0..0ec10f773 100644 --- a/gstools/normalizer/tools.py +++ b/gstools/normalizer/tools.py @@ -90,7 +90,7 @@ def apply_mean_norm_trend( pos, shape, dim = format_unstruct_pos_shape( pos, field.shape, check_stacked_shape=stacked ) - field = np.array(field, dtype=np.double).reshape(shape) + field = np.asarray(field, dtype=np.double).reshape(shape) else: dim = len(pos) if not stacked: @@ -169,7 +169,7 @@ def remove_trend_norm_mean( pos, shape, dim = format_unstruct_pos_shape( pos, field.shape, check_stacked_shape=stacked ) - field = np.array(field, dtype=np.double).reshape(shape) + field = np.asarray(field, dtype=np.double).reshape(shape) else: dim = len(pos) if not stacked: diff --git a/gstools/random/__init__.py b/gstools/random/__init__.py index ec78b43c0..6c02da31f 100644 --- a/gstools/random/__init__.py +++ b/gstools/random/__init__.py @@ -8,18 +8,24 @@ ^^^^^^^^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + RNG Seed Generator ^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + MasterRNG Distribution factory ^^^^^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + dist_gen ---- diff --git a/gstools/tools/__init__.py b/gstools/tools/__init__.py index e03529c3b..15e3c2c3f 100644 --- a/gstools/tools/__init__.py +++ b/gstools/tools/__init__.py @@ -8,6 +8,8 @@ ^^^^^^ .. autosummary:: + :toctree: generated + vtk_export vtk_export_structured vtk_export_unstructured @@ -19,8 +21,11 @@ ^^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + confidence_scaling inc_gamma + inc_gamma_low exp_int inc_beta tplstable_cor @@ -31,6 +36,8 @@ ^^^^^^^^^ .. autosummary:: + :toctree: generated + rotated_main_axes set_angles set_anis @@ -68,6 +75,7 @@ from gstools.tools.special import ( confidence_scaling, inc_gamma, + inc_gamma_low, exp_int, inc_beta, tplstable_cor, @@ -107,6 +115,7 @@ "to_vtk_unstructured", "confidence_scaling", "inc_gamma", + "inc_gamma_low", "exp_int", "inc_beta", "tplstable_cor", diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index 204df5aea..843114479 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -79,7 +79,7 @@ def set_angles(dim, angles): ----- If too few angles are given, they are filled up with `0`. """ - out_angles = np.array(angles, dtype=np.double) + out_angles = np.asarray(angles, dtype=np.double) out_angles = np.atleast_1d(out_angles)[: no_of_angles(dim)] # fill up the rotation angle array with zeros out_angles = np.pad( @@ -110,7 +110,7 @@ def set_anis(dim, anis): ----- If too few anisotropy ratios are given, they are filled up with `1`. """ - out_anis = np.array(anis, dtype=np.double) + out_anis = np.asarray(anis, dtype=np.double) out_anis = np.atleast_1d(out_anis)[: dim - 1] if len(out_anis) < dim - 1: # fill up the anisotropies with ones, such that len()==dim-1 @@ -351,9 +351,9 @@ def generate_grid(pos): :class:`numpy.ndarray` Unstructured position tuple. """ - return np.array(np.meshgrid(*pos, indexing="ij"), dtype=np.double).reshape( - (len(pos), -1) - ) + return np.asarray( + np.meshgrid(*pos, indexing="ij"), dtype=np.double + ).reshape((len(pos), -1)) def generate_st_grid(pos, time, mesh_type="unstructured"): @@ -379,14 +379,14 @@ def generate_st_grid(pos, time, mesh_type="unstructured"): ----- Time dimension will be the last one. """ - time = np.array(time, dtype=np.double).reshape(-1) + time = np.asarray(time, dtype=np.double).reshape(-1) if mesh_type != "unstructured": pos = generate_grid(pos) else: - pos = np.array(pos, dtype=np.double, ndmin=2) + pos = np.array(pos, dtype=np.double, ndmin=2, copy=False) out = [np.repeat(p.reshape(-1), np.size(time)) for p in pos] out.append(np.tile(time, np.size(pos[0]))) - return np.array(out, dtype=np.double) + return np.asarray(out, dtype=np.double) # conversion ################################################################## @@ -416,12 +416,12 @@ def format_struct_pos_dim(pos, dim): Shape of the resulting field. """ if dim == 1: - pos = (np.array(pos, dtype=np.double).reshape(-1),) + pos = (np.asarray(pos, dtype=np.double).reshape(-1),) elif len(pos) != dim: raise ValueError("Formatting: position tuple doesn't match dimension.") else: - pos = tuple(np.array(p_i, dtype=np.double).reshape(-1) for p_i in pos) - shape = tuple(len(p_i) for p_i in pos) + pos = tuple(np.asarray(p, dtype=np.double).reshape(-1) for p in pos) + shape = tuple(len(p) for p in pos) return pos, shape @@ -479,7 +479,7 @@ def format_struct_pos_shape(pos, shape, check_stacked_shape=False): else: wrong_shape = True else: - struct_size = np.prod([p_i.size for p_i in check_pos]) + struct_size = np.prod([p.size for p in check_pos]) # case: 1D unstacked if check_pos.size == shape_size: dim = 1 @@ -552,7 +552,7 @@ def format_unstruct_pos_shape(pos, shape, check_stacked_shape=False): # now we try to be smart pre_len = len(np.atleast_1d(pos)) # care about 1D: pos can be given as 1D array here -> convert to 2D array - pos = np.array(pos, dtype=np.double, ndmin=2) + pos = np.array(pos, dtype=np.double, ndmin=2, copy=False) post_len = len(pos) # first array dimension should be spatial dimension (1D is special case) dim = post_len if pre_len == post_len else 1 @@ -606,7 +606,7 @@ def ang2dir(angles, dtype=np.double, dim=None): the array of direction vectors """ pre_dim = np.asanyarray(angles).ndim - angles = np.array(angles, ndmin=2, dtype=dtype) + angles = np.array(angles, ndmin=2, dtype=dtype, copy=False) if len(angles.shape) > 2: raise ValueError(f"Can't interpret angles array {angles}") dim = angles.shape[1] + 1 if dim is None else dim @@ -644,7 +644,7 @@ def latlon2pos(latlon, radius=1.0, dtype=np.double): :class:`numpy.ndarray` the 3D position array """ - latlon = np.array(latlon, dtype=dtype).reshape((2, -1)) + latlon = np.asarray(latlon, dtype=dtype).reshape((2, -1)) lat, lon = np.deg2rad(latlon) return np.array( ( @@ -675,7 +675,7 @@ def pos2latlon(pos, radius=1.0, dtype=np.double): :class:`numpy.ndarray` the 3D position array """ - pos = np.array(pos, dtype=dtype).reshape((3, -1)) + pos = np.asarray(pos, dtype=dtype).reshape((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]) diff --git a/gstools/tools/misc.py b/gstools/tools/misc.py index 53e3c9cd6..c3e994684 100755 --- a/gstools/tools/misc.py +++ b/gstools/tools/misc.py @@ -104,7 +104,7 @@ def eval_func( # care about scalar inputs func_val = 0 if func_val is None else func_val if broadcast and not callable(func_val) and np.size(func_val) == 1: - return np.array(func_val, dtype=np.double).item() + return np.asarray(func_val, dtype=np.double).item() if not callable(func_val): func_val = _func_from_single_val(func_val, dim, value_type=value_type) # care about mesh and function call @@ -112,7 +112,7 @@ def eval_func( pos, shape = format_struct_pos_dim(pos, dim) pos = generate_grid(pos) else: - pos = np.array(pos, dtype=np.double).reshape(dim, -1) + pos = np.asarray(pos, dtype=np.double).reshape(dim, -1) shape = np.shape(pos[0]) # prepend dimension if we have a vector field if value_type == "vector": @@ -125,7 +125,7 @@ def _func_from_single_val(value, dim=None, value_type="scalar"): v_d = dim if value_type == "vector" else 1 # value dim if v_d is None: # pragma: no cover raise ValueError("_func_from_single_val: dim needed for vector value.") - value = np.array(value, dtype=np.double).ravel()[:v_d] + value = np.asarray(value, dtype=np.double).ravel()[:v_d] # fill up vector valued output to dimension with last value value = np.pad( value, (0, v_d - len(value)), "constant", constant_values=value[-1] diff --git a/gstools/tools/special.py b/gstools/tools/special.py index 63940c08f..f71ae12da 100644 --- a/gstools/tools/special.py +++ b/gstools/tools/special.py @@ -8,6 +8,7 @@ .. autosummary:: inc_gamma + inc_gamma_low exp_int inc_beta tplstable_cor @@ -21,6 +22,7 @@ __all__ = [ "confidence_scaling", "inc_gamma", + "inc_gamma_low", "exp_int", "inc_beta", "tplstable_cor", @@ -64,12 +66,31 @@ def inc_gamma(s, x): if np.isclose(s, 0): return sps.exp1(x) if np.isclose(s, np.around(s)) and s < -0.5: - return x ** (s - 1) * sps.expn(int(1 - np.around(s)), x) + return x ** s * sps.expn(int(1 - np.around(s)), x) if s < 0: return (inc_gamma(s + 1, x) - x ** s * np.exp(-x)) / s return sps.gamma(s) * sps.gammaincc(s, x) +def inc_gamma_low(s, x): + r"""Calculate the lower incomplete gamma function. + + Given by: :math:`\gamma(s,x) = \int_0^x t^{s-1}\,e^{-t}\,{\rm d}t` + + Parameters + ---------- + s : :class:`float` + exponent in the integral + x : :class:`numpy.ndarray` + input values + """ + if np.isclose(s, np.around(s)) and s < 0.5: + return np.full_like(x, np.inf, dtype=np.double) + if s < 0: + return (inc_gamma_low(s + 1, x) + x ** s * np.exp(-x)) / s + return sps.gamma(s) * sps.gammainc(s, x) + + def exp_int(s, x): r"""Calculate the exponential integral :math:`E_s(x)`. @@ -86,7 +107,7 @@ def exp_int(s, x): return sps.exp1(x) if np.isclose(s, np.around(s)) and s > -0.5: return sps.expn(int(np.around(s)), x) - x = np.array(x, dtype=np.double) + x = np.asarray(x, dtype=np.double) x_neg = x < 0 x = np.abs(x) x_compare = x ** min((10, max(((1 - s), 1)))) @@ -146,7 +167,7 @@ def tplstable_cor(r, len_scale, hurst, alpha): alpha : :class:`float`, optional Shape parameter of the stable model. """ - r = np.array(np.abs(r / len_scale), dtype=np.double) + r = np.asarray(np.abs(r / len_scale), dtype=np.double) r[np.isclose(r, 0)] = 0 # hack to prevent numerical errors res = np.ones_like(r) res[r > 0] = (2 * hurst / alpha) * exp_int( @@ -179,7 +200,7 @@ def tpl_exp_spec_dens(k, dim, len_scale, hurst, len_low=0.0): spectal density of the TPLExponential model """ if np.isclose(len_low, 0.0): - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) z = (k * len_scale) ** 2 a = hurst + dim / 2.0 b = hurst + 0.5 @@ -218,14 +239,14 @@ def tpl_gau_spec_dens(k, dim, len_scale, hurst, len_low=0.0): spectal density of the TPLExponential model """ if np.isclose(len_low, 0.0): - k = np.array(k, dtype=np.double) + k = np.asarray(k, dtype=np.double) z = np.array((k * len_scale / 2.0) ** 2) res = np.empty_like(z) z_gz = z > 0.1 # greater zero z_nz = np.logical_not(z_gz) # near zero a = hurst + dim / 2.0 fac = (len_scale / 2.0) ** dim * hurst / np.pi ** (dim / 2.0) - res[z_gz] = fac * (sps.gamma(a) - inc_gamma(a, z[z_gz])) / z[z_gz] ** a + res[z_gz] = fac * inc_gamma_low(a, z[z_gz]) / z[z_gz] ** a # first order approximation for z near zero res[z_nz] = fac * (1.0 / a - z[z_nz] / (a + 1.0)) return res diff --git a/gstools/transform/__init__.py b/gstools/transform/__init__.py index 0c52fd6b2..0000a5c73 100644 --- a/gstools/transform/__init__.py +++ b/gstools/transform/__init__.py @@ -4,10 +4,20 @@ .. currentmodule:: gstools.transform -Field-Transformations +Wrapper +^^^^^^^ + +.. autosummary:: + :toctree: generated + + apply + +Field Transformations ^^^^^^^^^^^^^^^^^^^^^ .. autosummary:: + :toctree: generated + binary discrete boxcox @@ -17,11 +27,29 @@ normal_to_uniform normal_to_arcsin normal_to_uquad + apply_function + +Array Transformations +^^^^^^^^^^^^^^^^^^^^^ + +.. autosummary:: + :toctree: generated + + array_discrete + array_boxcox + array_zinnharvey + array_force_moments + array_to_lognormal + array_to_uniform + array_to_arcsin + array_to_uquad ---- """ from gstools.transform.field import ( + apply, + apply_function, binary, discrete, boxcox, @@ -32,8 +60,20 @@ normal_to_arcsin, normal_to_uquad, ) +from gstools.transform.array import ( + array_discrete, + array_boxcox, + array_zinnharvey, + array_force_moments, + array_to_lognormal, + array_to_uniform, + array_to_arcsin, + array_to_uquad, +) __all__ = [ + "apply", + "apply_function", "binary", "discrete", "boxcox", @@ -43,4 +83,12 @@ "normal_to_uniform", "normal_to_arcsin", "normal_to_uquad", + "array_discrete", + "array_boxcox", + "array_zinnharvey", + "array_force_moments", + "array_to_lognormal", + "array_to_uniform", + "array_to_arcsin", + "array_to_uquad", ] diff --git a/gstools/transform/array.py b/gstools/transform/array.py new file mode 100644 index 000000000..6ec9fd132 --- /dev/null +++ b/gstools/transform/array.py @@ -0,0 +1,351 @@ +# -*- coding: utf-8 -*- +""" +GStools subpackage providing array transformations. + +.. currentmodule:: gstools.transform.array + +The following functions are provided + +Transformations +^^^^^^^^^^^^^^^ + +.. autosummary:: + array_discrete + array_boxcox + array_zinnharvey + array_force_moments + array_to_lognormal + array_to_uniform + array_to_arcsin + array_to_uquad +""" +# pylint: disable=C0103, C0123, R0911 +from warnings import warn +import numpy as np +from scipy.special import erf, erfinv + +__all__ = [ + "array_discrete", + "array_boxcox", + "array_zinnharvey", + "array_force_moments", + "array_to_lognormal", + "array_to_uniform", + "array_to_arcsin", + "array_to_uquad", +] + + +def array_discrete( + field, values, thresholds="arithmetic", mean=None, var=None +): + """ + Discrete transformation. + + After this transformation, the field has only `len(values)` discrete + values. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + values : :any:`numpy.ndarray` + The discrete values the field will take + thresholds : :class:`str` or :any:`numpy.ndarray`, optional + the thresholds, where the value classes are separated + possible values are: + * "arithmetic": the mean of the 2 neighbouring values + * "equal": devide the field into equal parts + * an array of explicitly given thresholds + Default: "arithmetic" + mean : :class:`float`or :any:`None` + Mean of the field for "equal" thresholds. Default: np.mean(field) + var : :class:`float`or :any:`None` + Variance of the field for "equal" thresholds. Default: np.var(field) + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + if thresholds == "arithmetic": + # just in case, sort the values + values = np.sort(values) + thresholds = (values[1:] + values[:-1]) / 2 + elif thresholds == "equal": + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + values = np.asarray(values) + n = len(values) + p = np.arange(1, n) / n # n-1 equal subdivisions of [0, 1] + rescale = np.sqrt(var * 2) + # use quantile of the normal distribution to get equal ratios + thresholds = mean + rescale * erfinv(2 * p - 1) + else: + if len(values) != len(thresholds) + 1: + raise ValueError( + "discrete transformation: len(values) != len(thresholds) + 1" + ) + values = np.asarray(values) + thresholds = np.asarray(thresholds) + # check thresholds + if not np.all(thresholds[:-1] < thresholds[1:]): + raise ValueError( + "discrete transformation: thresholds need to be ascending" + ) + # use a separate result so the intermediate results are not affected + result = np.empty_like(field) + # handle edge cases + result[field <= thresholds[0]] = values[0] + result[field > thresholds[-1]] = values[-1] + for i, value in enumerate(values[1:-1]): + result[ + np.logical_and(thresholds[i] < field, field <= thresholds[i + 1]) + ] = value + return result + + +def array_boxcox(field, lmbda=1, shift=0): + """ + (Inverse) Box-Cox transformation to denormalize data. + + After this transformation, the again Box-Cox transformed field is normal + distributed. + + See: https://en.wikipedia.org/wiki/Power_transform#Box%E2%80%93Cox_transformation + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + lmbda : :class:`float`, optional + The lambda parameter of the Box-Cox transformation. + For ``lmbda=0`` one obtains the log-normal transformation. + Default: ``1`` + shift : :class:`float`, optional + The shift parameter from the two-parametric Box-Cox transformation. + The field will be shifted by that value before transformation. + Default: ``0`` + """ + field = np.asarray(field) + result = field + shift + if np.isclose(lmbda, 0): + return array_to_lognormal(result) + if np.min(result) < -1 / lmbda: + warn("Box-Cox: Some values will be cut off!") + return (np.maximum(lmbda * result + 1, 0)) ** (1 / lmbda) + + +def array_zinnharvey(field, conn="high", mean=None, var=None): + """ + Zinn and Harvey transformation to connect low or high values. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + conn : :class:`str`, optional + Desired connectivity. Either "low" or "high". + Default: "high" + mean : :class:`float` or :any:`None`, optional + Mean of the given field. If None is given, the mean will be calculated. + Default: :any:`None` + var : :class:`float` or :any:`None`, optional + Variance of the given field. + If None is given, the variance will be calculated. + Default: :any:`None` + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + result = np.abs((field - mean) / np.sqrt(var)) + result = np.sqrt(2) * erfinv(2 * erf(result / np.sqrt(2)) - 1) + if conn == "high": + result = -result + return result * np.sqrt(var) + mean + + +def array_force_moments(field, mean=0, var=1): + """ + Force moments of a normal distributed field. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + mean : :class:`float`, optional + Desired mean of the field. + Default: 0 + var : :class:`float` or :any:`None`, optional + Desired variance of the field. + Default: 1 + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + var_in = np.var(field) + mean_in = np.mean(field) + rescale = np.sqrt(var / var_in) + return rescale * (field - mean_in) + mean + + +def array_to_lognormal(field): + """ + Transform normal distribution to log-normal distribution. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + return np.exp(field) + + +def array_to_uniform(field, mean=None, var=None): + """ + Transform normal distribution to uniform distribution on [0, 1]. + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + mean : :class:`float` or :any:`None`, optional + Mean of the given field. If None is given, the mean will be calculated. + Default: :any:`None` + var : :class:`float` or :any:`None`, optional + Variance of the given field. + If None is given, the variance will be calculated. + Default: :any:`None` + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + return 0.5 * (1 + erf((field - mean) / np.sqrt(2 * var))) + + +def array_to_arcsin(field, mean=None, var=None, a=None, b=None): + """ + Transform normal distribution to arcsin distribution. + + See: https://en.wikipedia.org/wiki/Arcsine_distribution + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + mean : :class:`float` or :any:`None`, optional + Mean of the given field. If None is given, the mean will be calculated. + Default: :any:`None` + var : :class:`float` or :any:`None`, optional + Variance of the given field. + If None is given, the mean will be calculated. + Default: :any:`None` + a : :class:`float`, optional + Parameter a of the arcsin distribution (lower bound). + Default: keep mean and variance + b : :class:`float`, optional + Parameter b of the arcsin distribution (upper bound). + Default: keep mean and variance + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + a = mean - np.sqrt(2.0 * var) if a is None else float(a) + b = mean + np.sqrt(2.0 * var) if b is None else float(b) + return _uniform_to_arcsin(array_to_uniform(field, mean, var), a, b) + + +def array_to_uquad(field, mean=None, var=None, a=None, b=None): + """ + Transform normal distribution to U-quadratic distribution. + + See: https://en.wikipedia.org/wiki/U-quadratic_distribution + + Parameters + ---------- + field : :class:`numpy.ndarray` + Normal distributed values. + mean : :class:`float` or :any:`None`, optional + Mean of the given field. If None is given, the mean will be calculated. + Default: :any:`None` + var : :class:`float` or :any:`None`, optional + Variance of the given field. + If None is given, the variance will be calculated. + Default: :any:`None` + a : :class:`float`, optional + Parameter a of the U-quadratic distribution (lower bound). + Default: keep mean and variance + b : :class:`float`, optional + Parameter b of the U-quadratic distribution (upper bound). + Default: keep mean and variance + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + field = np.asarray(field) + mean = np.mean(field) if mean is None else float(mean) + var = np.var(field) if var is None else float(var) + a = mean - np.sqrt(5.0 / 3.0 * var) if a is None else float(a) + b = mean + np.sqrt(5.0 / 3.0 * var) if b is None else float(b) + return _uniform_to_uquad(array_to_uniform(field, mean, var), a, b) + + +def _uniform_to_arcsin(field, a=0, b=1): + """ + PPF of your desired distribution. + + The PPF is the inverse of the CDF and is used to sample a distribution + from uniform distributed values on [0, 1] + + in this case: the arcsin distribution + See: https://en.wikipedia.org/wiki/Arcsine_distribution + """ + field = np.asarray(field) + return (b - a) * np.sin(np.pi * 0.5 * field) ** 2 + a + + +def _uniform_to_uquad(field, a=0, b=1): + """ + PPF of your desired distribution. + + The PPF is the inverse of the CDF and is used to sample a distribution + from uniform distributed values on [0, 1] + + in this case: the U-quadratic distribution + See: https://en.wikipedia.org/wiki/U-quadratic_distribution + """ + field = np.asarray(field) + al = 12 / (b - a) ** 3 + be = (a + b) / 2 + ga = (a - b) ** 3 / 8 + y_raw = 3 * field / al + ga + result = np.zeros_like(y_raw) + result[y_raw > 0] = y_raw[y_raw > 0] ** (1 / 3) + result[y_raw < 0] = -((-y_raw[y_raw < 0]) ** (1 / 3)) + return result + be diff --git a/gstools/transform/field.py b/gstools/transform/field.py index bc0e201ee..6b9076543 100644 --- a/gstools/transform/field.py +++ b/gstools/transform/field.py @@ -6,7 +6,17 @@ The following functions are provided +Wrapper +^^^^^^^ + .. autosummary:: + apply + +Transformations +^^^^^^^^^^^^^^^ + +.. autosummary:: + apply_function binary discrete boxcox @@ -17,13 +27,27 @@ normal_to_arcsin normal_to_uquad """ -# pylint: disable=C0103 -from warnings import warn +# pylint: disable=C0103, C0123, R0911 import numpy as np -from scipy.special import erf, erfinv - +from gstools.normalizer import ( + Normalizer, + remove_trend_norm_mean, + apply_mean_norm_trend, +) +from gstools.transform.array import ( + array_discrete, + array_boxcox, + array_zinnharvey, + array_force_moments, + array_to_lognormal, + array_to_uniform, + array_to_arcsin, + array_to_uquad, +) __all__ = [ + "apply", + "apply_function", "binary", "discrete", "boxcox", @@ -36,7 +60,170 @@ ] -def binary(fld, divide=None, upper=None, lower=None): +def _pre_process(fld, data, keep_mean): + return remove_trend_norm_mean( + pos=fld.pos, + field=data, + mean=None if keep_mean else fld.mean, + normalizer=fld.normalizer, + trend=fld.trend, + mesh_type=fld.mesh_type, + value_type=fld.value_type, + check_shape=False, + ) + + +def _post_process(fld, data, keep_mean): + return apply_mean_norm_trend( + pos=fld.pos, + field=data, + mean=None if keep_mean else fld.mean, + normalizer=fld.normalizer, + trend=fld.trend, + mesh_type=fld.mesh_type, + value_type=fld.value_type, + check_shape=False, + ) + + +def _check_for_default_normal(fld): + if not type(fld.normalizer) == Normalizer: + raise ValueError( + "transform: need a normal field but there is a normalizer defined" + ) + if fld.trend is not None: + raise ValueError( + "transform: need a normal field but there is a trend defined" + ) + if callable(fld.mean) or fld.mean is None: + raise ValueError( + "transform: need a normal field but mean is not constant" + ) + + +def apply(fld, method, field="field", store=True, process=False, **kwargs): + """ + Apply field transformation. + + Parameters + ---------- + fld : :any:`Field` + Field class containing a generated field. + method : :class:`str` + Method to use. + See :any:`gstools.transform` for available transformations. + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or with a specified name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + **kwargs + Keyword arguments forwarded to selected method. + + Raises + ------ + ValueError + When method is unknown. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + kwargs["field"] = field + kwargs["store"] = store + kwargs["process"] = process + method = str(method) # ensure method is a string + if method == "binary": + return binary(fld, **kwargs) + if method == "discrete": + return discrete(fld, **kwargs) + if method == "boxcox": + return boxcox(fld, **kwargs) + if method == "zinnharvey": + return zinnharvey(fld, **kwargs) + if method.endswith("force_moments"): + return normal_force_moments(fld, **kwargs) + if method.endswith("lognormal"): + return normal_to_lognormal(fld, **kwargs) + if method.endswith("uniform"): + return normal_to_uniform(fld, **kwargs) + if method.endswith("arcsin"): + return normal_to_arcsin(fld, **kwargs) + if method.endswith("uquad"): + return normal_to_uquad(fld, **kwargs) + if method.endswith("function"): + return apply_function(fld, **kwargs) + raise ValueError(f"transform.apply: unknown method '{method}'") + + +def apply_function( + fld, + function, + field="field", + store=True, + process=False, + keep_mean=True, + **kwargs, +): + """ + Apply function as field transformation. + + Parameters + ---------- + fld : :any:`Field` + Field class containing a generated field. + function : :any:`callable` + Function to use. + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. + **kwargs + Keyword arguments forwarded to given function. + + Raises + ------ + ValueError + When function is not callable. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not callable(function): + raise ValueError("transform.apply_function: function not a 'callable'") + data = fld[field] + name, save = fld.get_store_config(store, default=field) + if process: + data = _pre_process(fld, data, keep_mean=keep_mean) + data = function(data, **kwargs) + if process: + data = _post_process(fld, data, keep_mean=keep_mean) + return fld.post_field(data, name=name, process=False, save=save) + + +def binary( + fld, + divide=None, + upper=None, + lower=None, + field="field", + store=True, + process=False, + keep_mean=True, +): """ Binary transformation. @@ -45,8 +232,7 @@ def binary(fld, divide=None, upper=None, lower=None): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. divide : :class:`float`, optional The dividing value. Default: ``fld.mean`` @@ -56,17 +242,54 @@ def binary(fld, divide=None, upper=None, lower=None): lower : :class:`float`, optional The resulting lower value of the field. Default: ``mean - sqrt(fld.model.sill)`` - """ - if fld.field is None: - print("binary: no field stored in SRF class.") - else: - divide = fld.mean if divide is None else divide - upper = fld.mean + np.sqrt(fld.model.sill) if upper is None else upper - lower = fld.mean - np.sqrt(fld.model.sill) if lower is None else lower - discrete(fld, [lower, upper], thresholds=[divide]) + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. - -def discrete(fld, values, thresholds="arithmetic"): + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not process and divide is None: + _check_for_default_normal(fld) + if divide is None: + mean = 0.0 if process and not keep_mean else fld.mean + divide = mean if divide is None else divide + upper = mean + np.sqrt(fld.model.sill) if upper is None else upper + lower = mean - np.sqrt(fld.model.sill) if lower is None else lower + kw = dict( + values=[lower, upper], + thresholds=[divide], + ) + return apply_function( + fld=fld, + function=array_discrete, + field=field, + store=store, + process=process, + keep_mean=keep_mean, + **kw, + ) + + +def discrete( + fld, + values, + thresholds="arithmetic", + field="field", + store=True, + process=False, + keep_mean=True, +): """ Discrete transformation. @@ -76,8 +299,7 @@ def discrete(fld, values, thresholds="arithmetic"): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. values : :any:`numpy.ndarray` The discrete values the field will take thresholds : :class:`str` or :any:`numpy.ndarray`, optional @@ -87,50 +309,51 @@ def discrete(fld, values, thresholds="arithmetic"): * "equal": devide the field into equal parts * an array of explicitly given thresholds Default: "arithmetic" - """ - if fld.field is None: - print("discrete: no field stored in SRF class.") - else: - if thresholds == "arithmetic": - # just in case, sort the values - values = np.sort(values) - thresholds = (values[1:] + values[:-1]) / 2 - elif thresholds == "equal": - values = np.array(values) - n = len(values) - p = np.arange(1, n) / n # n-1 equal subdivisions of [0, 1] - rescale = np.sqrt(fld.model.sill * 2) - # use quantile of the normal distribution to get equal ratios - thresholds = fld.mean + rescale * erfinv(2 * p - 1) - else: - if len(values) != len(thresholds) + 1: - raise ValueError( - "discrete transformation: " - "len(values) != len(thresholds) + 1" - ) - values = np.array(values) - thresholds = np.array(thresholds) - # check thresholds - if not np.all(thresholds[:-1] < thresholds[1:]): - raise ValueError( - "discrete transformation: thresholds need to be ascending." - ) - # use a separate result so the intermediate results are not affected - result = np.empty_like(fld.field) - # handle edge cases - result[fld.field <= thresholds[0]] = values[0] - result[fld.field > thresholds[-1]] = values[-1] - for i, value in enumerate(values[1:-1]): - result[ - np.logical_and( - thresholds[i] < fld.field, fld.field <= thresholds[i + 1] - ) - ] = value - # overwrite the field - fld.field = result - - -def boxcox(fld, lmbda=1, shift=0): + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not process and thresholds == "equal": + _check_for_default_normal(fld) + kw = dict( + values=values, + thresholds=thresholds, + mean=0.0 if process and not keep_mean else fld.mean, + var=fld.model.sill, + ) + return apply_function( + fld=fld, + function=array_discrete, + field=field, + store=store, + process=process, + keep_mean=keep_mean, + **kw, + ) + + +def boxcox( + fld, + lmbda=1, + shift=0, + field="field", + store=True, + process=False, + keep_mean=True, +): """ (Inverse) Box-Cox transformation to denormalize data. @@ -142,8 +365,7 @@ def boxcox(fld, lmbda=1, shift=0): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. lmbda : :class:`float`, optional The lambda parameter of the Box-Cox transformation. For ``lmbda=0`` one obtains the log-normal transformation. @@ -152,20 +374,43 @@ def boxcox(fld, lmbda=1, shift=0): The shift parameter from the two-parametric Box-Cox transformation. The field will be shifted by that value before transformation. Default: ``0`` - """ - if fld.field is None: - print("Box-Cox: no field stored in SRF class.") - else: - fld.mean += shift - fld.field += shift - if np.isclose(lmbda, 0): - normal_to_lognormal(fld) - if np.min(fld.field) < -1 / lmbda: - warn("Box-Cox: Some values will be cut off!") - fld.field = (np.maximum(lmbda * fld.field + 1, 0)) ** (1 / lmbda) - - -def zinnharvey(fld, conn="high"): + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. + + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + kw = dict(lmbda=lmbda, shift=shift) + return apply_function( + fld=fld, + function=array_boxcox, + field=field, + store=store, + process=process, + keep_mean=keep_mean, + **kw, + ) + + +def zinnharvey( + fld, + conn="high", + field="field", + store=True, + process=False, + keep_mean=True, +): """ Zinn and Harvey transformation to connect low or high values. @@ -174,19 +419,52 @@ def zinnharvey(fld, conn="high"): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. conn : :class:`str`, optional Desired connectivity. Either "low" or "high". Default: "high" - """ - if fld.field is None: - print("zinnharvey: no field stored in SRF class.") - else: - fld.field = _zinnharvey(fld.field, conn, fld.mean, fld.model.sill) - + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. -def normal_force_moments(fld): + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not process: + _check_for_default_normal(fld) + kw = dict( + conn=conn, + mean=0.0 if process and not keep_mean else fld.mean, + var=fld.model.sill, + ) + return apply_function( + fld=fld, + function=array_zinnharvey, + field=field, + store=store, + process=process, + keep_mean=keep_mean, + **kw, + ) + + +def normal_force_moments( + fld, + field="field", + store=True, + process=False, + keep_mean=True, +): """ Force moments of a normal distributed field. @@ -195,16 +473,43 @@ def normal_force_moments(fld): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. - """ - if fld.field is None: - print("normal_force_moments: no field stored in SRF class.") - else: - fld.field = _normal_force_moments(fld.field, fld.mean, fld.model.sill) + Field class containing a generated field. + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. - -def normal_to_lognormal(fld): + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not process: + _check_for_default_normal(fld) + kw = dict( + mean=0.0 if process and not keep_mean else fld.mean, var=fld.model.sill + ) + return apply_function( + fld=fld, + function=array_force_moments, + field=field, + store=store, + process=process, + keep_mean=keep_mean, + **kw, + ) + + +def normal_to_lognormal( + fld, field="field", store=True, process=False, keep_mean=True +): """ Transform normal distribution to log-normal distribution. @@ -213,16 +518,41 @@ def normal_to_lognormal(fld): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. - """ - if fld.field is None: - print("normal_to_lognormal: no field stored in SRF class.") - else: - fld.field = _normal_to_lognormal(fld.field) + Field class containing a generated field. + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. - -def normal_to_uniform(fld): + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + return apply_function( + fld=fld, + function=array_to_lognormal, + field=field, + store=store, + process=process, + keep_mean=keep_mean, + ) + + +def normal_to_uniform( + fld, + field="field", + store=True, + process=False, + keep_mean=True, +): """ Transform normal distribution to uniform distribution on [0, 1]. @@ -231,16 +561,36 @@ def normal_to_uniform(fld): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. - """ - if fld.field is None: - print("normal_to_uniform: no field stored in SRF class.") - else: - fld.field = _normal_to_uniform(fld.field, fld.mean, fld.model.sill) - - -def normal_to_arcsin(fld, a=None, b=None): + Field class containing a generated field. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. + """ + if not process: + _check_for_default_normal(fld) + kw = dict( + mean=0.0 if process and not keep_mean else fld.mean, var=fld.model.sill + ) + return apply_function( + fld=fld, + function=array_to_uniform, + field=field, + store=store, + process=process, + keep_mean=keep_mean, + **kw, + ) + + +def normal_to_arcsin( + fld, + a=None, + b=None, + field="field", + store=True, + process=False, + keep_mean=True, +): """ Transform normal distribution to the bimodal arcsin distribution. @@ -251,27 +601,58 @@ def normal_to_arcsin(fld, a=None, b=None): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. a : :class:`float`, optional Parameter a of the arcsin distribution (lower bound). Default: keep mean and variance b : :class:`float`, optional Parameter b of the arcsin distribution (upper bound). Default: keep mean and variance - """ - if fld.field is None: - print("normal_to_arcsin: no field stored in SRF class.") - else: - a = fld.mean - np.sqrt(2.0 * fld.model.sill) if a is None else a - b = fld.mean + np.sqrt(2.0 * fld.model.sill) if b is None else b - fld.field = _normal_to_arcsin( - fld.field, fld.mean, fld.model.sill, a, b - ) - fld.mean = (a + b) / 2.0 - + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. -def normal_to_uquad(fld, a=None, b=None): + Returns + ------- + :class:`numpy.ndarray` + Transformed field. + """ + if not process: + _check_for_default_normal(fld) + kw = dict( + mean=0.0 if process and not keep_mean else fld.mean, + var=fld.model.sill, + a=a, + b=b, + ) + return apply_function( + fld=fld, + function=array_to_arcsin, + field=field, + store=store, + process=process, + keep_mean=keep_mean, + **kw, + ) + + +def normal_to_uquad( + fld, + a=None, + b=None, + field="field", + store=True, + process=False, + keep_mean=True, +): """ Transform normal distribution to U-quadratic distribution. @@ -282,228 +663,44 @@ def normal_to_uquad(fld, a=None, b=None): Parameters ---------- fld : :any:`Field` - Spatial Random Field class containing a generated field. - Field will be transformed inplace. + Field class containing a generated field. a : :class:`float`, optional Parameter a of the U-quadratic distribution (lower bound). Default: keep mean and variance b : :class:`float`, optional Parameter b of the U-quadratic distribution (upper bound). Default: keep mean and variance - """ - if fld.field is None: - print("normal_to_uquad: no field stored in SRF class.") - else: - a = fld.mean - np.sqrt(5.0 / 3.0 * fld.model.sill) if a is None else a - b = fld.mean + np.sqrt(5.0 / 3.0 * fld.model.sill) if b is None else b - fld.field = _normal_to_uquad(fld.field, fld.mean, fld.model.sill, a, b) - fld.mean = (a + b) / 2.0 - - -# low level functions - - -def _zinnharvey(field, conn="high", mean=None, var=None): - """ - Zinn and Harvey transformation to connect low or high values. - - Parameters - ---------- - field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. - conn : :class:`str`, optional - Desired connectivity. Either "low" or "high". - Default: "high" - mean : :class:`float` or :any:`None`, optional - Mean of the given field. If None is given, the mean will be calculated. - Default: :any:`None` - var : :class:`float` or :any:`None`, optional - Variance of the given field. - If None is given, the mean will be calculated. - Default: :any:`None` + field : :class:`str`, optional + Name of field to be transformed. The default is "field". + store : :class:`str` or :class:`bool`, optional + Whether to store field inplace (True/False) or under a given name. + The default is True. + process : :class:`bool`, optional + Whether to process in/out fields with trend, normalizer and mean + of given Field instance. The default is False. + keep_mean : :class:`bool`, optional + Whether to keep the mean of the field if process=True. + The default is True. Returns ------- - :class:`numpy.ndarray` - Transformed field. - """ - if mean is None: - mean = np.mean(field) - if var is None: - var = np.var(field) - field = np.abs((field - mean) / np.sqrt(var)) - field = 2 * erf(field / np.sqrt(2)) - 1 - field = np.sqrt(2) * erfinv(field) - if conn == "high": - field = -field - return field * np.sqrt(var) + mean - - -def _normal_force_moments(field, mean=0, var=1): - """ - Force moments of a normal distributed field. - - Parameters - ---------- - field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. - mean : :class:`float`, optional - Desired mean of the field. - Default: 0 - var : :class:`float` or :any:`None`, optional - Desired variance of the field. - Default: 1 - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - var_in = np.var(field) - mean_in = np.mean(field) - rescale = np.sqrt(var / var_in) - return rescale * (field - mean_in) + mean - - -def _normal_to_lognormal(field): - """ - Transform normal distribution to log-normal distribution. - - Parameters - ---------- - field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - return np.exp(field) - - -def _normal_to_uniform(field, mean=None, var=None): - """ - Transform normal distribution to uniform distribution on [0, 1]. - - Parameters - ---------- - field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. - mean : :class:`float` or :any:`None`, optional - Mean of the given field. If None is given, the mean will be calculated. - Default: :any:`None` - var : :class:`float` or :any:`None`, optional - Variance of the given field. - If None is given, the mean will be calculated. - Default: :any:`None` - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - if mean is None: - mean = np.mean(field) - if var is None: - var = np.var(field) - return 0.5 * (1 + erf((field - mean) / np.sqrt(2 * var))) - - -def _normal_to_arcsin(field, mean=None, var=None, a=0, b=1): - """ - Transform normal distribution to arcsin distribution. - - See: https://en.wikipedia.org/wiki/Arcsine_distribution - - Parameters - ---------- - field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. - mean : :class:`float` or :any:`None`, optional - Mean of the given field. If None is given, the mean will be calculated. - Default: :any:`None` - var : :class:`float` or :any:`None`, optional - Variance of the given field. - If None is given, the mean will be calculated. - Default: :any:`None` - a : :class:`float`, optional - Parameter a of the arcsin distribution. Default: 0 - b : :class:`float`, optional - Parameter b of the arcsin distribution. Default: 1 - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - return _uniform_to_arcsin(_normal_to_uniform(field, mean, var), a, b) - - -def _normal_to_uquad(field, mean=None, var=None, a=0, b=1): - """ - Transform normal distribution to U-quadratic distribution. - - See: https://en.wikipedia.org/wiki/U-quadratic_distribution - - Parameters - ---------- - field : :class:`numpy.ndarray` - Spatial Random Field with normal distributed values. - As returned by SRF. - mean : :class:`float` or :any:`None`, optional - Mean of the given field. If None is given, the mean will be calculated. - Default: :any:`None` - var : :class:`float` or :any:`None`, optional - Variance of the given field. - If None is given, the mean will be calculated. - Default: :any:`None` - a : :class:`float`, optional - Parameter a of the U-quadratic distribution. Default: 0 - b : :class:`float`, optional - Parameter b of the U-quadratic distribution. Default: 1 - - Returns - ------- - :class:`numpy.ndarray` - Transformed field. - """ - return _uniform_to_uquad(_normal_to_uniform(field, mean, var), a, b) - - -def _uniform_to_arcsin(field, a=0, b=1): - """ - PPF of your desired distribution. - - The PPF is the inverse of the CDF and is used to sample a distribution - from uniform distributed values on [0, 1] - - in this case: the arcsin distribution - See: https://en.wikipedia.org/wiki/Arcsine_distribution - """ - return (b - a) * np.sin(np.pi * 0.5 * field) ** 2 + a - - -def _uniform_to_uquad(field, a=0, b=1): - """ - PPF of your desired distribution. - - The PPF is the inverse of the CDF and is used to sample a distribution - from uniform distributed values on [0, 1] - - in this case: the U-quadratic distribution - See: https://en.wikipedia.org/wiki/U-quadratic_distribution - """ - al = 12 / (b - a) ** 3 - be = (a + b) / 2 - ga = (a - b) ** 3 / 8 - y_raw = 3 * field / al + ga - out = np.zeros_like(y_raw) - out[y_raw > 0] = y_raw[y_raw > 0] ** (1 / 3) - out[y_raw < 0] = -((-y_raw[y_raw < 0]) ** (1 / 3)) - return out + be + :class:`numpy.ndarray` + Transformed field. + """ + if not process: + _check_for_default_normal(fld) + kw = dict( + mean=0.0 if process and not keep_mean else fld.mean, + var=fld.model.sill, + a=a, + b=b, + ) + return apply_function( + fld=fld, + function=array_to_uquad, + field=field, + store=store, + process=process, + keep_mean=keep_mean, + **kw, + ) diff --git a/gstools/variogram/binning.py b/gstools/variogram/binning.py index c7be0eca0..c2ea7a6c2 100644 --- a/gstools/variogram/binning.py +++ b/gstools/variogram/binning.py @@ -80,13 +80,13 @@ def standard_bins( if mesh_type != "unstructured": pos = generate_grid(format_struct_pos_dim(pos, dim)[0]) else: - pos = np.array(pos, dtype=np.double).reshape(dim, -1) + pos = np.asarray(pos, dtype=np.double).reshape(dim, -1) pos = latlon2pos(pos) if latlon else pos pnt_cnt = len(pos[0]) box = [] for axis in pos: box.append([np.min(axis), np.max(axis)]) - box = np.array(box) + 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 diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index c6d98e9fc..b25dacb47 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -239,7 +239,7 @@ def vario_estimate( John Wiley & Sons. (2007) """ if bin_edges is not None: - bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double) + bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double, copy=False) bin_centres = (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 @@ -292,7 +292,7 @@ def vario_estimate( # set directions dir_no = 0 if direction is not None and dim > 1: - direction = np.array(direction, ndmin=2, dtype=np.double) + direction = np.array(direction, ndmin=2, dtype=np.double, copy=False) if len(direction.shape) > 2: raise ValueError(f"Can't interpret directions: {direction}") if direction.shape[1] != dim: @@ -441,9 +441,9 @@ def vario_estimate_axis( field = np.ma.array(field, ndmin=1, dtype=np.double) if missing: field.mask = np.logical_or(field.mask, missing_mask) - mask = np.array(np.ma.getmaskarray(field), dtype=np.int32) + mask = np.asarray(np.ma.getmaskarray(field), dtype=np.int32) else: - field = np.array(field, ndmin=1, dtype=np.double) + field = np.array(field, ndmin=1, dtype=np.double, copy=False) missing_mask = None # free space axis_to_swap = AXIS_DIR[direction] if direction in AXIS else int(direction) diff --git a/tests/test_condition.py b/tests/test_condition.py index 587469d13..edb5b9669 100644 --- a/tests/test_condition.py +++ b/tests/test_condition.py @@ -1,8 +1,6 @@ # -*- coding: utf-8 -*- -""" -This is the unittest of CovModel class. -""" - +"""This is the unittest of CondSRF class.""" +from copy import copy import numpy as np import unittest import gstools as gs @@ -62,6 +60,10 @@ def test_simple(self): crf = gs.CondSRF(krige, seed=19970221) field_1 = crf.unstructured(self.pos[:dim]) field_2 = crf.structured(self.pos[:dim]) + # check reuse + raw_kr2 = copy(crf["raw_krige"]) + crf(seed=19970222) + self.assertTrue(np.allclose(raw_kr2, crf["raw_krige"])) for i, val in enumerate(self.cond_val): self.assertAlmostEqual(val, field_1[i], places=2) self.assertAlmostEqual(val, field_2[dim * (i,)], places=2) diff --git a/tests/test_covmodel.py b/tests/test_covmodel.py index 15123af20..f62fccaf3 100644 --- a/tests/test_covmodel.py +++ b/tests/test_covmodel.py @@ -260,7 +260,7 @@ def test_fitting(self): model.fit_variogram( self.gamma_x, self.gamma_y, weights=lambda x: 1 / (1 + x) ) - self.assertAlmostEqual(model.len_scale, len_save) + self.assertAlmostEqual(model.len_scale, len_save, places=6) # check ValueErrors with self.assertRaises(ValueError): model.fit_variogram(self.gamma_x, self.gamma_y, sill=2, var=3) diff --git a/tests/test_field.py b/tests/test_field.py index 53e7c3367..a5209dcb9 100644 --- a/tests/test_field.py +++ b/tests/test_field.py @@ -46,6 +46,75 @@ def test_raise(self): with self.assertRaises(ValueError): gs.field.Field(dim=3, mean=[1, 2]) + def test_pos_compare(self): + fld = gs.field.Field(dim=1) + fld.set_pos([1, 2]) + fld._dim = 2 + info = fld.set_pos([[1], [2]], info=True) + self.assertTrue(info["deleted"]) + info = fld.set_pos([[2], [3]], info=True) + self.assertTrue(info["deleted"]) + + def test_magic(self): + fld = gs.field.Field(dim=1) + f1 = np.array([0, 0], dtype=np.double) + f2 = np.array([2, 3], dtype=np.double) + fld([1, 2], store="f1") # default field with zeros + fld([1, 2], f2, store="f2") + fields1 = fld[:] + fields2 = fld[[0, 1]] + fields3 = fld[["f1", "f2"]] + fields4 = fld.all_fields + self.assertTrue(np.allclose([f1, f2], fields1)) + self.assertTrue(np.allclose([f1, f2], fields2)) + self.assertTrue(np.allclose([f1, f2], fields3)) + self.assertTrue(np.allclose([f1, f2], fields4)) + self.assertEqual(len(fld), 2) + self.assertTrue("f1" in fld) + self.assertTrue("f2" in fld) + self.assertFalse("f3" in fld) + # subscription + with self.assertRaises(KeyError): + fld["f3"] + with self.assertRaises(KeyError): + del fld["f3"] + with self.assertRaises(KeyError): + del fld[["f3"]] + del fld["f1"] + self.assertFalse("f1" in fld) + fld([1, 2], f1, store="f1") + del fld[-1] + self.assertFalse("f1" in fld) + fld([1, 2], f1, store="f1") + del fld[:] + self.assertEqual(len(fld), 0) + fld([1, 2], f1, store="f1") + del fld.field_names + self.assertEqual(len(fld), 0) + # store config (missing check) + name, save = fld.get_store_config(store="fld", fld_cnt=1) + self.assertEqual(name, ["fld"]) + self.assertTrue(save[0]) + + def test_reuse(self): + fld = gs.field.Field(dim=1) + # no pos tuple + with self.assertRaises(ValueError): + fld() + # no field shape + with self.assertRaises(ValueError): + fld.post_field([1, 2]) + # bad name + fld.set_pos([1, 2]) + with self.assertRaises(ValueError): + fld.post_field([1, 2], process=False, name=0) + # incompatible reuse + with self.assertRaises(ValueError): + fld.structured() + fld.set_pos([1, 2], "structured") + with self.assertRaises(ValueError): + fld.unstructured() + if __name__ == "__main__": unittest.main() diff --git a/tests/test_srf.py b/tests/test_srf.py index 7cca60bf4..0a83067de 100644 --- a/tests/test_srf.py +++ b/tests/test_srf.py @@ -256,46 +256,6 @@ def test_mesh_pyvista(self): _ = srf.mesh(pv_mesh, seed=self.seed, points="points") self.assertTrue(np.allclose(field, pv_mesh["field"])) - def test_transform(self): - self.cov_model.dim = 2 - srf = gs.SRF(self.cov_model, mean=self.mean, mode_no=self.mode_no) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.normal_force_moments(srf) # force ergodicity of the given field - self.assertAlmostEqual(srf.field.mean(), srf.mean) - self.assertAlmostEqual(srf.field.var(), srf.model.var) - tf.zinnharvey(srf) # make high values mostly connected - tf.normal_force_moments(srf) # force ergodicity of the given field - tf.normal_to_lognormal(srf) # log-normal - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.normal_to_arcsin(srf) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.normal_to_uquad(srf) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.normal_to_uniform(srf) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.binary(srf) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - tf.boxcox(srf) - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - values = np.linspace(np.min(srf.field), np.max(srf.field), 3) - tf.discrete(srf, values) - - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - values = [-1, 0, 1] - thresholds = [-0.9, 0.1] - tf.discrete(srf, values, thresholds) - np.testing.assert_array_equal(np.unique(srf.field), [-1, 0, 1]) - - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - values = [-1, 0, 1] - tf.discrete(srf, values, thresholds="arithmetic") - np.testing.assert_array_equal(np.unique(srf.field), [-1.0, 0.0, 1.0]) - - srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") - values = [-1, 0, 0.5, 1] - tf.discrete(srf, values, thresholds="equal") - np.testing.assert_array_equal(np.unique(srf.field), values) - def test_incomprrandmeth(self): self.cov_model = gs.Gaussian(dim=2, var=0.5, len_scale=1.0) srf = gs.SRF( diff --git a/tests/test_transform.py b/tests/test_transform.py new file mode 100644 index 000000000..641f7d013 --- /dev/null +++ b/tests/test_transform.py @@ -0,0 +1,189 @@ +# -*- coding: utf-8 -*- +"""This is the unittest of the transform submodule.""" +import unittest +import numpy as np +import gstools as gs + + +class TestSRF(unittest.TestCase): + def setUp(self): + self.cov_model = gs.Gaussian(dim=2, var=1.5, len_scale=4.0) + self.mean = 0.3 + self.mode_no = 100 + + self.seed = 825718662 + self.x_grid = np.linspace(0.0, 12.0, 48) + self.y_grid = np.linspace(0.0, 10.0, 46) + + self.methods = [ + "binary", + "boxcox", + "zinnharvey", + "normal_force_moments", + "normal_to_lognormal", + "normal_to_uniform", + "normal_to_arcsin", + "normal_to_uquad", + ] + + def test_transform_normal(self): + srf = gs.SRF(self.cov_model, mean=self.mean, mode_no=self.mode_no) + srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") + for method in self.methods: + srf.transform(method, store=method) + std = np.sqrt(srf.model.var) + self.assertTrue(set(self.methods) == set(srf.field_names[1:])) + # force moments + self.assertAlmostEqual(srf["normal_force_moments"].mean(), srf.mean) + self.assertAlmostEqual(srf["normal_force_moments"].var(), std ** 2) + # binary + np.testing.assert_allclose( + np.unique(srf.binary), srf.mean + np.array([-std, std]) + ) + # boxcox + np.testing.assert_allclose( + srf.field, gs.normalizer.BoxCox().normalize(srf.boxcox) + ) + with self.assertWarns(Warning): + srf.transform("boxcox", store="boxcox_warn", lmbda=2) + # lognormal + np.testing.assert_allclose(srf.field, np.log(srf.normal_to_lognormal)) + srf.transform("boxcox", store="boxcox2", lmbda=0) + np.testing.assert_allclose(srf.boxcox2, srf.normal_to_lognormal) + # unifrom + self.assertTrue(np.all(srf.normal_to_uniform < 1)) + self.assertTrue(np.all(srf.normal_to_uniform > 0)) + # how to test arcsin and uquad?! + + # discrete + values = [-1, 0, 1] + thresholds = [-0.9, 0.1] + srf.transform( + "discrete", values=values, thresholds=thresholds, store="f1" + ) + np.testing.assert_array_equal(np.unique(srf.f1), [-1, 0, 1]) + + values = [-1, 0, 1] + srf.transform( + "discrete", values=values, thresholds="arithmetic", store="f2" + ) + np.testing.assert_array_equal(np.unique(srf.f2), [-1.0, 0.0, 1.0]) + + values = [-1, 0, 0.5, 1] + srf.transform( + "discrete", values=values, thresholds="equal", store="f3" + ) + np.testing.assert_array_equal(np.unique(srf.f3), values) + # checks + with self.assertRaises(ValueError): + srf.transform("discrete", values=values, thresholds=[1]) + with self.assertRaises(ValueError): + srf.transform("discrete", values=values, thresholds=[1, 3, 2]) + + # function + srf.transform("function", function=lambda x: 2 * x, store="f4") + np.testing.assert_array_equal(2 * srf.field, srf.f4) + with self.assertRaises(ValueError): + srf.transform("function", function=None) + + # unknown method + with self.assertRaises(ValueError): + srf.transform("foobar") + + def test_transform_denormal(self): + srf_fail = gs.SRF( + model=self.cov_model, + mean=self.mean, + mode_no=self.mode_no, + trend=lambda x, y: x, + ) + srf_fail((self.x_grid, self.y_grid), mesh_type="structured") + with self.assertRaises(ValueError): + srf_fail.transform("zinnharvey") + + srf_fail = gs.SRF( + model=self.cov_model, + mean=lambda x, y: x, + mode_no=self.mode_no, + ) + srf_fail((self.x_grid, self.y_grid), mesh_type="structured") + with self.assertRaises(ValueError): + srf_fail.transform("zinnharvey") + + srf = gs.SRF( + model=self.cov_model, + mean=self.mean, + mode_no=self.mode_no, + normalizer=gs.normalizer.LogNormal, + ) + srf((self.x_grid, self.y_grid), seed=self.seed, mesh_type="structured") + + for method in self.methods: + if method in ("normal_to_lognormal", "boxcox"): + continue + with self.assertRaises(ValueError): + srf.transform(method, store=method) + + for method in self.methods: + srf.transform(method, store=method, process=True) + std = np.sqrt(srf.model.var) + self.assertTrue(set(self.methods) == set(srf.field_names[1:])) + # force moments + self.assertAlmostEqual( + np.log(srf["normal_force_moments"]).mean(), srf.mean + ) + self.assertAlmostEqual( + np.log(srf["normal_force_moments"]).var(), std ** 2 + ) + # binary + np.testing.assert_allclose( + np.unique(np.log(srf.binary)), srf.mean + np.array([-std, std]) + ) + # boxcox + np.testing.assert_allclose( + np.log(srf.field), + gs.normalizer.BoxCox().normalize(np.log(srf.boxcox)), + ) + # lognormal + np.testing.assert_allclose(srf.field, np.log(srf.normal_to_lognormal)) + # unifrom + self.assertTrue(np.all(np.log(srf.normal_to_uniform) < 1)) + self.assertTrue(np.all(np.log(srf.normal_to_uniform) > 0)) + + # discrete + values = [-1, 0, 1] + thresholds = [-0.9, 0.1] + srf.transform( + "discrete", + values=values, + thresholds=thresholds, + store="f1", + process=True, + ) + np.testing.assert_array_equal(np.unique(np.log(srf.f1)), [-1, 0, 1]) + + values = [-1, 0, 1] + srf.transform( + "discrete", + values=values, + thresholds="arithmetic", + store="f2", + process=True, + ) + np.testing.assert_array_equal( + np.unique(np.log(srf.f2)), [-1.0, 0.0, 1.0] + ) + + values = [-1, 0, 0.5, 1] + srf.transform( + "discrete", + values=values, + thresholds="equal", + store="f3", + process=True, + ) + np.testing.assert_array_equal(np.unique(np.log(srf.f3)), values) + + +if __name__ == "__main__": + unittest.main()