From 95b14311e4333c73d3c229c9e5a494d62aa93235 Mon Sep 17 00:00:00 2001 From: Tobias Glaubach <34371519+TobiasGlaubach@users.noreply.github.com> Date: Fri, 6 Nov 2020 16:35:27 +0100 Subject: [PATCH 01/44] Implemented directional variograms for `unstructured` (#87) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * implemented a prototype for the unstructured function with angles. * minor fix with variable name * added docstring and arguments for angle estimation * added a working version which tests the basic functionality * docstring adaption * changed for automatic testcase data generation and split up the test cases * changed default angle tolerance to 25deg * added option to also return the counts (number of pairs) from unstructured * implemented a meaningful test for 2d variogram estimation * bugfix for 3d case when elevation is 90° or 270° * implemented some basic 3d test cases * vario: cleanup cython routines; use greate-circle for tolerance in 3D; check both directions between point pairs * vario: doc update; correct intervals for angles; general formatting of angles array * vario: better handling of angle ranges * vario: fix wrong assumption about hemisphere for angles Co-authored-by: MuellerSeb --- gstools/variogram/estimator.pyx | 97 ++++++++- gstools/variogram/variogram.py | 81 ++++++-- tests/test_variogram_unstructured.py | 296 ++++++++++++++++++++++++++- 3 files changed, 449 insertions(+), 25 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 41414f6c..5c05a32d 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -10,7 +10,7 @@ import numpy as np cimport cython from cython.parallel import prange, parallel from libcpp.vector cimport vector -from libc.math cimport fabs, sqrt +from libc.math cimport fabs, sqrt, atan2, acos, asin, sin, cos, M_PI cimport numpy as np @@ -47,6 +47,68 @@ cdef inline double _distance_3d( (y[i] - y[j]) * (y[i] - y[j]) + (z[i] - z[j]) * (z[i] - z[j])) +cdef inline bint _angle_test_1d( + const double[:] x, + const double[:] y, + const double[:] z, + const double[:] angles, + const double angles_tol, + const int i, + const int j +) nogil: + return True + +cdef inline bint _angle_test_2d( + const double[:] x, + const double[:] y, + const double[:] z, + const double[:] angles, + const double angles_tol, + const int i, + const int j +) nogil: + cdef double dx = x[i] - x[j] + cdef double dy = y[i] - y[j] + # azimuth + cdef double phi1 = atan2(dy,dx) % (2.0 * M_PI) + cdef double phi2 = atan2(-dy,-dx) % (2.0 * M_PI) + # check both directions (+/-) + cdef bint dir1 = fabs(phi1 - angles[0]) <= angles_tol + cdef bint dir2 = fabs(phi2 - angles[0]) <= angles_tol + return dir1 or dir2 + +cdef inline bint _angle_test_3d( + const double[:] x, + const double[:] y, + const double[:] z, + const double[:] angles, + const double angles_tol, + const int i, + const int j +) nogil: + cdef double dx = x[i] - x[j] + cdef double dy = y[i] - y[j] + cdef double dz = z[i] - z[j] + cdef double dr = sqrt(dx**2 + dy**2 + dz**2) + # azimuth + cdef double phi1 = atan2(dy, dx) % (2.0 * M_PI) + cdef double phi2 = atan2(-dy, -dx) % (2.0 * M_PI) + # elevation + cdef double theta1 = acos(dz / dr) + cdef double theta2 = acos(-dz / dr) + # definitions for great-circle distance (for tolerance check) + cdef double dx1 = sin(theta1) * cos(phi1) - sin(angles[1]) * cos(angles[0]) + cdef double dy1 = sin(theta1) * sin(phi1) - sin(angles[1]) * sin(angles[0]) + cdef double dz1 = cos(theta1) - cos(angles[1]) + cdef double dx2 = sin(theta2) * cos(phi2) - sin(angles[1]) * cos(angles[0]) + cdef double dy2 = sin(theta2) * sin(phi2) - sin(angles[1]) * sin(angles[0]) + cdef double dz2 = cos(theta2) - cos(angles[1]) + cdef double dist1 = 2.0 * asin(sqrt(dx1**2 + dy1**2 + dz1**2) * 0.5) + cdef double dist2 = 2.0 * asin(sqrt(dx2**2 + dy2**2 + dz2**2) * 0.5) + # check both directions (+/-) + cdef bint dir1 = dist1 <= angles_tol + cdef bint dir2 = dist2 <= angles_tol + return dir1 or dir2 cdef inline double estimator_matheron(const double f_diff) nogil: return f_diff * f_diff @@ -110,6 +172,16 @@ ctypedef double (*_dist_func)( const int ) nogil +ctypedef bint (*_angle_test_func)( + const double[:], + const double[:], + const double[:], + const double[:], + const double, + const int, + const int +) nogil + def unstructured( const double[:] f, @@ -117,6 +189,8 @@ def unstructured( const double[:] x, const double[:] y=None, const double[:] z=None, + const double[:] angles=None, + const double angles_tol=0.436332, str estimator_type='m' ): if x.shape[0] != f.shape[0]: @@ -126,21 +200,35 @@ def unstructured( raise ValueError('len(bin_edges) too small') cdef _dist_func distance + cdef _angle_test_func angle_test + # 3d if z is not None: if z.shape[0] != f.shape[0]: raise ValueError('len(z) = {0} != len(f) = {1} '. format(z.shape[0], f.shape[0])) distance = _distance_3d + angle_test = _angle_test_3d # 2d elif y is not None: if y.shape[0] != f.shape[0]: raise ValueError('len(y) = {0} != len(f) = {1} '. format(y.shape[0], f.shape[0])) distance = _distance_2d + angle_test = _angle_test_2d # 1d else: distance = _distance_1d + angle_test = _angle_test_1d + + if angles is not None: + if z is not None and angles.size < 2: + raise ValueError('3d requested but only one angle given') + if y is not None and angles.size < 1: + raise ValueError('2d with angle requested but no angle given') + + if angles_tol <= 0: + raise ValueError('tolerance for angle search masks must be > 0') cdef _estimator_func estimator_func = choose_estimator_func(estimator_type) cdef _normalization_func normalization_func = ( @@ -160,11 +248,12 @@ def unstructured( for k in range(j+1, k_max): dist = distance(x, y, z, k, j) if dist >= bin_edges[i] and dist < bin_edges[i+1]: - counts[i] += 1 - variogram[i] += estimator_func(f[k] - f[j]) + if angles is None or angle_test(x, y, z, angles, angles_tol, k, j): + counts[i] += 1 + variogram[i] += estimator_func(f[k] - f[j]) normalization_func(variogram, counts) - return np.asarray(variogram) + return np.asarray(variogram), np.asarray(counts) def structured(const double[:,:,:] f, str estimator_type='m'): diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index e13ddf32..353e7736 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -37,9 +37,12 @@ def vario_estimate_unstructured( pos, field, bin_edges, + angles=None, + angles_tol=0.436332, sampling_size=None, sampling_seed=None, estimator="matheron", + return_counts=False, ): r""" Estimates the variogram on a unstructured grid. @@ -64,10 +67,6 @@ def vario_estimate_unstructured( being the bins. The Cressie estimator is more robust to outliers. - Notes - ----- - Internally uses double precision and also returns doubles. - Parameters ---------- pos : :class:`list` @@ -77,6 +76,18 @@ def vario_estimate_unstructured( the spatially distributed data bin_edges : :class:`numpy.ndarray` the bins on which the variogram will be calculated + angles : :class:`numpy.ndarray` + the angles of the main axis to calculate the variogram for in radians + angle definitions from ISO standard 80000-2:2009 + for 1d this parameter will have no effect at all + for 2d supply one angle which is azimuth φ (ccw from +x in xy plane) + for 3d supply two angles which are azimuth φ (ccw from +x in xy plane) + and inclination θ (cw from +z) + angles_tol : class:`float` + the tolerance around the variogram angle to count a point as being + within this direction from another point (the angular tolerance around + the directional vector given by angles) + Default: 25°≈0.436332 sampling_size : :class:`int` or :any:`None`, optional for large input data, this method can take a long time to compute the variogram, therefore this argument specifies @@ -92,18 +103,43 @@ def vario_estimate_unstructured( * "cressie": an estimator more robust to outliers Default: "matheron" + return_counts: class:`bool`, optional + if set to true, this function will also return the number of data + points found at each lag distance as a third return value + Default: False Returns ------- :class:`tuple` of :class:`numpy.ndarray` - the estimated variogram and the bin centers + 1. the bin centers + 2. the estimated variogram values at bin centers + 3. (optional) the number of points found at each bin center + (see argument return_counts) + + Notes + ----- + Internally uses double precision and also returns doubles. """ # TODO check_mesh field = np.array(field, ndmin=1, dtype=np.double) bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double) x, y, z, dim = pos2xyz(pos, calc_dim=True, dtype=np.double) + bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + if angles is not None and dim > 1: + # there are at most (dim-1) angles + angles = np.array(angles, ndmin=1, dtype=np.double) # type convert + angles = angles.ravel()[: (dim - 1)] # cutoff + angles = np.pad( + angles, (0, dim - angles.size - 1), "constant", constant_values=0.0 + ) # fill with 0 if too less given + # correct intervalls for angles + angles[0] = angles[0] % (2 * np.pi) + angles[1:] = angles[1:] % np.pi + elif dim == 1: + angles = None + if sampling_size is not None and sampling_size < len(field): sampled_idx = np.random.RandomState(sampling_seed).choice( np.arange(len(field)), sampling_size, replace=False @@ -117,13 +153,22 @@ def vario_estimate_unstructured( cython_estimator = _set_estimator(estimator) - return ( - bin_centres, - unstructured( - field, bin_edges, x, y, z, estimator_type=cython_estimator - ), + estimates, counts = unstructured( + field, + bin_edges, + x, + y, + z, + angles, + angles_tol, + estimator_type=cython_estimator, ) + if return_counts: + return bin_centres, estimates, counts + else: + return bin_centres, estimates + def vario_estimate_structured(field, direction="x", estimator="matheron"): r"""Estimates the variogram on a regular grid. @@ -149,14 +194,6 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"): being the bins. The Cressie estimator is more robust to outliers. - Warnings - -------- - It is assumed that the field is defined on an equidistant Cartesian grid. - - Notes - ----- - Internally uses double precision and also returns doubles. - Parameters ---------- field : :class:`numpy.ndarray` @@ -175,6 +212,14 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"): ------- :class:`numpy.ndarray` the estimated variogram along the given direction. + + Warnings + -------- + It is assumed that the field is defined on an equidistant Cartesian grid. + + Notes + ----- + Internally uses double precision and also returns doubles. """ try: mask = np.array(field.mask, dtype=np.int32) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index c5a59a85..c219ec22 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -5,12 +5,38 @@ import unittest import numpy as np -from gstools import vario_estimate_unstructured - +from gstools import vario_estimate_unstructured, Exponential, SRF +import gstools as gs class TestVariogramUnstructured(unittest.TestCase): def setUp(self): - pass + + # this code generates the testdata for the 2D rotation test cases + x = np.random.RandomState(19970221).rand(20) * 10.0 - 5 + y = np.zeros_like(x) + model = gs.Exponential(dim=2, var=2, len_scale=8) + srf = gs.SRF(model, mean=0, seed=19970221) + field = srf((x, y)) + + bins = np.arange(10) + bin_center, gamma = gs.vario_estimate_unstructured((x, ), field, bins) + idx = np.argsort(x) + self.test_data_rotation_1 = {'gamma': gamma, 'x': x[idx], 'field': field[idx], 'bins': bins, 'bin_center': bin_center} + + # CODE ABOVE SHOULD GENERATE THIS DATA + # x = np.array([ + # -4.86210059, -4.1984934 , -3.9246953 , -3.28490663, -2.16332379, + # -1.87553275, -1.74125124, -1.27224687, -1.20931578, -0.2413368 , + # 0.03200921, 1.17099773, 1.53863105, 1.64478688, 2.75252136, + # 3.3556915 , 3.89828775, 4.21485964, 4.5364357 , 4.79236969]), + # field = np.array([ + # -1.10318365, -0.53566629, -0.41789049, -1.06167529, 0.38449961, + # -0.36550477, -0.98905552, -0.19352766, 0.16264266, 0.26920833, + # 0.05379665, 0.71275006, 0.36651935, 0.17366865, 1.20022343, + # 0.79385446, 0.69456069, 1.0733393 , 0.71191592, 0.71969766]) + # gamma_exp = np.array([ + # 0.14260989, 0.18301197, 0.25855841, 0.29990083, 0.67914526, + # 0.60136535, 0.92875492, 1.46910435, 1.10165104]) def test_doubles(self): x = np.arange(1, 11, 1, dtype=np.double) @@ -217,6 +243,270 @@ def test_assertions(self): ValueError, vario_estimate_unstructured, [x], field_e, bins ) + def test_angles_2D_x2x(self): + + x = self.test_data_rotation_1['x'] + field = self.test_data_rotation_1['field'] + gamma_exp = self.test_data_rotation_1['gamma'] + bins = self.test_data_rotation_1['bins'] + y = np.zeros_like(x) + + # test case 1.) + # all along x axis on x axis + + bin_centres, gamma = vario_estimate_unstructured( + (x, y), + field, + bins, + angles=[0] + ) + + for i in range(gamma.size): + self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) + + def test_angles_2D_y2x(self): + + x = self.test_data_rotation_1['x'] + field = self.test_data_rotation_1['field'] + gamma_exp = self.test_data_rotation_1['gamma'] + bins = self.test_data_rotation_1['bins'] + y = np.zeros_like(x) + + # test case 2.) + # all along y axis on y axis but calculation for x axis + + bin_centres, gamma = vario_estimate_unstructured( + (y, x), + field, + bins, + angles=[0] + ) + + for i in range(gamma.size): + self.assertAlmostEqual(0, gamma[i], places=3) + + def test_angles_2D_y2y(self): + + x = self.test_data_rotation_1['x'] + field = self.test_data_rotation_1['field'] + gamma_exp = self.test_data_rotation_1['gamma'] + bins = self.test_data_rotation_1['bins'] + y = np.zeros_like(x) + + # test case 3.) + # all along y axis on y axis and calculation for y axis + + bin_centres, gamma = vario_estimate_unstructured( + (y, x), + field, + bins, + angles=[np.pi/2.] + ) + + for i in range(gamma.size): + self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) + + def test_angles_2D_xy2x(self): + + x = self.test_data_rotation_1['x'] + field = self.test_data_rotation_1['field'] + gamma_exp = self.test_data_rotation_1['gamma'] + bins = self.test_data_rotation_1['bins'] + y = np.zeros_like(x) + + + # test case 4.) + # data along 45deg axis but calculation for x axis + + ccos, csin = np.cos(np.pi/4.), np.sin(np.pi/4.) + + xr = [xx * ccos - yy * csin for xx, yy in zip(x, y)] + yr = [xx * csin + yy * ccos for xx, yy in zip(x, y)] + + bin_centres, gamma = vario_estimate_unstructured( + (xr, yr), + field, + bins, + angles=[0] + ) + + for i in range(gamma.size): + self.assertAlmostEqual(0, gamma[i], places=3) + + def test_angles_2D_estim(self): + + seed = gs.random.MasterRNG(19970221) + rng = np.random.RandomState(seed()) + rng = np.random + x = rng.randint(0, 100, size=3000) + y = rng.randint(0, 100, size=3000) + + model = gs.Exponential(dim=2, var=1, len_scale=[12, 3], angles=np.pi / 8) + model_maj = gs.Exponential(dim=1, var=1, len_scale=[12]) + model_min = gs.Exponential(dim=1, var=1, len_scale=[3]) + + srf = gs.SRF(model, seed=20170519) + field = srf((x, y)) + + + bins = np.arange(0, 50, 2.5) + angle_mask = 22.5 + angle_tol = 22.5 + + bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured( + (x, y), + field, + bins, + angles=[np.deg2rad(angle_mask)], + angles_tol=np.deg2rad(angle_tol) + ) + + bin_centers_min, gamma_min = gs.vario_estimate_unstructured( + (x, y), + field, + bins, + angles=[np.deg2rad(angle_mask + 90.0)], + angles_tol=np.deg2rad(angle_tol) + ) + + gamma_maj_real = model_maj.variogram(bin_centers_maj) + gamma_min_real = model_min.variogram(bin_centers_min) + + # we have no real way of testing values, but we can test some basic properties which definitelly need to be true + # test that the major estimate aligns better with major real than minor real + self.assertTrue(np.sum((gamma_maj_real - gamma_maj)**2) < np.sum((gamma_min_real - gamma_maj)**2)) + # test that the minor estimate aligns better with minor real than major real + self.assertTrue(np.sum((gamma_min_real - gamma_min)**2) < np.sum((gamma_maj_real - gamma_min)**2)) + # test that both variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance + self.assertTrue((np.mean(gamma_min[-5:]) - np.var(field)) / np.var(field) < 0.1) + self.assertTrue((np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1) + + def test_angles_line_3D(self): + + x = self.test_data_rotation_1['x'] + field = self.test_data_rotation_1['field'] + gamma_exp = self.test_data_rotation_1['gamma'] + bins = self.test_data_rotation_1['bins'] + + def test_xyz(x, y, z, angles, gamma_exp): + bin_centres, gamma = vario_estimate_unstructured( + (x, y, z), + field, + bins, + angles=angles + ) + + if np.ndim(gamma_exp) == 0: + gamma_exp = np.ones_like(gamma) * gamma_exp + + for i in range(gamma.size): + self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) + + # all along x axis and calculation for x axis + test_xyz(x, np.zeros_like(x), np.zeros_like(x), [0], gamma_exp) + # all along y axis and calculation for x axis + test_xyz(np.zeros_like(x), x, np.zeros_like(x), [0], 0) + # all along z axis and calculation for x axis + test_xyz(np.zeros_like(x), np.zeros_like(x), x, [0], 0) + + angles_rot_azim_90 = [np.pi/2] + # all along x axis and calculation for y axis + test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90, 0) + # all along y axis and calculation for y axis + test_xyz(np.zeros_like(x), x, np.zeros_like(x), angles_rot_azim_90, gamma_exp) + # all along z axis and calculation for y axis + test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90, 0) + + # for elevation it is important to check, that IF elevation is 90° or 270° it does + # not matter how we rotated before, since any rotation around z (in XY plane) + # followed by a rotation around x' (in YZ' plane) by 90° will result in the same + # coordinates, (when the structure is two points with zero extend) + + # test with [0, 90] + angles_rot_azim_90_elev_90 = [0, np.pi/2] + # all along x axis and calculation for z axis + test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90_elev_90, 0) + # all along y axis and calculation for z axis + test_xyz(np.zeros_like(x), x, np.zeros_like(x), angles_rot_azim_90_elev_90, 0) + # all along z axis and calculation for z axis + test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90_elev_90, gamma_exp) + + # test with [90, 90] + angles_rot_azim_90_elev_90 = [np.pi/2, np.pi/2] + # all along x axis and calculation for z axis + test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90_elev_90, 0) + # all along y axis and calculation for z axis + test_xyz(np.zeros_like(x), x, np.zeros_like(x), angles_rot_azim_90_elev_90, 0) + # all along z axis and calculation for z axis + test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90_elev_90, gamma_exp) + + def test_angles_3D_estim(self): + + seed = gs.random.MasterRNG(19970221) + rng = np.random.RandomState(seed()) + rng = np.random + x = rng.randint(0, 50, size=1000) + y = rng.randint(0, 50, size=1000) + z = rng.randint(0, 50, size=1000) + + model = gs.Exponential(dim=3, var=1, len_scale=[10, 5, 2], angles=[np.pi / 8, np.pi / 16]) + model_maj = gs.Exponential(dim=1, var=1, len_scale=[10]) + model_min1 = gs.Exponential(dim=1, var=1, len_scale=[5]) + model_min2 = gs.Exponential(dim=1, var=1, len_scale=[2]) + + srf = gs.SRF(model, seed=20170519) + field = srf((x, y, z)) + + bins = np.arange(0, 25, 5) + angle_mask = [np.pi / 8, np.pi / 16] + angle_tol = 22.5 + + # in x' + bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured( + (x, y, z), + field, + bins, + angles=[angle_mask[0], angle_mask[1]], + angles_tol=np.deg2rad(angle_tol) + ) + + # in y' + bin_centers_min1, gamma_min1 = gs.vario_estimate_unstructured( + (x, y, z), + field, + bins, + angles=[angle_mask[0] + 0.5*np.pi, angle_mask[1]], + angles_tol=np.deg2rad(angle_tol) + ) + + # in z' + bin_centers_min2, gamma_min2 = gs.vario_estimate_unstructured( + (x, y, z), + field, + bins, + angles=[angle_mask[0], angle_mask[1] + 0.5*np.pi], + angles_tol=np.deg2rad(angle_tol) + ) + + gamma_maj_real = model_maj.variogram(bin_centers_maj) + gamma_min1_real = model_min1.variogram(bin_centers_min1) + gamma_min2_real = model_min2.variogram(bin_centers_min2) + + # we have no real way of testing values, but we can test some basic properties which definitelly need to be true + # test that the major estimate aligns better with major real than the minor reals + self.assertTrue(np.sum((gamma_maj_real - gamma_maj)**2) < np.sum((gamma_min1_real - gamma_maj)**2)) + self.assertTrue(np.sum((gamma_maj_real - gamma_maj)**2) < np.sum((gamma_min2_real - gamma_maj)**2)) + # test that the minor1 estimate aligns better with minor1 real than major real and minor2 real + self.assertTrue(np.sum((gamma_min1_real - gamma_min1)**2) < np.sum((gamma_maj_real - gamma_min1)**2)) + self.assertTrue(np.sum((gamma_min1_real - gamma_min1)**2) < np.sum((gamma_min2_real - gamma_min1)**2)) + # test that the minor2 estimate aligns better with minor2 real than major real and minor1 real + self.assertTrue(np.sum((gamma_min2_real - gamma_min2)**2) < np.sum((gamma_maj_real - gamma_min2)**2)) + self.assertTrue(np.sum((gamma_min2_real - gamma_min2)**2) < np.sum((gamma_min1_real - gamma_min2)**2)) + # test that all variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance + self.assertTrue((np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1) + self.assertTrue((np.mean(gamma_min1[-5:]) - np.var(field)) / np.var(field) < 0.1) + self.assertTrue((np.mean(gamma_min2[-5:]) - np.var(field)) / np.var(field) < 0.1) + if __name__ == "__main__": unittest.main() From 2ab750431db5dfa611009ec598cf4f3aa9bebb80 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 6 Nov 2020 19:42:19 +0100 Subject: [PATCH 02/44] blackened --- gstools/covmodel/base.py | 4 +--- gstools/krige/base.py | 4 +++- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/gstools/covmodel/base.py b/gstools/covmodel/base.py index 16ff4845..cbaa1d30 100644 --- a/gstools/covmodel/base.py +++ b/gstools/covmodel/base.py @@ -1219,9 +1219,7 @@ def name(self): @property def do_rotation(self): """:any:`bool`: State if a rotation is performed.""" - return ( - not np.all(np.isclose(self.angles, 0.0)) - ) + return not np.all(np.isclose(self.angles, 0.0)) @property def is_isotropic(self): diff --git a/gstools/krige/base.py b/gstools/krige/base.py index 8f04d83b..ff4b6a41 100755 --- a/gstools/krige/base.py +++ b/gstools/krige/base.py @@ -634,7 +634,9 @@ def __repr__(self): """Return String representation.""" return ( "{0}(model={1}, cond_no={2}".format( - self.name, self.model.name, self.cond_no, + self.name, + self.model.name, + self.cond_no, ) + ")" ) From 0ed1c15251411601aed70d2694c56ad8f0809682 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Fri, 6 Nov 2020 19:42:42 +0100 Subject: [PATCH 03/44] tests: blackened --- tests/test_variogram_unstructured.py | 304 ++++++++++++++++----------- 1 file changed, 187 insertions(+), 117 deletions(-) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index c219ec22..fee201e9 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -8,6 +8,7 @@ from gstools import vario_estimate_unstructured, Exponential, SRF import gstools as gs + class TestVariogramUnstructured(unittest.TestCase): def setUp(self): @@ -19,9 +20,15 @@ def setUp(self): field = srf((x, y)) bins = np.arange(10) - bin_center, gamma = gs.vario_estimate_unstructured((x, ), field, bins) + bin_center, gamma = gs.vario_estimate_unstructured((x,), field, bins) idx = np.argsort(x) - self.test_data_rotation_1 = {'gamma': gamma, 'x': x[idx], 'field': field[idx], 'bins': bins, 'bin_center': bin_center} + self.test_data_rotation_1 = { + "gamma": gamma, + "x": x[idx], + "field": field[idx], + "bins": bins, + "bin_center": bin_center, + } # CODE ABOVE SHOULD GENERATE THIS DATA # x = np.array([ @@ -245,94 +252,81 @@ def test_assertions(self): def test_angles_2D_x2x(self): - x = self.test_data_rotation_1['x'] - field = self.test_data_rotation_1['field'] - gamma_exp = self.test_data_rotation_1['gamma'] - bins = self.test_data_rotation_1['bins'] + x = self.test_data_rotation_1["x"] + field = self.test_data_rotation_1["field"] + gamma_exp = self.test_data_rotation_1["gamma"] + bins = self.test_data_rotation_1["bins"] y = np.zeros_like(x) - + # test case 1.) # all along x axis on x axis - + bin_centres, gamma = vario_estimate_unstructured( - (x, y), - field, - bins, - angles=[0] + (x, y), field, bins, angles=[0] ) - + for i in range(gamma.size): self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) def test_angles_2D_y2x(self): - x = self.test_data_rotation_1['x'] - field = self.test_data_rotation_1['field'] - gamma_exp = self.test_data_rotation_1['gamma'] - bins = self.test_data_rotation_1['bins'] + x = self.test_data_rotation_1["x"] + field = self.test_data_rotation_1["field"] + gamma_exp = self.test_data_rotation_1["gamma"] + bins = self.test_data_rotation_1["bins"] y = np.zeros_like(x) # test case 2.) # all along y axis on y axis but calculation for x axis - + bin_centres, gamma = vario_estimate_unstructured( - (y, x), - field, - bins, - angles=[0] + (y, x), field, bins, angles=[0] ) - + for i in range(gamma.size): self.assertAlmostEqual(0, gamma[i], places=3) - + def test_angles_2D_y2y(self): - x = self.test_data_rotation_1['x'] - field = self.test_data_rotation_1['field'] - gamma_exp = self.test_data_rotation_1['gamma'] - bins = self.test_data_rotation_1['bins'] + x = self.test_data_rotation_1["x"] + field = self.test_data_rotation_1["field"] + gamma_exp = self.test_data_rotation_1["gamma"] + bins = self.test_data_rotation_1["bins"] y = np.zeros_like(x) # test case 3.) # all along y axis on y axis and calculation for y axis - + bin_centres, gamma = vario_estimate_unstructured( - (y, x), - field, - bins, - angles=[np.pi/2.] + (y, x), field, bins, angles=[np.pi / 2.0] ) - + for i in range(gamma.size): self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) def test_angles_2D_xy2x(self): - x = self.test_data_rotation_1['x'] - field = self.test_data_rotation_1['field'] - gamma_exp = self.test_data_rotation_1['gamma'] - bins = self.test_data_rotation_1['bins'] + x = self.test_data_rotation_1["x"] + field = self.test_data_rotation_1["field"] + gamma_exp = self.test_data_rotation_1["gamma"] + bins = self.test_data_rotation_1["bins"] y = np.zeros_like(x) - # test case 4.) # data along 45deg axis but calculation for x axis - - ccos, csin = np.cos(np.pi/4.), np.sin(np.pi/4.) - + + ccos, csin = np.cos(np.pi / 4.0), np.sin(np.pi / 4.0) + xr = [xx * ccos - yy * csin for xx, yy in zip(x, y)] yr = [xx * csin + yy * ccos for xx, yy in zip(x, y)] - + bin_centres, gamma = vario_estimate_unstructured( - (xr, yr), - field, - bins, - angles=[0] + (xr, yr), field, bins, angles=[0] ) - + for i in range(gamma.size): self.assertAlmostEqual(0, gamma[i], places=3) - + def test_angles_2D_estim(self): seed = gs.random.MasterRNG(19970221) @@ -341,59 +335,67 @@ def test_angles_2D_estim(self): x = rng.randint(0, 100, size=3000) y = rng.randint(0, 100, size=3000) - model = gs.Exponential(dim=2, var=1, len_scale=[12, 3], angles=np.pi / 8) + model = gs.Exponential( + dim=2, var=1, len_scale=[12, 3], angles=np.pi / 8 + ) model_maj = gs.Exponential(dim=1, var=1, len_scale=[12]) model_min = gs.Exponential(dim=1, var=1, len_scale=[3]) srf = gs.SRF(model, seed=20170519) field = srf((x, y)) - bins = np.arange(0, 50, 2.5) angle_mask = 22.5 angle_tol = 22.5 bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured( - (x, y), - field, - bins, - angles=[np.deg2rad(angle_mask)], - angles_tol=np.deg2rad(angle_tol) - ) + (x, y), + field, + bins, + angles=[np.deg2rad(angle_mask)], + angles_tol=np.deg2rad(angle_tol), + ) bin_centers_min, gamma_min = gs.vario_estimate_unstructured( - (x, y), - field, - bins, - angles=[np.deg2rad(angle_mask + 90.0)], - angles_tol=np.deg2rad(angle_tol) - ) + (x, y), + field, + bins, + angles=[np.deg2rad(angle_mask + 90.0)], + angles_tol=np.deg2rad(angle_tol), + ) gamma_maj_real = model_maj.variogram(bin_centers_maj) gamma_min_real = model_min.variogram(bin_centers_min) # we have no real way of testing values, but we can test some basic properties which definitelly need to be true # test that the major estimate aligns better with major real than minor real - self.assertTrue(np.sum((gamma_maj_real - gamma_maj)**2) < np.sum((gamma_min_real - gamma_maj)**2)) + self.assertTrue( + np.sum((gamma_maj_real - gamma_maj) ** 2) + < np.sum((gamma_min_real - gamma_maj) ** 2) + ) # test that the minor estimate aligns better with minor real than major real - self.assertTrue(np.sum((gamma_min_real - gamma_min)**2) < np.sum((gamma_maj_real - gamma_min)**2)) + self.assertTrue( + np.sum((gamma_min_real - gamma_min) ** 2) + < np.sum((gamma_maj_real - gamma_min) ** 2) + ) # test that both variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance - self.assertTrue((np.mean(gamma_min[-5:]) - np.var(field)) / np.var(field) < 0.1) - self.assertTrue((np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1) + self.assertTrue( + (np.mean(gamma_min[-5:]) - np.var(field)) / np.var(field) < 0.1 + ) + self.assertTrue( + (np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1 + ) def test_angles_line_3D(self): - x = self.test_data_rotation_1['x'] - field = self.test_data_rotation_1['field'] - gamma_exp = self.test_data_rotation_1['gamma'] - bins = self.test_data_rotation_1['bins'] + x = self.test_data_rotation_1["x"] + field = self.test_data_rotation_1["field"] + gamma_exp = self.test_data_rotation_1["gamma"] + bins = self.test_data_rotation_1["bins"] def test_xyz(x, y, z, angles, gamma_exp): bin_centres, gamma = vario_estimate_unstructured( - (x, y, z), - field, - bins, - angles=angles + (x, y, z), field, bins, angles=angles ) if np.ndim(gamma_exp) == 0: @@ -409,36 +411,78 @@ def test_xyz(x, y, z, angles, gamma_exp): # all along z axis and calculation for x axis test_xyz(np.zeros_like(x), np.zeros_like(x), x, [0], 0) - angles_rot_azim_90 = [np.pi/2] + angles_rot_azim_90 = [np.pi / 2] # all along x axis and calculation for y axis test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90, 0) # all along y axis and calculation for y axis - test_xyz(np.zeros_like(x), x, np.zeros_like(x), angles_rot_azim_90, gamma_exp) + test_xyz( + np.zeros_like(x), + x, + np.zeros_like(x), + angles_rot_azim_90, + gamma_exp, + ) # all along z axis and calculation for y axis test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90, 0) - # for elevation it is important to check, that IF elevation is 90° or 270° it does - # not matter how we rotated before, since any rotation around z (in XY plane) + # for elevation it is important to check, that IF elevation is 90° or 270° it does + # not matter how we rotated before, since any rotation around z (in XY plane) # followed by a rotation around x' (in YZ' plane) by 90° will result in the same # coordinates, (when the structure is two points with zero extend) # test with [0, 90] - angles_rot_azim_90_elev_90 = [0, np.pi/2] + angles_rot_azim_90_elev_90 = [0, np.pi / 2] # all along x axis and calculation for z axis - test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90_elev_90, 0) + test_xyz( + x, + np.zeros_like(x), + np.zeros_like(x), + angles_rot_azim_90_elev_90, + 0, + ) # all along y axis and calculation for z axis - test_xyz(np.zeros_like(x), x, np.zeros_like(x), angles_rot_azim_90_elev_90, 0) + test_xyz( + np.zeros_like(x), + x, + np.zeros_like(x), + angles_rot_azim_90_elev_90, + 0, + ) # all along z axis and calculation for z axis - test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90_elev_90, gamma_exp) + test_xyz( + np.zeros_like(x), + np.zeros_like(x), + x, + angles_rot_azim_90_elev_90, + gamma_exp, + ) - # test with [90, 90] - angles_rot_azim_90_elev_90 = [np.pi/2, np.pi/2] + # test with [90, 90] + angles_rot_azim_90_elev_90 = [np.pi / 2, np.pi / 2] # all along x axis and calculation for z axis - test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90_elev_90, 0) + test_xyz( + x, + np.zeros_like(x), + np.zeros_like(x), + angles_rot_azim_90_elev_90, + 0, + ) # all along y axis and calculation for z axis - test_xyz(np.zeros_like(x), x, np.zeros_like(x), angles_rot_azim_90_elev_90, 0) + test_xyz( + np.zeros_like(x), + x, + np.zeros_like(x), + angles_rot_azim_90_elev_90, + 0, + ) # all along z axis and calculation for z axis - test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90_elev_90, gamma_exp) + test_xyz( + np.zeros_like(x), + np.zeros_like(x), + x, + angles_rot_azim_90_elev_90, + gamma_exp, + ) def test_angles_3D_estim(self): @@ -449,7 +493,9 @@ def test_angles_3D_estim(self): y = rng.randint(0, 50, size=1000) z = rng.randint(0, 50, size=1000) - model = gs.Exponential(dim=3, var=1, len_scale=[10, 5, 2], angles=[np.pi / 8, np.pi / 16]) + model = gs.Exponential( + dim=3, var=1, len_scale=[10, 5, 2], angles=[np.pi / 8, np.pi / 16] + ) model_maj = gs.Exponential(dim=1, var=1, len_scale=[10]) model_min1 = gs.Exponential(dim=1, var=1, len_scale=[5]) model_min2 = gs.Exponential(dim=1, var=1, len_scale=[2]) @@ -463,30 +509,30 @@ def test_angles_3D_estim(self): # in x' bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured( - (x, y, z), - field, - bins, - angles=[angle_mask[0], angle_mask[1]], - angles_tol=np.deg2rad(angle_tol) - ) + (x, y, z), + field, + bins, + angles=[angle_mask[0], angle_mask[1]], + angles_tol=np.deg2rad(angle_tol), + ) # in y' bin_centers_min1, gamma_min1 = gs.vario_estimate_unstructured( - (x, y, z), - field, - bins, - angles=[angle_mask[0] + 0.5*np.pi, angle_mask[1]], - angles_tol=np.deg2rad(angle_tol) - ) + (x, y, z), + field, + bins, + angles=[angle_mask[0] + 0.5 * np.pi, angle_mask[1]], + angles_tol=np.deg2rad(angle_tol), + ) # in z' bin_centers_min2, gamma_min2 = gs.vario_estimate_unstructured( - (x, y, z), - field, - bins, - angles=[angle_mask[0], angle_mask[1] + 0.5*np.pi], - angles_tol=np.deg2rad(angle_tol) - ) + (x, y, z), + field, + bins, + angles=[angle_mask[0], angle_mask[1] + 0.5 * np.pi], + angles_tol=np.deg2rad(angle_tol), + ) gamma_maj_real = model_maj.variogram(bin_centers_maj) gamma_min1_real = model_min1.variogram(bin_centers_min1) @@ -494,19 +540,43 @@ def test_angles_3D_estim(self): # we have no real way of testing values, but we can test some basic properties which definitelly need to be true # test that the major estimate aligns better with major real than the minor reals - self.assertTrue(np.sum((gamma_maj_real - gamma_maj)**2) < np.sum((gamma_min1_real - gamma_maj)**2)) - self.assertTrue(np.sum((gamma_maj_real - gamma_maj)**2) < np.sum((gamma_min2_real - gamma_maj)**2)) + self.assertTrue( + np.sum((gamma_maj_real - gamma_maj) ** 2) + < np.sum((gamma_min1_real - gamma_maj) ** 2) + ) + self.assertTrue( + np.sum((gamma_maj_real - gamma_maj) ** 2) + < np.sum((gamma_min2_real - gamma_maj) ** 2) + ) # test that the minor1 estimate aligns better with minor1 real than major real and minor2 real - self.assertTrue(np.sum((gamma_min1_real - gamma_min1)**2) < np.sum((gamma_maj_real - gamma_min1)**2)) - self.assertTrue(np.sum((gamma_min1_real - gamma_min1)**2) < np.sum((gamma_min2_real - gamma_min1)**2)) + self.assertTrue( + np.sum((gamma_min1_real - gamma_min1) ** 2) + < np.sum((gamma_maj_real - gamma_min1) ** 2) + ) + self.assertTrue( + np.sum((gamma_min1_real - gamma_min1) ** 2) + < np.sum((gamma_min2_real - gamma_min1) ** 2) + ) # test that the minor2 estimate aligns better with minor2 real than major real and minor1 real - self.assertTrue(np.sum((gamma_min2_real - gamma_min2)**2) < np.sum((gamma_maj_real - gamma_min2)**2)) - self.assertTrue(np.sum((gamma_min2_real - gamma_min2)**2) < np.sum((gamma_min1_real - gamma_min2)**2)) + self.assertTrue( + np.sum((gamma_min2_real - gamma_min2) ** 2) + < np.sum((gamma_maj_real - gamma_min2) ** 2) + ) + self.assertTrue( + np.sum((gamma_min2_real - gamma_min2) ** 2) + < np.sum((gamma_min1_real - gamma_min2) ** 2) + ) # test that all variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance - self.assertTrue((np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1) - self.assertTrue((np.mean(gamma_min1[-5:]) - np.var(field)) / np.var(field) < 0.1) - self.assertTrue((np.mean(gamma_min2[-5:]) - np.var(field)) / np.var(field) < 0.1) - + self.assertTrue( + (np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1 + ) + self.assertTrue( + (np.mean(gamma_min1[-5:]) - np.var(field)) / np.var(field) < 0.1 + ) + self.assertTrue( + (np.mean(gamma_min2[-5:]) - np.var(field)) / np.var(field) < 0.1 + ) + if __name__ == "__main__": unittest.main() From 04c0a71af99e8c3e5e718b996629ee6924c6d4f3 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 14:10:53 +0100 Subject: [PATCH 04/44] tools.geometric.ang2dir: new converter from angles to directional vector --- gstools/tools/geometric.py | 40 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index 6edcd758..1faa5f56 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -161,3 +161,43 @@ def xyz2pos(x, y=None, z=None, dtype=None, max_dim=3): if z is not None and max_dim > 2: pos.append(np.array(z, dtype=dtype).reshape(-1)) return tuple(pos) + + +def ang2dir(angles, dtype=np.double, dim=None): + """Convert n-D spherical coordinates to Euclidean direction vectors. + + Parameters + ---------- + angles : :class:`list` of :class:`numpy.ndarray` + spherical coordinates given as angles. + dtype : data-type, optional + The desired data-type for the array. + If not given, then the type will be determined as the minimum type + required to hold the objects in the sequence. Default: None + dim : :class:`int`, optional + Cut of information above the given dimension. + Otherwise, dimension is determined by number of angles + Default: None + + Returns + ------- + :class:`numpy.ndarray` + the list of direction vectors + """ + angles = np.array(angles, ndmin=2, dtype=dtype) + if len(angles.shape) > 2: + raise ValueError("Can't interpret angles array {}".format(angles)) + dim = angles.shape[1] + 1 if dim is None else dim + if dim == 2 and angles.shape[0] == 1: + # fix for 2D where only one angle per direction is given + angles = angles.T # can't be interpreted if dim=None is given + if dim != angles.shape[1] + 1 or dim == 1: + raise ValueError("Wrong dim. ({}) for angles {}".format(dim, angles)) + vec = np.empty((angles.shape[0], dim), dtype=dtype) + vec[:, 0] = np.prod(np.sin(angles), axis=1) + for i in range(1, dim): + vec[:, i] = np.prod(np.sin(angles[:, i:]), axis=1) # empty prod = 1 + vec[:, i] *= np.cos(angles[:, (i - 1)]) + if dim in [2, 3]: + vec[:, [0, 1]] = vec[:, [1, 0]] # to match convention in 2D and 3D + return vec From 73f1b6800f71558698efc6bdb4da39c09d2acc07 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 14:13:04 +0100 Subject: [PATCH 05/44] vario: introduce directional vectors+bandwith; simplify cython routines (pos array; nD dist; angle checks) --- gstools/variogram/estimator.pyx | 285 ++++++++++++++++---------------- gstools/variogram/variogram.py | 128 +++++++++----- 2 files changed, 226 insertions(+), 187 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 5c05a32d..447926df 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -10,7 +10,7 @@ import numpy as np cimport cython from cython.parallel import prange, parallel from libcpp.vector cimport vector -from libc.math cimport fabs, sqrt, atan2, acos, asin, sin, cos, M_PI +from libc.math cimport fabs, sqrt, acos, M_PI cimport numpy as np @@ -18,97 +18,54 @@ DTYPE = np.double ctypedef np.double_t DTYPE_t -cdef inline double _distance_1d( - const double[:] x, - const double[:] y, - const double[:] z, +cdef inline double distance( + const int dim, + const double[:,:] pos, const int i, const int j ) nogil: - return sqrt((x[i] - x[j]) * (x[i] - x[j])) + cdef int d + cdef double dist_squared = 0.0 + for d in range(dim): + dist_squared += ((pos[d,i] - pos[d,j]) * (pos[d,i] - pos[d,j])) + return sqrt(dist_squared) -cdef inline double _distance_2d( - const double[:] x, - const double[:] y, - const double[:] z, - const int i, - const int j -) nogil: - return sqrt((x[i] - x[j]) * (x[i] - x[j]) + (y[i] - y[j]) * (y[i] - y[j])) -cdef inline double _distance_3d( - const double[:] x, - const double[:] y, - const double[:] z, - const int i, - const int j -) nogil: - return sqrt((x[i] - x[j]) * (x[i] - x[j]) + - (y[i] - y[j]) * (y[i] - y[j]) + - (z[i] - z[j]) * (z[i] - z[j])) - -cdef inline bint _angle_test_1d( - const double[:] x, - const double[:] y, - const double[:] z, - const double[:] angles, +cdef inline bint dir_test( + const int dim, + const double[:,:] pos, + const double[:,:] direction, const double angles_tol, + const double bandwith, const int i, - const int j + const int j, + const int d ) nogil: - return True + cdef double s_prod = 0.0 + cdef double p_norm = 0.0 + cdef double b_dist = 0.0 + cdef int k + cdef bint in_band = True + cdef bint in_angle + + # scalar-product calculation for bandwith projection and angle calculation + for k in range(dim): + s_prod += (pos[k,i] - pos[k,j]) * direction[d,k] + p_norm += (pos[k,i] - pos[k,j])**2 + p_norm = sqrt(p_norm) + + # calculate band-distance by projection of point-pair-vec to direction line + if bandwith > 0.0: + for k in range(dim): + b_dist += ((pos[k,i] - pos[k,j]) - s_prod * direction[d,k]) ** 2 + b_dist = sqrt(b_dist) + in_band = b_dist < bandwith + + # use smallest angle by taking absolut value for arccos angle formula + in_angle = acos(fabs(s_prod) / p_norm) < angles_tol + + return in_band and in_angle -cdef inline bint _angle_test_2d( - const double[:] x, - const double[:] y, - const double[:] z, - const double[:] angles, - const double angles_tol, - const int i, - const int j -) nogil: - cdef double dx = x[i] - x[j] - cdef double dy = y[i] - y[j] - # azimuth - cdef double phi1 = atan2(dy,dx) % (2.0 * M_PI) - cdef double phi2 = atan2(-dy,-dx) % (2.0 * M_PI) - # check both directions (+/-) - cdef bint dir1 = fabs(phi1 - angles[0]) <= angles_tol - cdef bint dir2 = fabs(phi2 - angles[0]) <= angles_tol - return dir1 or dir2 - -cdef inline bint _angle_test_3d( - const double[:] x, - const double[:] y, - const double[:] z, - const double[:] angles, - const double angles_tol, - const int i, - const int j -) nogil: - cdef double dx = x[i] - x[j] - cdef double dy = y[i] - y[j] - cdef double dz = z[i] - z[j] - cdef double dr = sqrt(dx**2 + dy**2 + dz**2) - # azimuth - cdef double phi1 = atan2(dy, dx) % (2.0 * M_PI) - cdef double phi2 = atan2(-dy, -dx) % (2.0 * M_PI) - # elevation - cdef double theta1 = acos(dz / dr) - cdef double theta2 = acos(-dz / dr) - # definitions for great-circle distance (for tolerance check) - cdef double dx1 = sin(theta1) * cos(phi1) - sin(angles[1]) * cos(angles[0]) - cdef double dy1 = sin(theta1) * sin(phi1) - sin(angles[1]) * sin(angles[0]) - cdef double dz1 = cos(theta1) - cos(angles[1]) - cdef double dx2 = sin(theta2) * cos(phi2) - sin(angles[1]) * cos(angles[0]) - cdef double dy2 = sin(theta2) * sin(phi2) - sin(angles[1]) * sin(angles[0]) - cdef double dz2 = cos(theta2) - cos(angles[1]) - cdef double dist1 = 2.0 * asin(sqrt(dx1**2 + dy1**2 + dz1**2) * 0.5) - cdef double dist2 = 2.0 * asin(sqrt(dx2**2 + dy2**2 + dz2**2) * 0.5) - # check both directions (+/-) - cdef bint dir1 = dist1 <= angles_tol - cdef bint dir2 = dist2 <= angles_tol - return dir1 or dir2 cdef inline double estimator_matheron(const double f_diff) nogil: return f_diff * f_diff @@ -148,6 +105,38 @@ ctypedef void (*_normalization_func)( vector[long]& ) +cdef inline void normalization_matheron_vec( + double[:,:]& variogram, + long[:,:]& counts +): + cdef int d, i + for d in range(variogram.shape[0]): + for i in range(variogram.shape[1]): + # avoid division by zero + if counts[d, i] == 0: + counts[d, i] = 1 + variogram[d, i] /= (2. * counts[d, i]) + +cdef inline void normalization_cressie_vec( + double[:,:]& variogram, + long[:,:]& counts +): + cdef int d, i + for d in range(variogram.shape[0]): + for i in range(variogram.shape[1]): + # avoid division by zero + if counts[d, i] == 0: + counts[d, i] = 1 + variogram[d, i] = ( + 0.5 * (1./counts[d, i] * variogram[d, i])**4 / + (0.457 + 0.494 / counts[d, i] + 0.045 / counts[d, i]**2) + ) + +ctypedef void (*_normalization_func_vec)( + double[:,:]&, + long[:,:]& +) + cdef _estimator_func choose_estimator_func(str estimator_type): cdef _estimator_func estimator_func if estimator_type == 'm': @@ -164,93 +153,102 @@ cdef _normalization_func choose_estimator_normalization(str estimator_type): normalization_func = normalization_cressie return normalization_func -ctypedef double (*_dist_func)( - const double[:], - const double[:], - const double[:], - const int, - const int -) nogil - -ctypedef bint (*_angle_test_func)( - const double[:], - const double[:], - const double[:], - const double[:], - const double, - const int, - const int -) nogil +cdef _normalization_func_vec choose_estimator_normalization_vec(str estimator_type): + cdef _normalization_func_vec normalization_func_vec + if estimator_type == 'm': + normalization_func_vec = normalization_matheron_vec + elif estimator_type == 'c': + normalization_func_vec = normalization_cressie_vec + return normalization_func_vec -def unstructured( +def directional( + const int dim, const double[:] f, const double[:] bin_edges, - const double[:] x, - const double[:] y=None, - const double[:] z=None, - const double[:] angles=None, - const double angles_tol=0.436332, + const double[:,:] pos, + const double[:,:] direction, # should be normed + const double angles_tol=M_PI/8.0, + const double bandwith=-1.0, # negative values to turn of bandwith search str estimator_type='m' ): - if x.shape[0] != f.shape[0]: - raise ValueError('len(x) = {0} != len(f) = {1} '. - format(x.shape[0], f.shape[0])) + if pos.shape[1] != f.shape[0]: + raise ValueError('len(pos) = {0} != len(f) = {1} '. + format(pos.shape[1], f.shape[0])) + if bin_edges.shape[0] < 2: raise ValueError('len(bin_edges) too small') - cdef _dist_func distance - cdef _angle_test_func angle_test - - # 3d - if z is not None: - if z.shape[0] != f.shape[0]: - raise ValueError('len(z) = {0} != len(f) = {1} '. - format(z.shape[0], f.shape[0])) - distance = _distance_3d - angle_test = _angle_test_3d - # 2d - elif y is not None: - if y.shape[0] != f.shape[0]: - raise ValueError('len(y) = {0} != len(f) = {1} '. - format(y.shape[0], f.shape[0])) - distance = _distance_2d - angle_test = _angle_test_2d - # 1d - else: - distance = _distance_1d - angle_test = _angle_test_1d - - if angles is not None: - if z is not None and angles.size < 2: - raise ValueError('3d requested but only one angle given') - if y is not None and angles.size < 1: - raise ValueError('2d with angle requested but no angle given') - if angles_tol <= 0: raise ValueError('tolerance for angle search masks must be > 0') + cdef _estimator_func estimator_func = choose_estimator_func(estimator_type) + cdef _normalization_func_vec normalization_func_vec = ( + choose_estimator_normalization_vec(estimator_type) + ) + + cdef int d_max = direction.shape[0] + cdef int i_max = bin_edges.shape[0] - 1 + cdef int j_max = pos.shape[1] - 1 + cdef int k_max = pos.shape[1] + + cdef double[:,:] variogram = np.zeros((d_max, len(bin_edges)-1)) + cdef long[:,:] counts = np.zeros((d_max, len(bin_edges)-1), dtype=long) + cdef vector[double] pos1 = vector[double](dim, 0.0) + cdef vector[double] pos2 = vector[double](dim, 0.0) + cdef int i, j, k, d + cdef DTYPE_t dist + + for i in prange(i_max, nogil=True): + for j in range(j_max): + for k in range(j+1, k_max): + dist = distance(dim, pos, j, k) + if dist >= bin_edges[i] and dist < bin_edges[i+1]: + for d in range(d_max): + if dir_test(dim, pos, direction, angles_tol, bandwith, k, j, d): + counts[d, i] += 1 + variogram[d, i] += estimator_func(f[k] - f[j]) + + normalization_func_vec(variogram, counts) + return np.asarray(variogram), np.asarray(counts) + +def unstructured( + const int dim, + const double[:] f, + const double[:] bin_edges, + const double[:,:] pos, + str estimator_type='m' +): + if pos.shape[1] != f.shape[0]: + raise ValueError('len(pos) = {0} != len(f) = {1} '. + format(pos.shape[1], f.shape[0])) + + if bin_edges.shape[0] < 2: + raise ValueError('len(bin_edges) too small') + cdef _estimator_func estimator_func = choose_estimator_func(estimator_type) cdef _normalization_func normalization_func = ( choose_estimator_normalization(estimator_type) ) cdef int i_max = bin_edges.shape[0] - 1 - cdef int j_max = x.shape[0] - 1 - cdef int k_max = x.shape[0] + cdef int j_max = pos.shape[1] - 1 + cdef int k_max = pos.shape[1] cdef vector[double] variogram = vector[double](len(bin_edges)-1, 0.0) cdef vector[long] counts = vector[long](len(bin_edges)-1, 0) + cdef vector[double] pos1 = vector[double](dim, 0.0) + cdef vector[double] pos2 = vector[double](dim, 0.0) cdef int i, j, k cdef DTYPE_t dist + for i in prange(i_max, nogil=True): for j in range(j_max): for k in range(j+1, k_max): - dist = distance(x, y, z, k, j) + dist = distance(dim, pos, j, k) if dist >= bin_edges[i] and dist < bin_edges[i+1]: - if angles is None or angle_test(x, y, z, angles, angles_tol, k, j): - counts[i] += 1 - variogram[i] += estimator_func(f[k] - f[j]) + counts[i] += 1 + variogram[i] += estimator_func(f[k] - f[j]) normalization_func(variogram, counts) return np.asarray(variogram), np.asarray(counts) @@ -282,6 +280,7 @@ def structured(const double[:,:,:] f, str estimator_type='m'): normalization_func(variogram, counts) return np.asarray(variogram) + def ma_structured( const double[:,:,:] f, const bint[:,:,:] mask, diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 353e7736..0b5d8229 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -14,8 +14,13 @@ import numpy as np -from gstools.tools.geometric import pos2xyz -from gstools.variogram.estimator import unstructured, structured, ma_structured +from gstools.tools.geometric import pos2xyz, xyz2pos, ang2dir +from gstools.variogram.estimator import ( + unstructured, + structured, + ma_structured, + directional, +) __all__ = ["vario_estimate_unstructured", "vario_estimate_structured"] @@ -37,8 +42,10 @@ def vario_estimate_unstructured( pos, field, bin_edges, + direction=None, angles=None, - angles_tol=0.436332, + angles_tol=np.pi / 8, + bandwith=None, sampling_size=None, sampling_seed=None, estimator="matheron", @@ -76,18 +83,35 @@ def vario_estimate_unstructured( the spatially distributed data bin_edges : :class:`numpy.ndarray` the bins on which the variogram will be calculated - angles : :class:`numpy.ndarray` + direction : :class:`list` of :class:`numpy.ndarray`, optional + directions to evaluate a directional variogram. + Anglular tolerance is given by `angles_tol`. + Bandwith to cut off how wide the search for point pairs should be + is given by `bandwith`. + You can provide multiple directions at once to get one variogram + for each direction. + For a single direction you can also use the `angles` parameter, + to provide the direction by its spherical coordianates. + Default: :any:`None` + angles : :class:`numpy.ndarray`, optional the angles of the main axis to calculate the variogram for in radians angle definitions from ISO standard 80000-2:2009 for 1d this parameter will have no effect at all for 2d supply one angle which is azimuth φ (ccw from +x in xy plane) for 3d supply two angles which are azimuth φ (ccw from +x in xy plane) - and inclination θ (cw from +z) - angles_tol : class:`float` + and inclination θ (cw from +z). + Can be used instead of direction for a single main axis. + Default: :any:`None` + angles_tol : class:`float`, optional the tolerance around the variogram angle to count a point as being within this direction from another point (the angular tolerance around the directional vector given by angles) - Default: 25°≈0.436332 + Default: `np.pi/8` = 22.5° + bandwith : class:`float`, optional + Bandwith to cut off the angular tolerance for directional variograms. + If None is given, only the `angles_tol` parameter will control the + point selection. + Default: :any:`None` sampling_size : :class:`int` or :any:`None`, optional for large input data, this method can take a long time to compute the variogram, therefore this argument specifies @@ -124,50 +148,66 @@ def vario_estimate_unstructured( field = np.array(field, ndmin=1, dtype=np.double) bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double) x, y, z, dim = pos2xyz(pos, calc_dim=True, dtype=np.double) - bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 - - if angles is not None and dim > 1: - # there are at most (dim-1) angles - angles = np.array(angles, ndmin=1, dtype=np.double) # type convert - angles = angles.ravel()[: (dim - 1)] # cutoff - angles = np.pad( - angles, (0, dim - angles.size - 1), "constant", constant_values=0.0 - ) # fill with 0 if too less given - # correct intervalls for angles - angles[0] = angles[0] % (2 * np.pi) - angles[1:] = angles[1:] % np.pi - elif dim == 1: - angles = None - + # initialize number of directions + dir_no = 0 + if direction is not None and dim > 1: + direction = np.array(direction, ndmin=2, dtype=np.double) + if len(direction.shape) > 2: + raise ValueError("Can't interpret directions {}".format(direction)) + if direction.shape[1] != dim: + raise ValueError("Can't interpret directions {}".format(direction)) + dir_no = direction.shape[0] + # convert given angles to direction vector + if angles is not None and direction is None and dim > 1: + direction = ang2dir(angles=angles, dtype=np.double, dim=dim) + dir_no = direction.shape[0] + # prepare directional variogram + if dir_no > 0: + norms = np.linalg.norm(direction, axis=1) + if np.any(np.isclose(norms, 0)): + raise ValueError("Zero length direction {}".format(direction)) + # only unit-vectors for directions + direction = np.divide(direction, norms[:, np.newaxis]) + # negative bandwith to turn it off + bandwith = float(bandwith) if bandwith is not None else -1.0 + angles_tol = float(angles_tol) + # prepare positions + pos = np.array(xyz2pos(x, y, z, dtype=np.double, max_dim=dim)) + # prepare sampled variogram if sampling_size is not None and sampling_size < len(field): sampled_idx = np.random.RandomState(sampling_seed).choice( np.arange(len(field)), sampling_size, replace=False ) field = field[sampled_idx] - x = x[sampled_idx] - if dim > 1: - y = y[sampled_idx] - if dim > 2: - z = z[sampled_idx] - + pos = pos[:, sampled_idx] + # select variogram estimator cython_estimator = _set_estimator(estimator) - - estimates, counts = unstructured( - field, - bin_edges, - x, - y, - z, - angles, - angles_tol, - estimator_type=cython_estimator, - ) - + # run + if dir_no == 0: + estimates, counts = unstructured( + dim, + field, + bin_edges, + pos, + estimator_type=cython_estimator, + ) + else: + estimates, counts = directional( + dim, + field, + bin_edges, + pos, + direction, + angles_tol, + bandwith, + estimator_type=cython_estimator, + ) + if dir_no == 1: + estimates, counts = estimates[0], counts[0] if return_counts: return bin_centres, estimates, counts - else: - return bin_centres, estimates + return bin_centres, estimates def vario_estimate_structured(field, direction="x", estimator="matheron"): @@ -246,10 +286,10 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"): cython_estimator = _set_estimator(estimator) # fill up the field with empty dimensions up to a number of 3 - for i in range(3 - len(field.shape)): + for __ in range(3 - len(field.shape)): field = field[..., np.newaxis] if masked: - for i in range(3 - len(mask.shape)): + for __ in range(3 - len(mask.shape)): mask = mask[..., np.newaxis] if mask is None: From 90bdd995ebb0050038a3ef1af3a3ce6f52439662 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 14:15:49 +0100 Subject: [PATCH 06/44] tests/vario: skip directional tests for now --- tests/test_variogram_unstructured.py | 640 +++++++++++++-------------- 1 file changed, 320 insertions(+), 320 deletions(-) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index fee201e9..02fd30e5 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -250,332 +250,332 @@ def test_assertions(self): ValueError, vario_estimate_unstructured, [x], field_e, bins ) - def test_angles_2D_x2x(self): + # def test_angles_2D_x2x(self): - x = self.test_data_rotation_1["x"] - field = self.test_data_rotation_1["field"] - gamma_exp = self.test_data_rotation_1["gamma"] - bins = self.test_data_rotation_1["bins"] - y = np.zeros_like(x) - - # test case 1.) - # all along x axis on x axis - - bin_centres, gamma = vario_estimate_unstructured( - (x, y), field, bins, angles=[0] - ) - - for i in range(gamma.size): - self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) + # x = self.test_data_rotation_1["x"] + # field = self.test_data_rotation_1["field"] + # gamma_exp = self.test_data_rotation_1["gamma"] + # bins = self.test_data_rotation_1["bins"] + # y = np.zeros_like(x) - def test_angles_2D_y2x(self): - - x = self.test_data_rotation_1["x"] - field = self.test_data_rotation_1["field"] - gamma_exp = self.test_data_rotation_1["gamma"] - bins = self.test_data_rotation_1["bins"] - y = np.zeros_like(x) - - # test case 2.) - # all along y axis on y axis but calculation for x axis - - bin_centres, gamma = vario_estimate_unstructured( - (y, x), field, bins, angles=[0] - ) + # # test case 1.) + # # all along x axis on x axis - for i in range(gamma.size): - self.assertAlmostEqual(0, gamma[i], places=3) - - def test_angles_2D_y2y(self): - - x = self.test_data_rotation_1["x"] - field = self.test_data_rotation_1["field"] - gamma_exp = self.test_data_rotation_1["gamma"] - bins = self.test_data_rotation_1["bins"] - y = np.zeros_like(x) - - # test case 3.) - # all along y axis on y axis and calculation for y axis - - bin_centres, gamma = vario_estimate_unstructured( - (y, x), field, bins, angles=[np.pi / 2.0] - ) - - for i in range(gamma.size): - self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) - - def test_angles_2D_xy2x(self): - - x = self.test_data_rotation_1["x"] - field = self.test_data_rotation_1["field"] - gamma_exp = self.test_data_rotation_1["gamma"] - bins = self.test_data_rotation_1["bins"] - y = np.zeros_like(x) + # bin_centres, gamma = vario_estimate_unstructured( + # (x, y), field, bins, angles=[0] + # ) - # test case 4.) - # data along 45deg axis but calculation for x axis + # for i in range(gamma.size): + # self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) - ccos, csin = np.cos(np.pi / 4.0), np.sin(np.pi / 4.0) + # def test_angles_2D_y2x(self): - xr = [xx * ccos - yy * csin for xx, yy in zip(x, y)] - yr = [xx * csin + yy * ccos for xx, yy in zip(x, y)] - - bin_centres, gamma = vario_estimate_unstructured( - (xr, yr), field, bins, angles=[0] - ) - - for i in range(gamma.size): - self.assertAlmostEqual(0, gamma[i], places=3) - - def test_angles_2D_estim(self): - - seed = gs.random.MasterRNG(19970221) - rng = np.random.RandomState(seed()) - rng = np.random - x = rng.randint(0, 100, size=3000) - y = rng.randint(0, 100, size=3000) - - model = gs.Exponential( - dim=2, var=1, len_scale=[12, 3], angles=np.pi / 8 - ) - model_maj = gs.Exponential(dim=1, var=1, len_scale=[12]) - model_min = gs.Exponential(dim=1, var=1, len_scale=[3]) - - srf = gs.SRF(model, seed=20170519) - field = srf((x, y)) - - bins = np.arange(0, 50, 2.5) - angle_mask = 22.5 - angle_tol = 22.5 - - bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured( - (x, y), - field, - bins, - angles=[np.deg2rad(angle_mask)], - angles_tol=np.deg2rad(angle_tol), - ) - - bin_centers_min, gamma_min = gs.vario_estimate_unstructured( - (x, y), - field, - bins, - angles=[np.deg2rad(angle_mask + 90.0)], - angles_tol=np.deg2rad(angle_tol), - ) - - gamma_maj_real = model_maj.variogram(bin_centers_maj) - gamma_min_real = model_min.variogram(bin_centers_min) - - # we have no real way of testing values, but we can test some basic properties which definitelly need to be true - # test that the major estimate aligns better with major real than minor real - self.assertTrue( - np.sum((gamma_maj_real - gamma_maj) ** 2) - < np.sum((gamma_min_real - gamma_maj) ** 2) - ) - # test that the minor estimate aligns better with minor real than major real - self.assertTrue( - np.sum((gamma_min_real - gamma_min) ** 2) - < np.sum((gamma_maj_real - gamma_min) ** 2) - ) - # test that both variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance - self.assertTrue( - (np.mean(gamma_min[-5:]) - np.var(field)) / np.var(field) < 0.1 - ) - self.assertTrue( - (np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1 - ) - - def test_angles_line_3D(self): - - x = self.test_data_rotation_1["x"] - field = self.test_data_rotation_1["field"] - gamma_exp = self.test_data_rotation_1["gamma"] - bins = self.test_data_rotation_1["bins"] - - def test_xyz(x, y, z, angles, gamma_exp): - bin_centres, gamma = vario_estimate_unstructured( - (x, y, z), field, bins, angles=angles - ) - - if np.ndim(gamma_exp) == 0: - gamma_exp = np.ones_like(gamma) * gamma_exp - - for i in range(gamma.size): - self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) - - # all along x axis and calculation for x axis - test_xyz(x, np.zeros_like(x), np.zeros_like(x), [0], gamma_exp) - # all along y axis and calculation for x axis - test_xyz(np.zeros_like(x), x, np.zeros_like(x), [0], 0) - # all along z axis and calculation for x axis - test_xyz(np.zeros_like(x), np.zeros_like(x), x, [0], 0) - - angles_rot_azim_90 = [np.pi / 2] - # all along x axis and calculation for y axis - test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90, 0) - # all along y axis and calculation for y axis - test_xyz( - np.zeros_like(x), - x, - np.zeros_like(x), - angles_rot_azim_90, - gamma_exp, - ) - # all along z axis and calculation for y axis - test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90, 0) - - # for elevation it is important to check, that IF elevation is 90° or 270° it does - # not matter how we rotated before, since any rotation around z (in XY plane) - # followed by a rotation around x' (in YZ' plane) by 90° will result in the same - # coordinates, (when the structure is two points with zero extend) - - # test with [0, 90] - angles_rot_azim_90_elev_90 = [0, np.pi / 2] - # all along x axis and calculation for z axis - test_xyz( - x, - np.zeros_like(x), - np.zeros_like(x), - angles_rot_azim_90_elev_90, - 0, - ) - # all along y axis and calculation for z axis - test_xyz( - np.zeros_like(x), - x, - np.zeros_like(x), - angles_rot_azim_90_elev_90, - 0, - ) - # all along z axis and calculation for z axis - test_xyz( - np.zeros_like(x), - np.zeros_like(x), - x, - angles_rot_azim_90_elev_90, - gamma_exp, - ) - - # test with [90, 90] - angles_rot_azim_90_elev_90 = [np.pi / 2, np.pi / 2] - # all along x axis and calculation for z axis - test_xyz( - x, - np.zeros_like(x), - np.zeros_like(x), - angles_rot_azim_90_elev_90, - 0, - ) - # all along y axis and calculation for z axis - test_xyz( - np.zeros_like(x), - x, - np.zeros_like(x), - angles_rot_azim_90_elev_90, - 0, - ) - # all along z axis and calculation for z axis - test_xyz( - np.zeros_like(x), - np.zeros_like(x), - x, - angles_rot_azim_90_elev_90, - gamma_exp, - ) - - def test_angles_3D_estim(self): - - seed = gs.random.MasterRNG(19970221) - rng = np.random.RandomState(seed()) - rng = np.random - x = rng.randint(0, 50, size=1000) - y = rng.randint(0, 50, size=1000) - z = rng.randint(0, 50, size=1000) - - model = gs.Exponential( - dim=3, var=1, len_scale=[10, 5, 2], angles=[np.pi / 8, np.pi / 16] - ) - model_maj = gs.Exponential(dim=1, var=1, len_scale=[10]) - model_min1 = gs.Exponential(dim=1, var=1, len_scale=[5]) - model_min2 = gs.Exponential(dim=1, var=1, len_scale=[2]) - - srf = gs.SRF(model, seed=20170519) - field = srf((x, y, z)) - - bins = np.arange(0, 25, 5) - angle_mask = [np.pi / 8, np.pi / 16] - angle_tol = 22.5 - - # in x' - bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured( - (x, y, z), - field, - bins, - angles=[angle_mask[0], angle_mask[1]], - angles_tol=np.deg2rad(angle_tol), - ) - - # in y' - bin_centers_min1, gamma_min1 = gs.vario_estimate_unstructured( - (x, y, z), - field, - bins, - angles=[angle_mask[0] + 0.5 * np.pi, angle_mask[1]], - angles_tol=np.deg2rad(angle_tol), - ) - - # in z' - bin_centers_min2, gamma_min2 = gs.vario_estimate_unstructured( - (x, y, z), - field, - bins, - angles=[angle_mask[0], angle_mask[1] + 0.5 * np.pi], - angles_tol=np.deg2rad(angle_tol), - ) - - gamma_maj_real = model_maj.variogram(bin_centers_maj) - gamma_min1_real = model_min1.variogram(bin_centers_min1) - gamma_min2_real = model_min2.variogram(bin_centers_min2) - - # we have no real way of testing values, but we can test some basic properties which definitelly need to be true - # test that the major estimate aligns better with major real than the minor reals - self.assertTrue( - np.sum((gamma_maj_real - gamma_maj) ** 2) - < np.sum((gamma_min1_real - gamma_maj) ** 2) - ) - self.assertTrue( - np.sum((gamma_maj_real - gamma_maj) ** 2) - < np.sum((gamma_min2_real - gamma_maj) ** 2) - ) - # test that the minor1 estimate aligns better with minor1 real than major real and minor2 real - self.assertTrue( - np.sum((gamma_min1_real - gamma_min1) ** 2) - < np.sum((gamma_maj_real - gamma_min1) ** 2) - ) - self.assertTrue( - np.sum((gamma_min1_real - gamma_min1) ** 2) - < np.sum((gamma_min2_real - gamma_min1) ** 2) - ) - # test that the minor2 estimate aligns better with minor2 real than major real and minor1 real - self.assertTrue( - np.sum((gamma_min2_real - gamma_min2) ** 2) - < np.sum((gamma_maj_real - gamma_min2) ** 2) - ) - self.assertTrue( - np.sum((gamma_min2_real - gamma_min2) ** 2) - < np.sum((gamma_min1_real - gamma_min2) ** 2) - ) - # test that all variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance - self.assertTrue( - (np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1 - ) - self.assertTrue( - (np.mean(gamma_min1[-5:]) - np.var(field)) / np.var(field) < 0.1 - ) - self.assertTrue( - (np.mean(gamma_min2[-5:]) - np.var(field)) / np.var(field) < 0.1 - ) + # x = self.test_data_rotation_1["x"] + # field = self.test_data_rotation_1["field"] + # gamma_exp = self.test_data_rotation_1["gamma"] + # bins = self.test_data_rotation_1["bins"] + # y = np.zeros_like(x) + + # # test case 2.) + # # all along y axis on y axis but calculation for x axis + + # bin_centres, gamma = vario_estimate_unstructured( + # (y, x), field, bins, angles=[0] + # ) + + # for i in range(gamma.size): + # self.assertAlmostEqual(0, gamma[i], places=3) + + # def test_angles_2D_y2y(self): + + # x = self.test_data_rotation_1["x"] + # field = self.test_data_rotation_1["field"] + # gamma_exp = self.test_data_rotation_1["gamma"] + # bins = self.test_data_rotation_1["bins"] + # y = np.zeros_like(x) + + # # test case 3.) + # # all along y axis on y axis and calculation for y axis + + # bin_centres, gamma = vario_estimate_unstructured( + # (y, x), field, bins, angles=[np.pi / 2.0] + # ) + + # for i in range(gamma.size): + # self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) + + # def test_angles_2D_xy2x(self): + + # x = self.test_data_rotation_1["x"] + # field = self.test_data_rotation_1["field"] + # gamma_exp = self.test_data_rotation_1["gamma"] + # bins = self.test_data_rotation_1["bins"] + # y = np.zeros_like(x) + + # # test case 4.) + # # data along 45deg axis but calculation for x axis + + # ccos, csin = np.cos(np.pi / 4.0), np.sin(np.pi / 4.0) + + # xr = [xx * ccos - yy * csin for xx, yy in zip(x, y)] + # yr = [xx * csin + yy * ccos for xx, yy in zip(x, y)] + + # bin_centres, gamma = vario_estimate_unstructured( + # (xr, yr), field, bins, angles=[0] + # ) + + # for i in range(gamma.size): + # self.assertAlmostEqual(0, gamma[i], places=3) + + # def test_angles_2D_estim(self): + + # seed = gs.random.MasterRNG(19970221) + # rng = np.random.RandomState(seed()) + # rng = np.random + # x = rng.randint(0, 100, size=3000) + # y = rng.randint(0, 100, size=3000) + + # model = gs.Exponential( + # dim=2, var=1, len_scale=[12, 3], angles=np.pi / 8 + # ) + # model_maj = gs.Exponential(dim=1, var=1, len_scale=[12]) + # model_min = gs.Exponential(dim=1, var=1, len_scale=[3]) + + # srf = gs.SRF(model, seed=20170519) + # field = srf((x, y)) + + # bins = np.arange(0, 50, 2.5) + # angle_mask = 22.5 + # angle_tol = 22.5 + + # bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured( + # (x, y), + # field, + # bins, + # angles=[np.deg2rad(angle_mask)], + # angles_tol=np.deg2rad(angle_tol), + # ) + + # bin_centers_min, gamma_min = gs.vario_estimate_unstructured( + # (x, y), + # field, + # bins, + # angles=[np.deg2rad(angle_mask + 90.0)], + # angles_tol=np.deg2rad(angle_tol), + # ) + + # gamma_maj_real = model_maj.variogram(bin_centers_maj) + # gamma_min_real = model_min.variogram(bin_centers_min) + + # # we have no real way of testing values, but we can test some basic properties which definitelly need to be true + # # test that the major estimate aligns better with major real than minor real + # self.assertTrue( + # np.sum((gamma_maj_real - gamma_maj) ** 2) + # < np.sum((gamma_min_real - gamma_maj) ** 2) + # ) + # # test that the minor estimate aligns better with minor real than major real + # self.assertTrue( + # np.sum((gamma_min_real - gamma_min) ** 2) + # < np.sum((gamma_maj_real - gamma_min) ** 2) + # ) + # # test that both variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance + # self.assertTrue( + # (np.mean(gamma_min[-5:]) - np.var(field)) / np.var(field) < 0.1 + # ) + # self.assertTrue( + # (np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1 + # ) + + # def test_angles_line_3D(self): + + # x = self.test_data_rotation_1["x"] + # field = self.test_data_rotation_1["field"] + # gamma_exp = self.test_data_rotation_1["gamma"] + # bins = self.test_data_rotation_1["bins"] + + # def test_xyz(x, y, z, angles, gamma_exp): + # bin_centres, gamma = vario_estimate_unstructured( + # (x, y, z), field, bins, angles=angles + # ) + + # if np.ndim(gamma_exp) == 0: + # gamma_exp = np.ones_like(gamma) * gamma_exp + + # for i in range(gamma.size): + # self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) + + # # all along x axis and calculation for x axis + # test_xyz(x, np.zeros_like(x), np.zeros_like(x), [0], gamma_exp) + # # all along y axis and calculation for x axis + # test_xyz(np.zeros_like(x), x, np.zeros_like(x), [0], 0) + # # all along z axis and calculation for x axis + # test_xyz(np.zeros_like(x), np.zeros_like(x), x, [0], 0) + + # angles_rot_azim_90 = [np.pi / 2] + # # all along x axis and calculation for y axis + # test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90, 0) + # # all along y axis and calculation for y axis + # test_xyz( + # np.zeros_like(x), + # x, + # np.zeros_like(x), + # angles_rot_azim_90, + # gamma_exp, + # ) + # # all along z axis and calculation for y axis + # test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90, 0) + + # # for elevation it is important to check, that IF elevation is 90° or 270° it does + # # not matter how we rotated before, since any rotation around z (in XY plane) + # # followed by a rotation around x' (in YZ' plane) by 90° will result in the same + # # coordinates, (when the structure is two points with zero extend) + + # # test with [0, 90] + # angles_rot_azim_90_elev_90 = [0, np.pi / 2] + # # all along x axis and calculation for z axis + # test_xyz( + # x, + # np.zeros_like(x), + # np.zeros_like(x), + # angles_rot_azim_90_elev_90, + # 0, + # ) + # # all along y axis and calculation for z axis + # test_xyz( + # np.zeros_like(x), + # x, + # np.zeros_like(x), + # angles_rot_azim_90_elev_90, + # 0, + # ) + # # all along z axis and calculation for z axis + # test_xyz( + # np.zeros_like(x), + # np.zeros_like(x), + # x, + # angles_rot_azim_90_elev_90, + # gamma_exp, + # ) + + # # test with [90, 90] + # angles_rot_azim_90_elev_90 = [np.pi / 2, np.pi / 2] + # # all along x axis and calculation for z axis + # test_xyz( + # x, + # np.zeros_like(x), + # np.zeros_like(x), + # angles_rot_azim_90_elev_90, + # 0, + # ) + # # all along y axis and calculation for z axis + # test_xyz( + # np.zeros_like(x), + # x, + # np.zeros_like(x), + # angles_rot_azim_90_elev_90, + # 0, + # ) + # # all along z axis and calculation for z axis + # test_xyz( + # np.zeros_like(x), + # np.zeros_like(x), + # x, + # angles_rot_azim_90_elev_90, + # gamma_exp, + # ) + + # def test_angles_3D_estim(self): + + # seed = gs.random.MasterRNG(19970221) + # rng = np.random.RandomState(seed()) + # rng = np.random + # x = rng.randint(0, 50, size=1000) + # y = rng.randint(0, 50, size=1000) + # z = rng.randint(0, 50, size=1000) + + # model = gs.Exponential( + # dim=3, var=1, len_scale=[10, 5, 2], angles=[np.pi / 8, np.pi / 16] + # ) + # model_maj = gs.Exponential(dim=1, var=1, len_scale=[10]) + # model_min1 = gs.Exponential(dim=1, var=1, len_scale=[5]) + # model_min2 = gs.Exponential(dim=1, var=1, len_scale=[2]) + + # srf = gs.SRF(model, seed=20170519) + # field = srf((x, y, z)) + + # bins = np.arange(0, 25, 5) + # angle_mask = [np.pi / 8, np.pi / 16] + # angle_tol = 22.5 + + # # in x' + # bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured( + # (x, y, z), + # field, + # bins, + # angles=[angle_mask[0], angle_mask[1]], + # angles_tol=np.deg2rad(angle_tol), + # ) + + # # in y' + # bin_centers_min1, gamma_min1 = gs.vario_estimate_unstructured( + # (x, y, z), + # field, + # bins, + # angles=[angle_mask[0] + 0.5 * np.pi, angle_mask[1]], + # angles_tol=np.deg2rad(angle_tol), + # ) + + # # in z' + # bin_centers_min2, gamma_min2 = gs.vario_estimate_unstructured( + # (x, y, z), + # field, + # bins, + # angles=[angle_mask[0], angle_mask[1] + 0.5 * np.pi], + # angles_tol=np.deg2rad(angle_tol), + # ) + + # gamma_maj_real = model_maj.variogram(bin_centers_maj) + # gamma_min1_real = model_min1.variogram(bin_centers_min1) + # gamma_min2_real = model_min2.variogram(bin_centers_min2) + + # # we have no real way of testing values, but we can test some basic properties which definitelly need to be true + # # test that the major estimate aligns better with major real than the minor reals + # self.assertTrue( + # np.sum((gamma_maj_real - gamma_maj) ** 2) + # < np.sum((gamma_min1_real - gamma_maj) ** 2) + # ) + # self.assertTrue( + # np.sum((gamma_maj_real - gamma_maj) ** 2) + # < np.sum((gamma_min2_real - gamma_maj) ** 2) + # ) + # # test that the minor1 estimate aligns better with minor1 real than major real and minor2 real + # self.assertTrue( + # np.sum((gamma_min1_real - gamma_min1) ** 2) + # < np.sum((gamma_maj_real - gamma_min1) ** 2) + # ) + # self.assertTrue( + # np.sum((gamma_min1_real - gamma_min1) ** 2) + # < np.sum((gamma_min2_real - gamma_min1) ** 2) + # ) + # # test that the minor2 estimate aligns better with minor2 real than major real and minor1 real + # self.assertTrue( + # np.sum((gamma_min2_real - gamma_min2) ** 2) + # < np.sum((gamma_maj_real - gamma_min2) ** 2) + # ) + # self.assertTrue( + # np.sum((gamma_min2_real - gamma_min2) ** 2) + # < np.sum((gamma_min1_real - gamma_min2) ** 2) + # ) + # # test that all variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance + # self.assertTrue( + # (np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1 + # ) + # self.assertTrue( + # (np.mean(gamma_min1[-5:]) - np.var(field)) / np.var(field) < 0.1 + # ) + # self.assertTrue( + # (np.mean(gamma_min2[-5:]) - np.var(field)) / np.var(field) < 0.1 + # ) if __name__ == "__main__": From b8fb4a001b1313340f38bc798a209f608a9dc748 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 15:15:35 +0100 Subject: [PATCH 07/44] tools: add ang2dir to init --- gstools/tools/__init__.py | 11 ++++++++++- gstools/tools/geometric.py | 3 ++- 2 files changed, 12 insertions(+), 2 deletions(-) diff --git a/gstools/tools/__init__.py b/gstools/tools/__init__.py index 33e4c5b5..368cc900 100644 --- a/gstools/tools/__init__.py +++ b/gstools/tools/__init__.py @@ -33,6 +33,7 @@ .. autosummary:: xyz2pos pos2xyz + ang2dir r3d_x r3d_y r3d_z @@ -59,7 +60,14 @@ tpl_gau_spec_dens, ) -from gstools.tools.geometric import r3d_x, r3d_y, r3d_z, xyz2pos, pos2xyz +from gstools.tools.geometric import ( + r3d_x, + r3d_y, + r3d_z, + xyz2pos, + pos2xyz, + ang2dir, +) __all__ = [ "vtk_export", @@ -77,6 +85,7 @@ "tpl_gau_spec_dens", "xyz2pos", "pos2xyz", + "ang2dir", "r3d_x", "r3d_y", "r3d_z", diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index 1faa5f56..dfd5dd03 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -12,12 +12,13 @@ r3d_z pos2xyz xyz2pos + ang2dir """ # pylint: disable=C0103 import numpy as np -__all__ = ["r3d_x", "r3d_y", "r3d_z", "pos2xyz", "xyz2pos"] +__all__ = ["r3d_x", "r3d_y", "r3d_z", "pos2xyz", "xyz2pos", "ang2dir"] # Geometric functions ######################################################### From c044b620857fbe0666d942ad02e3960261db54e6 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 15:19:51 +0100 Subject: [PATCH 08/44] vario: add multi field option; add nodata option --- gstools/variogram/estimator.pyx | 34 ++++++++++++++++++++------------- gstools/variogram/variogram.py | 30 ++++++++++++++++++++++------- 2 files changed, 44 insertions(+), 20 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 447926df..0753b210 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -10,7 +10,7 @@ import numpy as np cimport cython from cython.parallel import prange, parallel from libcpp.vector cimport vector -from libc.math cimport fabs, sqrt, acos, M_PI +from libc.math cimport fabs, sqrt, isnan, acos, M_PI cimport numpy as np @@ -164,7 +164,7 @@ cdef _normalization_func_vec choose_estimator_normalization_vec(str estimator_ty def directional( const int dim, - const double[:] f, + const double[:,:] f, const double[:] bin_edges, const double[:,:] pos, const double[:,:] direction, # should be normed @@ -172,9 +172,9 @@ def directional( const double bandwith=-1.0, # negative values to turn of bandwith search str estimator_type='m' ): - if pos.shape[1] != f.shape[0]: + if pos.shape[1] != f.shape[1]: raise ValueError('len(pos) = {0} != len(f) = {1} '. - format(pos.shape[1], f.shape[0])) + format(pos.shape[1], f.shape[1])) if bin_edges.shape[0] < 2: raise ValueError('len(bin_edges) too small') @@ -191,12 +191,13 @@ def directional( cdef int i_max = bin_edges.shape[0] - 1 cdef int j_max = pos.shape[1] - 1 cdef int k_max = pos.shape[1] + cdef int f_max = f.shape[0] cdef double[:,:] variogram = np.zeros((d_max, len(bin_edges)-1)) cdef long[:,:] counts = np.zeros((d_max, len(bin_edges)-1), dtype=long) cdef vector[double] pos1 = vector[double](dim, 0.0) cdef vector[double] pos2 = vector[double](dim, 0.0) - cdef int i, j, k, d + cdef int i, j, k, m, d cdef DTYPE_t dist for i in prange(i_max, nogil=True): @@ -206,22 +207,25 @@ def directional( if dist >= bin_edges[i] and dist < bin_edges[i+1]: for d in range(d_max): if dir_test(dim, pos, direction, angles_tol, bandwith, k, j, d): - counts[d, i] += 1 - variogram[d, i] += estimator_func(f[k] - f[j]) + for m in range(f_max): + # skip no data values + if not (isnan(f[m,k]) or isnan(f[m,j])): + counts[d, i] += 1 + variogram[d, i] += estimator_func(f[m,k] - f[m,j]) normalization_func_vec(variogram, counts) return np.asarray(variogram), np.asarray(counts) def unstructured( const int dim, - const double[:] f, + const double[:,:] f, const double[:] bin_edges, const double[:,:] pos, str estimator_type='m' ): - if pos.shape[1] != f.shape[0]: + if pos.shape[1] != f.shape[1]: raise ValueError('len(pos) = {0} != len(f) = {1} '. - format(pos.shape[1], f.shape[0])) + format(pos.shape[1], f.shape[1])) if bin_edges.shape[0] < 2: raise ValueError('len(bin_edges) too small') @@ -234,12 +238,13 @@ def unstructured( cdef int i_max = bin_edges.shape[0] - 1 cdef int j_max = pos.shape[1] - 1 cdef int k_max = pos.shape[1] + cdef int f_max = f.shape[0] cdef vector[double] variogram = vector[double](len(bin_edges)-1, 0.0) cdef vector[long] counts = vector[long](len(bin_edges)-1, 0) cdef vector[double] pos1 = vector[double](dim, 0.0) cdef vector[double] pos2 = vector[double](dim, 0.0) - cdef int i, j, k + cdef int i, j, k, m cdef DTYPE_t dist for i in prange(i_max, nogil=True): @@ -247,8 +252,11 @@ def unstructured( for k in range(j+1, k_max): dist = distance(dim, pos, j, k) if dist >= bin_edges[i] and dist < bin_edges[i+1]: - counts[i] += 1 - variogram[i] += estimator_func(f[k] - f[j]) + for m in range(f_max): + # skip no data values + if not (isnan(f[m,k]) or isnan(f[m,j])): + counts[i] += 1 + variogram[i] += estimator_func(f[m,k] - f[m,j]) normalization_func(variogram, counts) return np.asarray(variogram), np.asarray(counts) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 0b5d8229..63727660 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -49,6 +49,7 @@ def vario_estimate_unstructured( sampling_size=None, sampling_seed=None, estimator="matheron", + no_data=np.nan, return_counts=False, ): r""" @@ -79,8 +80,11 @@ def vario_estimate_unstructured( pos : :class:`list` the position tuple, containing main direction and transversal directions - field : :class:`numpy.ndarray` - the spatially distributed data + field : :class:`numpy.ndarray` or :class:`list` of :class:`numpy.ndarray` + The spatially distributed data. + You can pass a list of fields, that will be used simultaneously. + This could be helpful, when there are multiple realizations at the + same points, with the same statistical properties. bin_edges : :class:`numpy.ndarray` the bins on which the variogram will be calculated direction : :class:`list` of :class:`numpy.ndarray`, optional @@ -127,6 +131,9 @@ def vario_estimate_unstructured( * "cressie": an estimator more robust to outliers Default: "matheron" + no_data : :class:`float`, optional + Value to identify missing data in the given field. + Default: `np.nan` return_counts: class:`bool`, optional if set to true, this function will also return the number of data points found at each lag distance as a third return value @@ -144,11 +151,20 @@ def vario_estimate_unstructured( ----- Internally uses double precision and also returns doubles. """ - # TODO check_mesh - field = np.array(field, ndmin=1, dtype=np.double) + # allow multiple fields at same positions (ndmin=2: first axis -> field ID) + field = np.array(field, ndmin=2, dtype=np.double) bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double) x, y, z, dim = pos2xyz(pos, calc_dim=True, dtype=np.double) bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 + # check_mesh shape + if len(field.shape) > 2 or field.shape[1] != len(x): + try: + field = field.reshape((-1, len(x))) + except ValueError: + raise ValueError("'field' has wrong shape") + # set no_data values + if not np.isnan(no_data): + field[np.isclose(field, float(no_data))] = np.nan # initialize number of directions dir_no = 0 if direction is not None and dim > 1: @@ -175,11 +191,11 @@ def vario_estimate_unstructured( # prepare positions pos = np.array(xyz2pos(x, y, z, dtype=np.double, max_dim=dim)) # prepare sampled variogram - if sampling_size is not None and sampling_size < len(field): + if sampling_size is not None and sampling_size < len(x): sampled_idx = np.random.RandomState(sampling_seed).choice( - np.arange(len(field)), sampling_size, replace=False + np.arange(len(x)), sampling_size, replace=False ) - field = field[sampled_idx] + field = field[:, sampled_idx] pos = pos[:, sampled_idx] # select variogram estimator cython_estimator = _set_estimator(estimator) From a4b4e6b66349478f21f99fb5dd1752b6f72bc375 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 15:20:17 +0100 Subject: [PATCH 09/44] example: add multi field vario estimation example --- examples/03_variogram/03_multi_vario.py | 43 +++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100755 examples/03_variogram/03_multi_vario.py diff --git a/examples/03_variogram/03_multi_vario.py b/examples/03_variogram/03_multi_vario.py new file mode 100755 index 00000000..e66030b2 --- /dev/null +++ b/examples/03_variogram/03_multi_vario.py @@ -0,0 +1,43 @@ +""" +Multi-field variogram estimation +-------------------------------- + +In this example, we demonstrate how to estimate a variogram from multiple +fields at the same point-set, that should have the same statistical properties. +""" +import numpy as np +import gstools as gs +import matplotlib.pyplot as plt + + +x = np.random.RandomState(19970221).rand(1000) * 100.0 +y = np.random.RandomState(20011012).rand(1000) * 100.0 +model = gs.Exponential(dim=2, var=2, len_scale=8) +srf = gs.SRF(model, mean=0) + +############################################################################### +# Generate two synthetic fields with an exponential model. + +field1 = srf((x, y), seed=19970221) +field2 = srf((x, y), seed=20011012) +fields = [field1, field2] + +############################################################################### +# Now we estimate the variograms for both fields individually and then again +# simultaneously with only one call. + +bins = np.arange(40) +bin_center, gamma1 = gs.vario_estimate_unstructured((x, y), field1, bins) +bin_center, gamma2 = gs.vario_estimate_unstructured((x, y), field2, bins) +bin_center, gamma = gs.vario_estimate_unstructured((x, y), fields, bins) + +############################################################################### +# Now we demonstrate, that the mean variogram from both fields coincides +# with the joined estimated one. + +plt.plot(bin_center, gamma1, label="field 1") +plt.plot(bin_center, gamma2, label="field 2") +plt.plot(bin_center, gamma, label="joined fields") +plt.plot(bin_center, 0.5 * (gamma1 + gamma2), ":", label="field 1+2 mean") +plt.legend() +plt.show() From e3461cf3abcf8d587347e792a1b92d32019757a6 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 15:20:47 +0100 Subject: [PATCH 10/44] tests/vario: add test for multi field and nodata --- tests/test_variogram_unstructured.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 02fd30e5..30de4cba 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -250,6 +250,34 @@ def test_assertions(self): ValueError, vario_estimate_unstructured, [x], field_e, bins ) + def test_multi_field(self): + x = np.random.RandomState(19970221).rand(100) * 100.0 + model = gs.Exponential(dim=1, var=2, len_scale=10) + srf = gs.SRF(model) + field1 = srf(x, seed=19970221) + field2 = srf(x, seed=20011012) + bins = np.arange(20) * 2 + bin_center, gamma1 = gs.vario_estimate_unstructured(x, field1, bins) + bin_center, gamma2 = gs.vario_estimate_unstructured(x, field2, bins) + bin_center, gamma = gs.vario_estimate_unstructured( + x, [field1, field2], bins + ) + gamma_mean = 0.5 * (gamma1 + gamma2) + for i in range(len(gamma)): + self.assertAlmostEqual(gamma[i], gamma_mean[i], places=2) + + def test_no_data(self): + x1 = np.random.RandomState(19970221).rand(100) * 100.0 + field1 = np.random.RandomState(20011012).rand(100) * 100.0 + field1[:10] = np.nan + x2 = x1[10:] + field2 = field1[10:] + bins = np.arange(20) * 2 + bin_center, gamma1 = gs.vario_estimate_unstructured(x1, field1, bins) + bin_center, gamma2 = gs.vario_estimate_unstructured(x2, field2, bins) + for i in range(len(gamma1)): + self.assertAlmostEqual(gamma1[i], gamma2[i], places=2) + # def test_angles_2D_x2x(self): # x = self.test_data_rotation_1["x"] From 89e83bb63701ca4282414ea4e9f7ef8b89458110 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 15:54:40 +0100 Subject: [PATCH 11/44] vario: add posibility to give structured field; rename unstructured variant to vario_estimate --- gstools/__init__.py | 8 +++++++- gstools/variogram/__init__.py | 8 +++++++- gstools/variogram/variogram.py | 23 ++++++++++++++++++++--- 3 files changed, 34 insertions(+), 5 deletions(-) diff --git a/gstools/__init__.py b/gstools/__init__.py index 1bafdfe4..ca0b6950 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -95,6 +95,7 @@ .. currentmodule:: gstools.variogram .. autosummary:: + vario_estimate vario_estimate_structured vario_estimate_unstructured """ @@ -110,6 +111,7 @@ to_vtk_unstructured, ) from gstools.variogram import ( + vario_estimate, vario_estimate_structured, vario_estimate_unstructured, ) @@ -154,7 +156,11 @@ "TPLStable", ] -__all__ += ["vario_estimate_structured", "vario_estimate_unstructured"] +__all__ += [ + "vario_estimate", + "vario_estimate_structured", + "vario_estimate_unstructured", +] __all__ += [ "SRF", diff --git a/gstools/variogram/__init__.py b/gstools/variogram/__init__.py index c5735c52..7f496713 100644 --- a/gstools/variogram/__init__.py +++ b/gstools/variogram/__init__.py @@ -8,6 +8,7 @@ ^^^^^^^^^^^^^^^^^^^^ .. autosummary:: + vario_estimate vario_estimate_unstructured vario_estimate_structured @@ -15,8 +16,13 @@ """ from gstools.variogram.variogram import ( + vario_estimate, vario_estimate_structured, vario_estimate_unstructured, ) -__all__ = ["vario_estimate_unstructured", "vario_estimate_structured"] +__all__ = [ + "vario_estimate", + "vario_estimate_unstructured", + "vario_estimate_structured", +] diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 63727660..e862cab5 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -7,6 +7,7 @@ The following functions are provided .. autosummary:: + vario_estimate vario_estimate_unstructured vario_estimate_structured """ @@ -14,6 +15,7 @@ import numpy as np +from gstools.field.tools import reshape_axis_from_struct_to_unstruct from gstools.tools.geometric import pos2xyz, xyz2pos, ang2dir from gstools.variogram.estimator import ( unstructured, @@ -22,7 +24,11 @@ directional, ) -__all__ = ["vario_estimate_unstructured", "vario_estimate_structured"] +__all__ = [ + "vario_estimate", + "vario_estimate_unstructured", + "vario_estimate_structured", +] def _set_estimator(estimator): @@ -38,7 +44,7 @@ def _set_estimator(estimator): return cython_estimator -def vario_estimate_unstructured( +def vario_estimate( pos, field, bin_edges, @@ -49,11 +55,12 @@ def vario_estimate_unstructured( sampling_size=None, sampling_seed=None, estimator="matheron", + mesh_type="unstructured", no_data=np.nan, return_counts=False, ): r""" - Estimates the variogram on a unstructured grid. + Estimates the empirical variogram. The algorithm calculates following equation: @@ -131,6 +138,10 @@ def vario_estimate_unstructured( * "cressie": an estimator more robust to outliers Default: "matheron" + mesh_type : :class:`str`, optional + 'structured' / 'unstructured', indicates whether the pos tuple + describes the axis or the point coordinates. + Default: `'unstructured'` no_data : :class:`float`, optional Value to identify missing data in the given field. Default: `np.nan` @@ -157,6 +168,8 @@ def vario_estimate_unstructured( x, y, z, dim = pos2xyz(pos, calc_dim=True, dtype=np.double) bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 # check_mesh shape + if mesh_type != "unstructured": + x, y, z, __ = reshape_axis_from_struct_to_unstruct(dim, x, y, z) if len(field.shape) > 2 or field.shape[1] != len(x): try: field = field.reshape((-1, len(x))) @@ -313,3 +326,7 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"): else: gamma = ma_structured(field, mask, cython_estimator) return gamma + + +# for backward compatibility +vario_estimate_unstructured = vario_estimate From c3a1eabbdc382cc68eb58dc1637a2352fcd6d12f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 16:11:37 +0100 Subject: [PATCH 12/44] vario: fix 'bandwidth' typo; add hints in documentation --- gstools/variogram/estimator.pyx | 12 ++++++------ gstools/variogram/variogram.py | 26 ++++++++++++++++++-------- 2 files changed, 24 insertions(+), 14 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 0753b210..34545664 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -36,7 +36,7 @@ cdef inline bint dir_test( const double[:,:] pos, const double[:,:] direction, const double angles_tol, - const double bandwith, + const double bandwidth, const int i, const int j, const int d @@ -48,18 +48,18 @@ cdef inline bint dir_test( cdef bint in_band = True cdef bint in_angle - # scalar-product calculation for bandwith projection and angle calculation + # scalar-product calculation for bandwidth projection and angle calculation for k in range(dim): s_prod += (pos[k,i] - pos[k,j]) * direction[d,k] p_norm += (pos[k,i] - pos[k,j])**2 p_norm = sqrt(p_norm) # calculate band-distance by projection of point-pair-vec to direction line - if bandwith > 0.0: + if bandwidth > 0.0: for k in range(dim): b_dist += ((pos[k,i] - pos[k,j]) - s_prod * direction[d,k]) ** 2 b_dist = sqrt(b_dist) - in_band = b_dist < bandwith + in_band = b_dist < bandwidth # use smallest angle by taking absolut value for arccos angle formula in_angle = acos(fabs(s_prod) / p_norm) < angles_tol @@ -169,7 +169,7 @@ def directional( const double[:,:] pos, const double[:,:] direction, # should be normed const double angles_tol=M_PI/8.0, - const double bandwith=-1.0, # negative values to turn of bandwith search + const double bandwidth=-1.0, # negative values to turn of bandwidth search str estimator_type='m' ): if pos.shape[1] != f.shape[1]: @@ -206,7 +206,7 @@ def directional( dist = distance(dim, pos, j, k) if dist >= bin_edges[i] and dist < bin_edges[i+1]: for d in range(d_max): - if dir_test(dim, pos, direction, angles_tol, bandwith, k, j, d): + if dir_test(dim, pos, direction, angles_tol, bandwidth, k, j, d): for m in range(f_max): # skip no data values if not (isnan(f[m,k]) or isnan(f[m,j])): diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index e862cab5..9cd44890 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -51,7 +51,7 @@ def vario_estimate( direction=None, angles=None, angles_tol=np.pi / 8, - bandwith=None, + bandwidth=None, sampling_size=None, sampling_seed=None, estimator="matheron", @@ -82,6 +82,16 @@ def vario_estimate( being the bins. The Cressie estimator is more robust to outliers. + By provding `direction` vector[s] or angles, a directional variogram + can be calculated. If multiple directions are given, a set of variograms + will be returned. + Directional bining is controled by a given angle tolerance (`angles_tol`) + and an optional `bandwidth`, that truncates the width of the search band + around the given direction[s]. + + To reduce the calcuation time, `sampling_size` could be passed to sample + down the number of field points. + Parameters ---------- pos : :class:`list` @@ -97,8 +107,8 @@ def vario_estimate( direction : :class:`list` of :class:`numpy.ndarray`, optional directions to evaluate a directional variogram. Anglular tolerance is given by `angles_tol`. - Bandwith to cut off how wide the search for point pairs should be - is given by `bandwith`. + bandwidth to cut off how wide the search for point pairs should be + is given by `bandwidth`. You can provide multiple directions at once to get one variogram for each direction. For a single direction you can also use the `angles` parameter, @@ -118,8 +128,8 @@ def vario_estimate( within this direction from another point (the angular tolerance around the directional vector given by angles) Default: `np.pi/8` = 22.5° - bandwith : class:`float`, optional - Bandwith to cut off the angular tolerance for directional variograms. + bandwidth : class:`float`, optional + bandwidth to cut off the angular tolerance for directional variograms. If None is given, only the `angles_tol` parameter will control the point selection. Default: :any:`None` @@ -198,8 +208,8 @@ def vario_estimate( raise ValueError("Zero length direction {}".format(direction)) # only unit-vectors for directions direction = np.divide(direction, norms[:, np.newaxis]) - # negative bandwith to turn it off - bandwith = float(bandwith) if bandwith is not None else -1.0 + # negative bandwidth to turn it off + bandwidth = float(bandwidth) if bandwidth is not None else -1.0 angles_tol = float(angles_tol) # prepare positions pos = np.array(xyz2pos(x, y, z, dtype=np.double, max_dim=dim)) @@ -229,7 +239,7 @@ def vario_estimate( pos, direction, angles_tol, - bandwith, + bandwidth, estimator_type=cython_estimator, ) if dir_no == 1: From 8f9c436fe44f651768602de0fa31ef62a9ba828c Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 16:57:53 +0100 Subject: [PATCH 13/44] vario: reorder input parameters for backward compatibility --- gstools/variogram/variogram.py | 46 +++++++++++++++++----------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 9cd44890..8e54ef84 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -48,15 +48,15 @@ def vario_estimate( pos, field, bin_edges, + sampling_size=None, + sampling_seed=None, + estimator="matheron", direction=None, angles=None, angles_tol=np.pi / 8, bandwidth=None, - sampling_size=None, - sampling_seed=None, - estimator="matheron", - mesh_type="unstructured", no_data=np.nan, + mesh_type="unstructured", return_counts=False, ): r""" @@ -104,6 +104,21 @@ def vario_estimate( same points, with the same statistical properties. bin_edges : :class:`numpy.ndarray` the bins on which the variogram will be calculated + sampling_size : :class:`int` or :any:`None`, optional + for large input data, this method can take a long + time to compute the variogram, therefore this argument specifies + the number of data points to sample randomly + Default: :any:`None` + sampling_seed : :class:`int` or :any:`None`, optional + seed for samples if sampling_size is given. + Default: :any:`None` + estimator : :class:`str`, optional + the estimator function, possible choices: + + * "matheron": the standard method of moments of Matheron + * "cressie": an estimator more robust to outliers + + Default: "matheron" direction : :class:`list` of :class:`numpy.ndarray`, optional directions to evaluate a directional variogram. Anglular tolerance is given by `angles_tol`. @@ -133,28 +148,13 @@ def vario_estimate( If None is given, only the `angles_tol` parameter will control the point selection. Default: :any:`None` - sampling_size : :class:`int` or :any:`None`, optional - for large input data, this method can take a long - time to compute the variogram, therefore this argument specifies - the number of data points to sample randomly - Default: :any:`None` - sampling_seed : :class:`int` or :any:`None`, optional - seed for samples if sampling_size is given. - Default: :any:`None` - estimator : :class:`str`, optional - the estimator function, possible choices: - - * "matheron": the standard method of moments of Matheron - * "cressie": an estimator more robust to outliers - - Default: "matheron" + no_data : :class:`float`, optional + Value to identify missing data in the given field. + Default: `np.nan` mesh_type : :class:`str`, optional 'structured' / 'unstructured', indicates whether the pos tuple describes the axis or the point coordinates. Default: `'unstructured'` - no_data : :class:`float`, optional - Value to identify missing data in the given field. - Default: `np.nan` return_counts: class:`bool`, optional if set to true, this function will also return the number of data points found at each lag distance as a third return value @@ -188,7 +188,7 @@ def vario_estimate( # set no_data values if not np.isnan(no_data): field[np.isclose(field, float(no_data))] = np.nan - # initialize number of directions + # set directions dir_no = 0 if direction is not None and dim > 1: direction = np.array(direction, ndmin=2, dtype=np.double) From eac5e87b20c63bfd3ca46f45aa81aab4c2f62b36 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 17:02:05 +0100 Subject: [PATCH 14/44] vario: angles doc fix --- gstools/variogram/variogram.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 8e54ef84..3eb185f4 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -136,7 +136,7 @@ def vario_estimate( for 2d supply one angle which is azimuth φ (ccw from +x in xy plane) for 3d supply two angles which are azimuth φ (ccw from +x in xy plane) and inclination θ (cw from +z). - Can be used instead of direction for a single main axis. + Can be used instead of direction. Default: :any:`None` angles_tol : class:`float`, optional the tolerance around the variogram angle to count a point as being From 92d888a1b4f3619e7cfe6b16e18054fd0f0718cb Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 18:52:29 +0100 Subject: [PATCH 15/44] tests: cleanup --- tests/test_variogram_unstructured.py | 381 ++------------------------- 1 file changed, 27 insertions(+), 354 deletions(-) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 30de4cba..13c7cff3 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -5,7 +5,7 @@ import unittest import numpy as np -from gstools import vario_estimate_unstructured, Exponential, SRF +from gstools import vario_estimate, Exponential, SRF import gstools as gs @@ -20,7 +20,7 @@ def setUp(self): field = srf((x, y)) bins = np.arange(10) - bin_center, gamma = gs.vario_estimate_unstructured((x,), field, bins) + bin_center, gamma = gs.vario_estimate((x,), field, bins) idx = np.argsort(x) self.test_data_rotation_1 = { "gamma": gamma, @@ -52,21 +52,21 @@ def test_doubles(self): dtype=np.double, ) bins = np.arange(1, 11, 1, dtype=np.double) - bin_centres, gamma = vario_estimate_unstructured([x], z, bins) + bin_centres, gamma = vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 0.4917, places=4) def test_ints(self): x = np.arange(1, 5, 1, dtype=int) z = np.array((10, 20, 30, 40), dtype=int) bins = np.arange(1, 11, 1, dtype=int) - bin_centres, gamma = vario_estimate_unstructured([x], z, bins) + bin_centres, gamma = vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 50.0, places=4) def test_np_int(self): x = np.arange(1, 5, 1, dtype=np.int) z = np.array((10, 20, 30, 40), dtype=np.int) bins = np.arange(1, 11, 1, dtype=np.int) - bin_centres, gamma = vario_estimate_unstructured([x], z, bins) + bin_centres, gamma = vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 50.0, places=4) def test_mixed(self): @@ -76,26 +76,26 @@ def test_mixed(self): dtype=np.double, ) bins = np.arange(1, 11, 1, dtype=int) - bin_centres, gamma = vario_estimate_unstructured([x], z, bins) + bin_centres, gamma = vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 0.4917, places=4) x = np.arange(1, 5, 1, dtype=np.double) z = np.array((10, 20, 30, 40), dtype=int) bins = np.arange(1, 11, 1, dtype=int) - bin_centres, gamma = vario_estimate_unstructured([x], z, bins) + bin_centres, gamma = vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 50.0, places=4) x = np.arange(1, 5, 1, dtype=np.double) z = np.array((10, 20, 30, 40), dtype=int) bins = np.arange(1, 11, 1, dtype=np.double) - bin_centres, gamma = vario_estimate_unstructured([x], z, bins) + bin_centres, gamma = vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 50.0, places=4) def test_list(self): x = np.arange(1, 11, 1, dtype=np.double) z = [41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3] bins = np.arange(1, 11, 1, dtype=np.double) - bin_centres, gamma = vario_estimate_unstructured([x], z, bins) + bin_centres, gamma = vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[1], 0.7625, places=4) def test_1d(self): @@ -106,7 +106,7 @@ def test_1d(self): dtype=np.double, ) bins = np.arange(1, 11, 1, dtype=np.double) - bin_centres, gamma = vario_estimate_unstructured([x], z, bins) + bin_centres, gamma = vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 0.4917, places=4) self.assertAlmostEqual(gamma[1], 0.7625, places=4) @@ -122,7 +122,7 @@ def test_uncorrelated_2d(self): bins = np.arange(0, 100, 10) - bin_centres, gamma = vario_estimate_unstructured((x, y), field, bins) + bin_centres, gamma = vario_estimate((x, y), field, bins) var = 1.0 / 12.0 self.assertAlmostEqual(gamma[0], var, places=2) @@ -143,7 +143,7 @@ def test_uncorrelated_3d(self): bins = np.arange(0, 100, 10) - bin_centres, gamma = vario_estimate_unstructured( + bin_centres, gamma = vario_estimate( (x, y, z), field, bins ) @@ -160,7 +160,7 @@ def test_sampling_1d(self): bins = np.arange(0, 100, 10) - bin_centres, gamma = vario_estimate_unstructured( + bin_centres, gamma = vario_estimate( [x], field, bins, sampling_size=5000, sampling_seed=1479373475 ) @@ -181,7 +181,7 @@ def test_sampling_2d(self): bins = np.arange(0, 100, 10) - bin_centres, gamma = vario_estimate_unstructured( + bin_centres, gamma = vario_estimate( (x, y), field, bins, sampling_size=2000, sampling_seed=1479373475 ) @@ -204,7 +204,7 @@ def test_sampling_3d(self): bins = np.arange(0, 100, 10) - bin_centres, gamma = vario_estimate_unstructured( + bin_centres, gamma = vario_estimate( (x, y, z), field, bins, @@ -229,25 +229,25 @@ def test_assertions(self): field_e = np.arange(0, 9) self.assertRaises( - ValueError, vario_estimate_unstructured, [x_e], field, bins + ValueError, vario_estimate, [x_e], field, bins ) self.assertRaises( - ValueError, vario_estimate_unstructured, (x, y_e), field, bins + ValueError, vario_estimate, (x, y_e), field, bins ) self.assertRaises( - ValueError, vario_estimate_unstructured, (x, y_e, z), field, bins + ValueError, vario_estimate, (x, y_e, z), field, bins ) self.assertRaises( - ValueError, vario_estimate_unstructured, (x, y, z_e), field, bins + ValueError, vario_estimate, (x, y, z_e), field, bins ) self.assertRaises( - ValueError, vario_estimate_unstructured, (x_e, y, z), field, bins + ValueError, vario_estimate, (x_e, y, z), field, bins ) self.assertRaises( - ValueError, vario_estimate_unstructured, (x, y, z), field_e, bins + ValueError, vario_estimate, (x, y, z), field_e, bins ) self.assertRaises( - ValueError, vario_estimate_unstructured, [x], field_e, bins + ValueError, vario_estimate, [x], field_e, bins ) def test_multi_field(self): @@ -257,9 +257,9 @@ def test_multi_field(self): field1 = srf(x, seed=19970221) field2 = srf(x, seed=20011012) bins = np.arange(20) * 2 - bin_center, gamma1 = gs.vario_estimate_unstructured(x, field1, bins) - bin_center, gamma2 = gs.vario_estimate_unstructured(x, field2, bins) - bin_center, gamma = gs.vario_estimate_unstructured( + bin_center, gamma1 = gs.vario_estimate(x, field1, bins) + bin_center, gamma2 = gs.vario_estimate(x, field2, bins) + bin_center, gamma = gs.vario_estimate( x, [field1, field2], bins ) gamma_mean = 0.5 * (gamma1 + gamma2) @@ -273,338 +273,11 @@ def test_no_data(self): x2 = x1[10:] field2 = field1[10:] bins = np.arange(20) * 2 - bin_center, gamma1 = gs.vario_estimate_unstructured(x1, field1, bins) - bin_center, gamma2 = gs.vario_estimate_unstructured(x2, field2, bins) + bin_center, gamma1 = gs.vario_estimate(x1, field1, bins) + bin_center, gamma2 = gs.vario_estimate(x2, field2, bins) for i in range(len(gamma1)): self.assertAlmostEqual(gamma1[i], gamma2[i], places=2) - # def test_angles_2D_x2x(self): - - # x = self.test_data_rotation_1["x"] - # field = self.test_data_rotation_1["field"] - # gamma_exp = self.test_data_rotation_1["gamma"] - # bins = self.test_data_rotation_1["bins"] - # y = np.zeros_like(x) - - # # test case 1.) - # # all along x axis on x axis - - # bin_centres, gamma = vario_estimate_unstructured( - # (x, y), field, bins, angles=[0] - # ) - - # for i in range(gamma.size): - # self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) - - # def test_angles_2D_y2x(self): - - # x = self.test_data_rotation_1["x"] - # field = self.test_data_rotation_1["field"] - # gamma_exp = self.test_data_rotation_1["gamma"] - # bins = self.test_data_rotation_1["bins"] - # y = np.zeros_like(x) - - # # test case 2.) - # # all along y axis on y axis but calculation for x axis - - # bin_centres, gamma = vario_estimate_unstructured( - # (y, x), field, bins, angles=[0] - # ) - - # for i in range(gamma.size): - # self.assertAlmostEqual(0, gamma[i], places=3) - - # def test_angles_2D_y2y(self): - - # x = self.test_data_rotation_1["x"] - # field = self.test_data_rotation_1["field"] - # gamma_exp = self.test_data_rotation_1["gamma"] - # bins = self.test_data_rotation_1["bins"] - # y = np.zeros_like(x) - - # # test case 3.) - # # all along y axis on y axis and calculation for y axis - - # bin_centres, gamma = vario_estimate_unstructured( - # (y, x), field, bins, angles=[np.pi / 2.0] - # ) - - # for i in range(gamma.size): - # self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) - - # def test_angles_2D_xy2x(self): - - # x = self.test_data_rotation_1["x"] - # field = self.test_data_rotation_1["field"] - # gamma_exp = self.test_data_rotation_1["gamma"] - # bins = self.test_data_rotation_1["bins"] - # y = np.zeros_like(x) - - # # test case 4.) - # # data along 45deg axis but calculation for x axis - - # ccos, csin = np.cos(np.pi / 4.0), np.sin(np.pi / 4.0) - - # xr = [xx * ccos - yy * csin for xx, yy in zip(x, y)] - # yr = [xx * csin + yy * ccos for xx, yy in zip(x, y)] - - # bin_centres, gamma = vario_estimate_unstructured( - # (xr, yr), field, bins, angles=[0] - # ) - - # for i in range(gamma.size): - # self.assertAlmostEqual(0, gamma[i], places=3) - - # def test_angles_2D_estim(self): - - # seed = gs.random.MasterRNG(19970221) - # rng = np.random.RandomState(seed()) - # rng = np.random - # x = rng.randint(0, 100, size=3000) - # y = rng.randint(0, 100, size=3000) - - # model = gs.Exponential( - # dim=2, var=1, len_scale=[12, 3], angles=np.pi / 8 - # ) - # model_maj = gs.Exponential(dim=1, var=1, len_scale=[12]) - # model_min = gs.Exponential(dim=1, var=1, len_scale=[3]) - - # srf = gs.SRF(model, seed=20170519) - # field = srf((x, y)) - - # bins = np.arange(0, 50, 2.5) - # angle_mask = 22.5 - # angle_tol = 22.5 - - # bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured( - # (x, y), - # field, - # bins, - # angles=[np.deg2rad(angle_mask)], - # angles_tol=np.deg2rad(angle_tol), - # ) - - # bin_centers_min, gamma_min = gs.vario_estimate_unstructured( - # (x, y), - # field, - # bins, - # angles=[np.deg2rad(angle_mask + 90.0)], - # angles_tol=np.deg2rad(angle_tol), - # ) - - # gamma_maj_real = model_maj.variogram(bin_centers_maj) - # gamma_min_real = model_min.variogram(bin_centers_min) - - # # we have no real way of testing values, but we can test some basic properties which definitelly need to be true - # # test that the major estimate aligns better with major real than minor real - # self.assertTrue( - # np.sum((gamma_maj_real - gamma_maj) ** 2) - # < np.sum((gamma_min_real - gamma_maj) ** 2) - # ) - # # test that the minor estimate aligns better with minor real than major real - # self.assertTrue( - # np.sum((gamma_min_real - gamma_min) ** 2) - # < np.sum((gamma_maj_real - gamma_min) ** 2) - # ) - # # test that both variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance - # self.assertTrue( - # (np.mean(gamma_min[-5:]) - np.var(field)) / np.var(field) < 0.1 - # ) - # self.assertTrue( - # (np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1 - # ) - - # def test_angles_line_3D(self): - - # x = self.test_data_rotation_1["x"] - # field = self.test_data_rotation_1["field"] - # gamma_exp = self.test_data_rotation_1["gamma"] - # bins = self.test_data_rotation_1["bins"] - - # def test_xyz(x, y, z, angles, gamma_exp): - # bin_centres, gamma = vario_estimate_unstructured( - # (x, y, z), field, bins, angles=angles - # ) - - # if np.ndim(gamma_exp) == 0: - # gamma_exp = np.ones_like(gamma) * gamma_exp - - # for i in range(gamma.size): - # self.assertAlmostEqual(gamma_exp[i], gamma[i], places=3) - - # # all along x axis and calculation for x axis - # test_xyz(x, np.zeros_like(x), np.zeros_like(x), [0], gamma_exp) - # # all along y axis and calculation for x axis - # test_xyz(np.zeros_like(x), x, np.zeros_like(x), [0], 0) - # # all along z axis and calculation for x axis - # test_xyz(np.zeros_like(x), np.zeros_like(x), x, [0], 0) - - # angles_rot_azim_90 = [np.pi / 2] - # # all along x axis and calculation for y axis - # test_xyz(x, np.zeros_like(x), np.zeros_like(x), angles_rot_azim_90, 0) - # # all along y axis and calculation for y axis - # test_xyz( - # np.zeros_like(x), - # x, - # np.zeros_like(x), - # angles_rot_azim_90, - # gamma_exp, - # ) - # # all along z axis and calculation for y axis - # test_xyz(np.zeros_like(x), np.zeros_like(x), x, angles_rot_azim_90, 0) - - # # for elevation it is important to check, that IF elevation is 90° or 270° it does - # # not matter how we rotated before, since any rotation around z (in XY plane) - # # followed by a rotation around x' (in YZ' plane) by 90° will result in the same - # # coordinates, (when the structure is two points with zero extend) - - # # test with [0, 90] - # angles_rot_azim_90_elev_90 = [0, np.pi / 2] - # # all along x axis and calculation for z axis - # test_xyz( - # x, - # np.zeros_like(x), - # np.zeros_like(x), - # angles_rot_azim_90_elev_90, - # 0, - # ) - # # all along y axis and calculation for z axis - # test_xyz( - # np.zeros_like(x), - # x, - # np.zeros_like(x), - # angles_rot_azim_90_elev_90, - # 0, - # ) - # # all along z axis and calculation for z axis - # test_xyz( - # np.zeros_like(x), - # np.zeros_like(x), - # x, - # angles_rot_azim_90_elev_90, - # gamma_exp, - # ) - - # # test with [90, 90] - # angles_rot_azim_90_elev_90 = [np.pi / 2, np.pi / 2] - # # all along x axis and calculation for z axis - # test_xyz( - # x, - # np.zeros_like(x), - # np.zeros_like(x), - # angles_rot_azim_90_elev_90, - # 0, - # ) - # # all along y axis and calculation for z axis - # test_xyz( - # np.zeros_like(x), - # x, - # np.zeros_like(x), - # angles_rot_azim_90_elev_90, - # 0, - # ) - # # all along z axis and calculation for z axis - # test_xyz( - # np.zeros_like(x), - # np.zeros_like(x), - # x, - # angles_rot_azim_90_elev_90, - # gamma_exp, - # ) - - # def test_angles_3D_estim(self): - - # seed = gs.random.MasterRNG(19970221) - # rng = np.random.RandomState(seed()) - # rng = np.random - # x = rng.randint(0, 50, size=1000) - # y = rng.randint(0, 50, size=1000) - # z = rng.randint(0, 50, size=1000) - - # model = gs.Exponential( - # dim=3, var=1, len_scale=[10, 5, 2], angles=[np.pi / 8, np.pi / 16] - # ) - # model_maj = gs.Exponential(dim=1, var=1, len_scale=[10]) - # model_min1 = gs.Exponential(dim=1, var=1, len_scale=[5]) - # model_min2 = gs.Exponential(dim=1, var=1, len_scale=[2]) - - # srf = gs.SRF(model, seed=20170519) - # field = srf((x, y, z)) - - # bins = np.arange(0, 25, 5) - # angle_mask = [np.pi / 8, np.pi / 16] - # angle_tol = 22.5 - - # # in x' - # bin_centers_maj, gamma_maj = gs.vario_estimate_unstructured( - # (x, y, z), - # field, - # bins, - # angles=[angle_mask[0], angle_mask[1]], - # angles_tol=np.deg2rad(angle_tol), - # ) - - # # in y' - # bin_centers_min1, gamma_min1 = gs.vario_estimate_unstructured( - # (x, y, z), - # field, - # bins, - # angles=[angle_mask[0] + 0.5 * np.pi, angle_mask[1]], - # angles_tol=np.deg2rad(angle_tol), - # ) - - # # in z' - # bin_centers_min2, gamma_min2 = gs.vario_estimate_unstructured( - # (x, y, z), - # field, - # bins, - # angles=[angle_mask[0], angle_mask[1] + 0.5 * np.pi], - # angles_tol=np.deg2rad(angle_tol), - # ) - - # gamma_maj_real = model_maj.variogram(bin_centers_maj) - # gamma_min1_real = model_min1.variogram(bin_centers_min1) - # gamma_min2_real = model_min2.variogram(bin_centers_min2) - - # # we have no real way of testing values, but we can test some basic properties which definitelly need to be true - # # test that the major estimate aligns better with major real than the minor reals - # self.assertTrue( - # np.sum((gamma_maj_real - gamma_maj) ** 2) - # < np.sum((gamma_min1_real - gamma_maj) ** 2) - # ) - # self.assertTrue( - # np.sum((gamma_maj_real - gamma_maj) ** 2) - # < np.sum((gamma_min2_real - gamma_maj) ** 2) - # ) - # # test that the minor1 estimate aligns better with minor1 real than major real and minor2 real - # self.assertTrue( - # np.sum((gamma_min1_real - gamma_min1) ** 2) - # < np.sum((gamma_maj_real - gamma_min1) ** 2) - # ) - # self.assertTrue( - # np.sum((gamma_min1_real - gamma_min1) ** 2) - # < np.sum((gamma_min2_real - gamma_min1) ** 2) - # ) - # # test that the minor2 estimate aligns better with minor2 real than major real and minor1 real - # self.assertTrue( - # np.sum((gamma_min2_real - gamma_min2) ** 2) - # < np.sum((gamma_maj_real - gamma_min2) ** 2) - # ) - # self.assertTrue( - # np.sum((gamma_min2_real - gamma_min2) ** 2) - # < np.sum((gamma_min1_real - gamma_min2) ** 2) - # ) - # # test that all variograms converge within reasonable closeness (less than 10% rel error) to the actual field variance - # self.assertTrue( - # (np.mean(gamma_maj[-5:]) - np.var(field)) / np.var(field) < 0.1 - # ) - # self.assertTrue( - # (np.mean(gamma_min1[-5:]) - np.var(field)) / np.var(field) < 0.1 - # ) - # self.assertTrue( - # (np.mean(gamma_min2[-5:]) - np.var(field)) / np.var(field) < 0.1 - # ) - if __name__ == "__main__": unittest.main() From 66bb632efa7d76632bf06932677266083b34f95f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 7 Nov 2020 18:53:31 +0100 Subject: [PATCH 16/44] doc: use 'vario_estimate' in all examples --- README.md | 2 +- docs/source/index.rst | 2 +- examples/03_variogram/00_fit_variogram.py | 2 +- examples/03_variogram/01_variogram_estimation.py | 2 +- examples/03_variogram/02_find_best_model.py | 2 +- examples/03_variogram/03_multi_vario.py | 6 +++--- 6 files changed, 8 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 81182352..15de8d9f 100644 --- a/README.md +++ b/README.md @@ -154,7 +154,7 @@ srf = gs.SRF(model, mean=0, seed=19970221) field = srf((x, y)) # estimate the variogram of the field with 40 bins bins = np.arange(40) -bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins) +bin_center, gamma = gs.vario_estimate((x, y), field, bins) # fit the variogram with a stable model. (no nugget fitted) fit_model = gs.Stable(dim=2) fit_model.fit_variogram(bin_center, gamma, nugget=False) diff --git a/docs/source/index.rst b/docs/source/index.rst index bc3d23e9..f76f299f 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -182,7 +182,7 @@ model again. field = srf((x, y)) # estimate the variogram of the field with 40 bins bins = np.arange(40) - bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins) + bin_center, gamma = gs.vario_estimate((x, y), field, bins) # fit the variogram with a stable model. (no nugget fitted) fit_model = gs.Stable(dim=2) fit_model.fit_variogram(bin_center, gamma, nugget=False) diff --git a/examples/03_variogram/00_fit_variogram.py b/examples/03_variogram/00_fit_variogram.py index 8564209f..c1ddae62 100644 --- a/examples/03_variogram/00_fit_variogram.py +++ b/examples/03_variogram/00_fit_variogram.py @@ -18,7 +18,7 @@ # Estimate the variogram of the field with 40 bins. bins = np.arange(40) -bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins) +bin_center, gamma = gs.vario_estimate((x, y), field, bins) ############################################################################### # Fit the variogram with a stable model (no nugget fitted). diff --git a/examples/03_variogram/01_variogram_estimation.py b/examples/03_variogram/01_variogram_estimation.py index 08619386..f0574adf 100644 --- a/examples/03_variogram/01_variogram_estimation.py +++ b/examples/03_variogram/01_variogram_estimation.py @@ -146,7 +146,7 @@ def generate_transmissivity(): bins = np.linspace(0, 10, 50) -bin_center, gamma = gs.vario_estimate_unstructured( +bin_center, gamma = gs.vario_estimate( (x_u, y_u), herten_log_trans.reshape(-1), bins, diff --git a/examples/03_variogram/02_find_best_model.py b/examples/03_variogram/02_find_best_model.py index dc5d0f0b..1e7b9f87 100755 --- a/examples/03_variogram/02_find_best_model.py +++ b/examples/03_variogram/02_find_best_model.py @@ -19,7 +19,7 @@ # Estimate the variogram of the field with 40 bins and plot the result. bins = np.arange(40) -bin_center, gamma = gs.vario_estimate_unstructured((x, y), field, bins) +bin_center, gamma = gs.vario_estimate((x, y), field, bins) ############################################################################### # Define a set of models to test. diff --git a/examples/03_variogram/03_multi_vario.py b/examples/03_variogram/03_multi_vario.py index e66030b2..08951c08 100755 --- a/examples/03_variogram/03_multi_vario.py +++ b/examples/03_variogram/03_multi_vario.py @@ -27,9 +27,9 @@ # simultaneously with only one call. bins = np.arange(40) -bin_center, gamma1 = gs.vario_estimate_unstructured((x, y), field1, bins) -bin_center, gamma2 = gs.vario_estimate_unstructured((x, y), field2, bins) -bin_center, gamma = gs.vario_estimate_unstructured((x, y), fields, bins) +bin_center, gamma1 = gs.vario_estimate((x, y), field1, bins) +bin_center, gamma2 = gs.vario_estimate((x, y), field2, bins) +bin_center, gamma = gs.vario_estimate((x, y), fields, bins) ############################################################################### # Now we demonstrate, that the mean variogram from both fields coincides From e1cc3c5ad2bb3e02ba65416aa43c3f9a0ae061d0 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 8 Nov 2020 13:01:13 +0100 Subject: [PATCH 17/44] vario: cleanup mask handling in structured vario-estimate --- gstools/variogram/variogram.py | 42 ++++++++++++++-------------------- 1 file changed, 17 insertions(+), 25 deletions(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 3eb185f4..9f935af6 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -31,6 +31,9 @@ ] +AXIS_DIR = {"x": 0, "y": 1, "z": 2} + + def _set_estimator(estimator): """Translate the verbose Python estimator identifier to single char.""" if estimator.lower() == "matheron": @@ -300,42 +303,31 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"): ----- Internally uses double precision and also returns doubles. """ - try: - mask = np.array(field.mask, dtype=np.int32) + masked = np.ma.is_masked(field) + if masked: + mask = np.array(np.ma.getmaskarray(field), dtype=np.int32) field = np.ma.array(field, ndmin=1, dtype=np.double) - masked = True - except AttributeError: - mask = None - field = np.array(field, ndmin=1, dtype=np.double) - masked = False - - if direction == "x": - axis_to_swap = 0 - elif direction == "y": - axis_to_swap = 1 - elif direction == "z": - axis_to_swap = 2 else: + field = np.array(field, ndmin=1, dtype=np.double) + + axis_to_swap = AXIS_DIR.get(direction) + if axis_to_swap is None: raise ValueError("Unknown direction {0}".format(direction)) + # desired axis first, fill up field with empty dimensions up to a no. of 3 field = field.swapaxes(0, axis_to_swap) - if masked: - mask = mask.swapaxes(0, axis_to_swap) - - cython_estimator = _set_estimator(estimator) - - # fill up the field with empty dimensions up to a number of 3 + mask = mask.swapaxes(0, axis_to_swap) if masked else None for __ in range(3 - len(field.shape)): field = field[..., np.newaxis] if masked: for __ in range(3 - len(mask.shape)): mask = mask[..., np.newaxis] - if mask is None: - gamma = structured(field, cython_estimator) - else: - gamma = ma_structured(field, mask, cython_estimator) - return gamma + cython_estimator = _set_estimator(estimator) + + if masked: + return ma_structured(field, mask, cython_estimator) + return structured(field, cython_estimator) # for backward compatibility From 73801c94df1400337021d5a12ac52f6936dbca0d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 8 Nov 2020 14:09:29 +0100 Subject: [PATCH 18/44] vario: allow masked array for field in vario_estimate; add mask parameter to provide an external mask --- gstools/variogram/variogram.py | 42 ++++++++++++++++++++++++++-------- 1 file changed, 33 insertions(+), 9 deletions(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 9f935af6..5c39669b 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -59,6 +59,7 @@ def vario_estimate( angles_tol=np.pi / 8, bandwidth=None, no_data=np.nan, + mask=np.ma.nomask, mesh_type="unstructured", return_counts=False, ): @@ -98,10 +99,11 @@ def vario_estimate( Parameters ---------- pos : :class:`list` - the position tuple, containing main direction and transversal - directions + the position tuple, containing either the point coordinates (x, y, ...) + or the axes descriptions (for mesh_type='structured') field : :class:`numpy.ndarray` or :class:`list` of :class:`numpy.ndarray` The spatially distributed data. + Can also be of type :class:`numpy.ma.MaskedArray` to use masked values. You can pass a list of fields, that will be used simultaneously. This could be helpful, when there are multiple realizations at the same points, with the same statistical properties. @@ -153,12 +155,15 @@ def vario_estimate( Default: :any:`None` no_data : :class:`float`, optional Value to identify missing data in the given field. - Default: `np.nan` + Default: `numpy.nan` + mask : :class:`numpy.ndarray` of :class:`bool`, optional + Mask to deselect data in the given field. + Default: :any:`numpy.ma.nomask` mesh_type : :class:`str`, optional 'structured' / 'unstructured', indicates whether the pos tuple describes the axis or the point coordinates. Default: `'unstructured'` - return_counts: class:`bool`, optional + return_counts: :class:`bool`, optional if set to true, this function will also return the number of data points found at each lag distance as a third return value Default: False @@ -176,7 +181,13 @@ def vario_estimate( Internally uses double precision and also returns doubles. """ # allow multiple fields at same positions (ndmin=2: first axis -> field ID) - field = np.array(field, ndmin=2, dtype=np.double) + # need to convert to ma.array, since list of ma.array is not recognised + field = np.ma.array(field, ndmin=2, dtype=np.double) + masked = np.ma.is_masked(field) or np.any(mask) + if masked and np.all(mask): + raise ValueError("Given mask, masks all values.") + elif not masked: + field = field.filled() bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double) x, y, z, dim = pos2xyz(pos, calc_dim=True, dtype=np.double) bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 @@ -188,6 +199,21 @@ def vario_estimate( field = field.reshape((-1, len(x))) except ValueError: raise ValueError("'field' has wrong shape") + # prepare positions + pos = np.array(xyz2pos(x, y, z, dtype=np.double, max_dim=dim)) + # apply mask if wanted + if masked: + # if fields have different masks, take the minimal common mask + # given mask will be applied in addition + if np.size(mask) > 1: # not only np.ma.nomask + mask = np.logical_or( + mask.reshape(len(x)), np.all(field.mask, axis=0) + ) + else: + mask = np.all(field.mask, axis=0) + pos = pos[:, ~mask] + field.fill_value = np.nan # use no-data val. for remaining masked vals + field = field[:, ~mask].filled() # convert to ndarray # set no_data values if not np.isnan(no_data): field[np.isclose(field, float(no_data))] = np.nan @@ -214,8 +240,6 @@ def vario_estimate( # negative bandwidth to turn it off bandwidth = float(bandwidth) if bandwidth is not None else -1.0 angles_tol = float(angles_tol) - # prepare positions - pos = np.array(xyz2pos(x, y, z, dtype=np.double, max_dim=dim)) # prepare sampled variogram if sampling_size is not None and sampling_size < len(x): sampled_idx = np.random.RandomState(sampling_seed).choice( @@ -278,8 +302,8 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"): Parameters ---------- - field : :class:`numpy.ndarray` - the spatially distributed data + field : :class:`numpy.ndarray` or :class:`numpy.ma.MaskedArray` + the spatially distributed data (can be masked) direction : :class:`str` the axis over which the variogram will be estimated (x, y, z) estimator : :class:`str`, optional From 0abedeadc88283eae518cf1af3b3efae2e9ada84 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 8 Nov 2020 14:59:11 +0100 Subject: [PATCH 19/44] vario: add no_data option to vario_estimate_structured (solves #83) --- gstools/variogram/variogram.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 5c39669b..1f4a0a04 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -276,7 +276,9 @@ def vario_estimate( return bin_centres, estimates -def vario_estimate_structured(field, direction="x", estimator="matheron"): +def vario_estimate_structured( + field, direction="x", estimator="matheron", no_data=np.nan +): r"""Estimates the variogram on a regular grid. The indices of the given direction are used for the bins. @@ -314,6 +316,10 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"): Default: "matheron" + no_data : :class:`float`, optional + Value to identify missing data in the given field. + Default: `numpy.nan` + Returns ------- :class:`numpy.ndarray` @@ -327,12 +333,19 @@ def vario_estimate_structured(field, direction="x", estimator="matheron"): ----- Internally uses double precision and also returns doubles. """ - masked = np.ma.is_masked(field) + missing_mask = ( + np.isnan(field) if np.isnan(no_data) else np.isclose(field, no_data) + ) + missing = np.any(missing_mask) + masked = np.ma.is_masked(field) or missing if masked: - mask = np.array(np.ma.getmaskarray(field), dtype=np.int32) 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) else: field = np.array(field, ndmin=1, dtype=np.double) + missing_mask = None # free space axis_to_swap = AXIS_DIR.get(direction) if axis_to_swap is None: From 5b3dae2db7ba661086d92be84a7c689386eaf40a Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 8 Nov 2020 20:51:16 +0100 Subject: [PATCH 20/44] vario: allow arbitrary dim. in structured estimation --- gstools/variogram/estimator.pyx | 40 +++++++++++++++------------------ gstools/variogram/variogram.py | 19 +++++++--------- 2 files changed, 26 insertions(+), 33 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 34545664..ba086dfd 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -262,7 +262,7 @@ def unstructured( return np.asarray(variogram), np.asarray(counts) -def structured(const double[:,:,:] f, str estimator_type='m'): +def structured(const double[:,:] f, str estimator_type='m'): cdef _estimator_func estimator_func = choose_estimator_func(estimator_type) cdef _normalization_func normalization_func = ( choose_estimator_normalization(estimator_type) @@ -270,28 +270,26 @@ def structured(const double[:,:,:] f, str estimator_type='m'): cdef int i_max = f.shape[0] - 1 cdef int j_max = f.shape[1] - cdef int k_max = f.shape[2] - cdef int l_max = i_max + 1 + cdef int k_max = i_max + 1 - cdef vector[double] variogram = vector[double](l_max, 0.0) - cdef vector[long] counts = vector[long](l_max, 0) - cdef int i, j, k, l + cdef vector[double] variogram = vector[double](k_max, 0.0) + cdef vector[long] counts = vector[long](k_max, 0) + cdef int i, j, k with nogil, parallel(): for i in range(i_max): for j in range(j_max): - for k in range(k_max): - for l in prange(1, l_max-i): - counts[l] += 1 - variogram[l] += estimator_func(f[i,j,k] - f[i+l,j,k]) + for k in prange(1, k_max-i): + counts[k] += 1 + variogram[k] += estimator_func(f[i,j] - f[i+k,j]) normalization_func(variogram, counts) return np.asarray(variogram) def ma_structured( - const double[:,:,:] f, - const bint[:,:,:] mask, + const double[:,:] f, + const bint[:,:] mask, str estimator_type='m' ): cdef _estimator_func estimator_func = choose_estimator_func(estimator_type) @@ -301,21 +299,19 @@ def ma_structured( cdef int i_max = f.shape[0] - 1 cdef int j_max = f.shape[1] - cdef int k_max = f.shape[2] - cdef int l_max = i_max + 1 + cdef int k_max = i_max + 1 - cdef vector[double] variogram = vector[double](l_max, 0.0) - cdef vector[long] counts = vector[long](l_max, 0) - cdef int i, j, k, l + cdef vector[double] variogram = vector[double](k_max, 0.0) + cdef vector[long] counts = vector[long](k_max, 0) + cdef int i, j, k with nogil, parallel(): for i in range(i_max): for j in range(j_max): - for k in range(k_max): - for l in prange(1, l_max-i): - if not mask[i,j,k] and not mask[i+l,j,k]: - counts[l] += 1 - variogram[l] += estimator_func(f[i,j,k] - f[i+l,j,k]) + for k in prange(1, k_max-i): + if not mask[i,j] and not mask[i+k,j]: + counts[k] += 1 + variogram[k] += estimator_func(f[i,j] - f[i+k,j]) normalization_func(variogram, counts) return np.asarray(variogram) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 1f4a0a04..74231fc8 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -31,6 +31,7 @@ ] +AXIS = ["x", "y", "z"] AXIS_DIR = {"x": 0, "y": 1, "z": 2} @@ -306,8 +307,9 @@ def vario_estimate_structured( ---------- field : :class:`numpy.ndarray` or :class:`numpy.ma.MaskedArray` the spatially distributed data (can be masked) - direction : :class:`str` + direction : :class:`str` or :class:`int` the axis over which the variogram will be estimated (x, y, z) + or (0, 1, 2, ...) estimator : :class:`str`, optional the estimator function, possible choices: @@ -347,18 +349,13 @@ def vario_estimate_structured( field = np.array(field, ndmin=1, dtype=np.double) missing_mask = None # free space - axis_to_swap = AXIS_DIR.get(direction) - if axis_to_swap is None: - raise ValueError("Unknown direction {0}".format(direction)) - - # desired axis first, fill up field with empty dimensions up to a no. of 3 + axis_to_swap = AXIS_DIR[direction] if direction in AXIS else int(direction) + # desired axis first, convert to 2D array afterwards field = field.swapaxes(0, axis_to_swap) - mask = mask.swapaxes(0, axis_to_swap) if masked else None - for __ in range(3 - len(field.shape)): - field = field[..., np.newaxis] + field = field.reshape((field.shape[0], -1)) if masked: - for __ in range(3 - len(mask.shape)): - mask = mask[..., np.newaxis] + mask = mask.swapaxes(0, axis_to_swap) + mask = mask.reshape((mask.shape[0], -1)) cython_estimator = _set_estimator(estimator) From 834356621fec38c7bbc06cb92a2d569b556c1a93 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 8 Nov 2020 23:30:51 +0100 Subject: [PATCH 21/44] vario: rename vario_estimate_structured to vario_estimate_axis --- .../03_variogram/01_variogram_estimation.py | 4 +- gstools/__init__.py | 5 +- gstools/variogram/__init__.py | 5 +- gstools/variogram/variogram.py | 11 +-- tests/test_variogram_structured.py | 82 +++++++++---------- 5 files changed, 55 insertions(+), 52 deletions(-) diff --git a/examples/03_variogram/01_variogram_estimation.py b/examples/03_variogram/01_variogram_estimation.py index f0574adf..93817260 100644 --- a/examples/03_variogram/01_variogram_estimation.py +++ b/examples/03_variogram/01_variogram_estimation.py @@ -209,8 +209,8 @@ def generate_transmissivity(): # With this much smaller data set, we can immediately estimate the variogram in # the x- and y-axis -gamma_x = gs.vario_estimate_structured(herten_trans_skip, direction="x") -gamma_y = gs.vario_estimate_structured(herten_trans_skip, direction="y") +gamma_x = gs.vario_estimate_axis(herten_trans_skip, direction="x") +gamma_y = gs.vario_estimate_axis(herten_trans_skip, direction="y") ############################################################################### # With these two estimated variograms, we can start fitting :any:`Exponential` diff --git a/gstools/__init__.py b/gstools/__init__.py index ca0b6950..d6eb19b4 100644 --- a/gstools/__init__.py +++ b/gstools/__init__.py @@ -96,8 +96,7 @@ .. autosummary:: vario_estimate - vario_estimate_structured - vario_estimate_unstructured + vario_estimate_axis """ from gstools import field, variogram, random, covmodel, tools, krige, transform @@ -112,6 +111,7 @@ ) from gstools.variogram import ( vario_estimate, + vario_estimate_axis, vario_estimate_structured, vario_estimate_unstructured, ) @@ -158,6 +158,7 @@ __all__ += [ "vario_estimate", + "vario_estimate_axis", "vario_estimate_structured", "vario_estimate_unstructured", ] diff --git a/gstools/variogram/__init__.py b/gstools/variogram/__init__.py index 7f496713..8fb0444d 100644 --- a/gstools/variogram/__init__.py +++ b/gstools/variogram/__init__.py @@ -9,20 +9,21 @@ .. autosummary:: vario_estimate - vario_estimate_unstructured - vario_estimate_structured + vario_estimate_axis ---- """ from gstools.variogram.variogram import ( vario_estimate, + vario_estimate_axis, vario_estimate_structured, vario_estimate_unstructured, ) __all__ = [ "vario_estimate", + "vario_estimate_axis", "vario_estimate_unstructured", "vario_estimate_structured", ] diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 74231fc8..5c64ce6d 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -8,8 +8,7 @@ .. autosummary:: vario_estimate - vario_estimate_unstructured - vario_estimate_structured + vario_estimate_axis """ # pylint: disable=C0103 @@ -26,6 +25,7 @@ __all__ = [ "vario_estimate", + "vario_estimate_axis", "vario_estimate_unstructured", "vario_estimate_structured", ] @@ -187,7 +187,7 @@ def vario_estimate( masked = np.ma.is_masked(field) or np.any(mask) if masked and np.all(mask): raise ValueError("Given mask, masks all values.") - elif not masked: + if not masked: field = field.filled() bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double) x, y, z, dim = pos2xyz(pos, calc_dim=True, dtype=np.double) @@ -277,10 +277,10 @@ def vario_estimate( return bin_centres, estimates -def vario_estimate_structured( +def vario_estimate_axis( field, direction="x", estimator="matheron", no_data=np.nan ): - r"""Estimates the variogram on a regular grid. + r"""Estimates the variogram along array axis. The indices of the given direction are used for the bins. The algorithm calculates following equation: @@ -366,3 +366,4 @@ def vario_estimate_structured( # for backward compatibility vario_estimate_unstructured = vario_estimate +vario_estimate_structured = vario_estimate_axis diff --git a/tests/test_variogram_structured.py b/tests/test_variogram_structured.py index e7ff0672..ca1f94f2 100644 --- a/tests/test_variogram_structured.py +++ b/tests/test_variogram_structured.py @@ -17,17 +17,17 @@ def test_doubles(self): (41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3), dtype=np.double, ) - gamma = variogram.vario_estimate_structured(z) + gamma = variogram.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 0.4917, places=4) def test_ints(self): z = np.array((10, 20, 30, 40), dtype=int) - gamma = variogram.vario_estimate_structured(z) + gamma = variogram.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 50.0, places=4) def test_np_int(self): z = np.array((10, 20, 30, 40), dtype=np.int) - gamma = variogram.vario_estimate_structured(z) + gamma = variogram.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 50.0, places=4) def test_mixed(self): @@ -35,22 +35,22 @@ def test_mixed(self): (41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3), dtype=np.double, ) - gamma = variogram.vario_estimate_structured(z) + gamma = variogram.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 0.4917, places=4) z = np.array((10, 20, 30, 40), dtype=int) - gamma = variogram.vario_estimate_structured(z) + gamma = variogram.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 50.0, places=4) def test_list(self): z = [41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3] - gamma = variogram.vario_estimate_structured(z) + gamma = variogram.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 0.4917, places=4) def test_cressie_1d(self): z = [41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3] - gamma = variogram.vario_estimate_structured(z, estimator="cressie") + gamma = variogram.vario_estimate_axis(z, estimator="cressie") self.assertAlmostEqual(gamma[1], 1.546 / 2.0, places=3) def test_1d(self): @@ -59,7 +59,7 @@ def test_1d(self): (41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3), dtype=np.double, ) - gamma = variogram.vario_estimate_structured(z) + gamma = variogram.vario_estimate_axis(z) self.assertAlmostEqual(gamma[0], 0.0000, places=4) self.assertAlmostEqual(gamma[1], 0.4917, places=4) self.assertAlmostEqual(gamma[2], 0.7625, places=4) @@ -71,12 +71,12 @@ def test_masked_1d(self): dtype=np.double, ) z_ma = np.ma.masked_array(z, mask=[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]) - gamma = variogram.vario_estimate_structured(z_ma) + gamma = variogram.vario_estimate_axis(z_ma) self.assertAlmostEqual(gamma[0], 0.0000, places=4) self.assertAlmostEqual(gamma[1], 0.4917, places=4) self.assertAlmostEqual(gamma[2], 0.7625, places=4) z_ma = np.ma.masked_array(z, mask=[1, 0, 0, 0, 0, 0, 0, 0, 0, 0]) - gamma = variogram.vario_estimate_structured(z_ma) + gamma = variogram.vario_estimate_axis(z_ma) self.assertAlmostEqual(gamma[0], 0.0000, places=4) self.assertAlmostEqual(gamma[1], 0.4906, places=4) self.assertAlmostEqual(gamma[2], 0.7107, places=4) @@ -87,8 +87,8 @@ def test_masked_2d(self): mask = np.zeros_like(field) field_ma = np.ma.masked_array(field, mask=mask) - gamma_x = variogram.vario_estimate_structured(field_ma, direction="x") - gamma_y = variogram.vario_estimate_structured(field_ma, direction="y") + gamma_x = variogram.vario_estimate_axis(field_ma, direction="x") + gamma_y = variogram.vario_estimate_axis(field_ma, direction="y") var = 1.0 / 12.0 self.assertAlmostEqual(gamma_x[0], 0.0, places=2) @@ -101,8 +101,8 @@ def test_masked_2d(self): mask = np.zeros_like(field) mask[0, 0] = 1 field = np.ma.masked_array(field, mask=mask) - gamma_x = variogram.vario_estimate_structured(field_ma, direction="x") - gamma_y = variogram.vario_estimate_structured(field_ma, direction="y") + gamma_x = variogram.vario_estimate_axis(field_ma, direction="x") + gamma_y = variogram.vario_estimate_axis(field_ma, direction="y") self.assertAlmostEqual(gamma_x[0], 0.0, places=2) self.assertAlmostEqual(gamma_y[0], 0.0, places=2) @@ -112,9 +112,9 @@ def test_masked_3d(self): mask = np.zeros_like(field) field_ma = np.ma.masked_array(field, mask=mask) - gamma_x = variogram.vario_estimate_structured(field_ma, direction="x") - gamma_y = variogram.vario_estimate_structured(field_ma, direction="y") - gamma_z = variogram.vario_estimate_structured(field_ma, direction="z") + gamma_x = variogram.vario_estimate_axis(field_ma, direction="x") + gamma_y = variogram.vario_estimate_axis(field_ma, direction="y") + gamma_z = variogram.vario_estimate_axis(field_ma, direction="z") var = 1.0 / 12.0 self.assertAlmostEqual(gamma_x[0], 0.0, places=2) @@ -130,9 +130,9 @@ def test_masked_3d(self): mask = np.zeros_like(field) mask[0, 0, 0] = 1 field = np.ma.masked_array(field, mask=mask) - gamma_x = variogram.vario_estimate_structured(field_ma, direction="x") - gamma_y = variogram.vario_estimate_structured(field_ma, direction="y") - gamma_z = variogram.vario_estimate_structured(field_ma, direction="z") + gamma_x = variogram.vario_estimate_axis(field_ma, direction="x") + gamma_y = variogram.vario_estimate_axis(field_ma, direction="y") + gamma_z = variogram.vario_estimate_axis(field_ma, direction="z") self.assertAlmostEqual(gamma_x[0], 0.0, places=2) self.assertAlmostEqual(gamma_y[0], 0.0, places=2) self.assertAlmostEqual(gamma_z[0], 0.0, places=2) @@ -144,8 +144,8 @@ def test_uncorrelated_2d(self): rng = np.random.RandomState(1479373475) field = rng.rand(len(x), len(y)) - gamma_x = variogram.vario_estimate_structured(field, direction="x") - gamma_y = variogram.vario_estimate_structured(field, direction="y") + gamma_x = variogram.vario_estimate_axis(field, direction="x") + gamma_y = variogram.vario_estimate_axis(field, direction="y") var = 1.0 / 12.0 self.assertAlmostEqual(gamma_x[0], 0.0, places=2) @@ -162,10 +162,10 @@ def test_uncorrelated_cressie_2d(self): rng = np.random.RandomState(1479373475) field = rng.rand(len(x), len(y)) - gamma_x = variogram.vario_estimate_structured( + gamma_x = variogram.vario_estimate_axis( field, direction="x", estimator="cressie" ) - gamma_y = variogram.vario_estimate_structured( + gamma_y = variogram.vario_estimate_axis( field, direction="y", estimator="cressie" ) @@ -183,9 +183,9 @@ def test_uncorrelated_3d(self): rng = np.random.RandomState(1479373475) field = rng.rand(len(x), len(y), len(z)) - gamma = variogram.vario_estimate_structured(field, "x") - gamma = variogram.vario_estimate_structured(field, "y") - gamma = variogram.vario_estimate_structured(field, "z") + gamma = variogram.vario_estimate_axis(field, "x") + gamma = variogram.vario_estimate_axis(field, "y") + gamma = variogram.vario_estimate_axis(field, "z") var = 1.0 / 12.0 self.assertAlmostEqual(gamma[0], 0.0, places=2) @@ -203,11 +203,11 @@ def test_directions_2d(self): # random values repeated along x-axis field_y = np.tile(y_rand, (len(x), 1)) - gamma_x_x = variogram.vario_estimate_structured(field_x, direction="x") - gamma_x_y = variogram.vario_estimate_structured(field_x, direction="y") + gamma_x_x = variogram.vario_estimate_axis(field_x, direction="x") + gamma_x_y = variogram.vario_estimate_axis(field_x, direction="y") - gamma_y_x = variogram.vario_estimate_structured(field_y, direction="x") - gamma_y_y = variogram.vario_estimate_structured(field_y, direction="y") + gamma_y_x = variogram.vario_estimate_axis(field_y, direction="x") + gamma_y_y = variogram.vario_estimate_axis(field_y, direction="y") self.assertAlmostEqual(gamma_x_y[1], 0.0) self.assertAlmostEqual(gamma_x_y[len(gamma_x_y) // 2], 0.0) @@ -229,17 +229,17 @@ def test_directions_3d(self): field_y = np.tile(y_rand.reshape((1, len(y), 1)), (len(x), 1, len(z))) field_z = np.tile(z_rand.reshape((1, 1, len(z))), (len(x), len(y), 1)) - gamma_x_x = variogram.vario_estimate_structured(field_x, direction="x") - gamma_x_y = variogram.vario_estimate_structured(field_x, direction="y") - gamma_x_z = variogram.vario_estimate_structured(field_x, direction="z") + gamma_x_x = variogram.vario_estimate_axis(field_x, direction="x") + gamma_x_y = variogram.vario_estimate_axis(field_x, direction="y") + gamma_x_z = variogram.vario_estimate_axis(field_x, direction="z") - gamma_y_x = variogram.vario_estimate_structured(field_y, direction="x") - gamma_y_y = variogram.vario_estimate_structured(field_y, direction="y") - gamma_y_z = variogram.vario_estimate_structured(field_y, direction="z") + gamma_y_x = variogram.vario_estimate_axis(field_y, direction="x") + gamma_y_y = variogram.vario_estimate_axis(field_y, direction="y") + gamma_y_z = variogram.vario_estimate_axis(field_y, direction="z") - gamma_z_x = variogram.vario_estimate_structured(field_z, direction="x") - gamma_z_y = variogram.vario_estimate_structured(field_z, direction="y") - gamma_z_z = variogram.vario_estimate_structured(field_z, direction="z") + gamma_z_x = variogram.vario_estimate_axis(field_z, direction="x") + gamma_z_y = variogram.vario_estimate_axis(field_z, direction="y") + gamma_z_z = variogram.vario_estimate_axis(field_z, direction="z") self.assertAlmostEqual(gamma_x_y[1], 0.0) self.assertAlmostEqual(gamma_x_y[len(gamma_x_y) // 2], 0.0) @@ -265,7 +265,7 @@ def test_exceptions(self): rng = np.random.RandomState(1479373475) x_rand = rng.rand(len(x)) self.assertRaises( - ValueError, variogram.vario_estimate_structured, x, "a" + ValueError, variogram.vario_estimate_axis, x, "a" ) From 036667d4878a40d2043b96372b17c18e7cbc43ec Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sun, 8 Nov 2020 23:46:41 +0100 Subject: [PATCH 22/44] vario: don't alter count array, if count==0 --- gstools/variogram/estimator.pyx | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index ba086dfd..3028a3b1 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -82,22 +82,20 @@ cdef inline void normalization_matheron( cdef int i for i in range(variogram.size()): # avoid division by zero - if counts[i] == 0: - counts[i] = 1 - variogram[i] /= (2. * counts[i]) + variogram[i] /= (2. * max(counts[i], 1)) cdef inline void normalization_cressie( vector[double]& variogram, vector[long]& counts ): cdef int i + cdef long cnt for i in range(variogram.size()): # avoid division by zero - if counts[i] == 0: - counts[i] = 1 + cnt = max(counts[i], 1) variogram[i] = ( - 0.5 * (1./counts[i] * variogram[i])**4 / - (0.457 + 0.494 / counts[i] + 0.045 / counts[i]**2) + 0.5 * (1./cnt * variogram[i])**4 / + (0.457 + 0.494 / cnt + 0.045 / cnt**2) ) ctypedef void (*_normalization_func)( @@ -113,23 +111,21 @@ cdef inline void normalization_matheron_vec( for d in range(variogram.shape[0]): for i in range(variogram.shape[1]): # avoid division by zero - if counts[d, i] == 0: - counts[d, i] = 1 - variogram[d, i] /= (2. * counts[d, i]) + variogram[d, i] /= (2. * max(counts[d, i], 1)) cdef inline void normalization_cressie_vec( double[:,:]& variogram, long[:,:]& counts ): cdef int d, i + cdef long cnt for d in range(variogram.shape[0]): for i in range(variogram.shape[1]): # avoid division by zero - if counts[d, i] == 0: - counts[d, i] = 1 + cnt = max(counts[d, i], 1) variogram[d, i] = ( - 0.5 * (1./counts[d, i] * variogram[d, i])**4 / - (0.457 + 0.494 / counts[d, i] + 0.045 / counts[d, i]**2) + 0.5 * (1./cnt * variogram[d, i])**4 / + (0.457 + 0.494 / cnt + 0.045 / cnt**2) ) ctypedef void (*_normalization_func_vec)( From 9050cab1c14740917e7f8a40b7948700e39d9ec3 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 9 Nov 2020 11:52:03 +0100 Subject: [PATCH 23/44] vario: use np.invert to create selection from mask; better Exception forwarding; doc update --- gstools/variogram/variogram.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 5c64ce6d..d27ed661 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -198,23 +198,25 @@ def vario_estimate( if len(field.shape) > 2 or field.shape[1] != len(x): try: field = field.reshape((-1, len(x))) - except ValueError: - raise ValueError("'field' has wrong shape") + except ValueError as exc: + raise ValueError("'field' has wrong shape") from exc # prepare positions pos = np.array(xyz2pos(x, y, z, dtype=np.double, max_dim=dim)) # apply mask if wanted if masked: # if fields have different masks, take the minimal common mask # given mask will be applied in addition + # selected region is the inverted masked (unmasked values) if np.size(mask) > 1: # not only np.ma.nomask - mask = np.logical_or( - mask.reshape(len(x)), np.all(field.mask, axis=0) + select = np.invert( + np.logical_or(mask.reshape(len(x)), np.all(field.mask, axis=0)) ) else: - mask = np.all(field.mask, axis=0) - pos = pos[:, ~mask] + select = np.invert(np.all(field.mask, axis=0)) + pos = pos[:, select] field.fill_value = np.nan # use no-data val. for remaining masked vals - field = field[:, ~mask].filled() # convert to ndarray + field = field[:, select].filled() # convert to ndarray + select = mask = None # free space # set no_data values if not np.isnan(no_data): field[np.isclose(field, float(no_data))] = np.nan @@ -283,6 +285,8 @@ def vario_estimate_axis( r"""Estimates the variogram along array axis. The indices of the given direction are used for the bins. + Uniform spacings along the given axis are assumed. + The algorithm calculates following equation: .. math:: From 47fe4382a81145f6e8e63efd13fa811522127a20 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 9 Nov 2020 12:14:44 +0100 Subject: [PATCH 24/44] tests: add test comparing struc/unstruc vario estimate along axis with missing data --- tests/test_variogram_structured.py | 86 +++++++++++---------- tests/test_variogram_unstructured.py | 109 +++++++++++---------------- 2 files changed, 86 insertions(+), 109 deletions(-) diff --git a/tests/test_variogram_structured.py b/tests/test_variogram_structured.py index ca1f94f2..73f566c4 100644 --- a/tests/test_variogram_structured.py +++ b/tests/test_variogram_structured.py @@ -5,7 +5,7 @@ import unittest import numpy as np -from gstools import variogram +import gstools as gs class TestVariogramstructured(unittest.TestCase): @@ -17,17 +17,17 @@ def test_doubles(self): (41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3), dtype=np.double, ) - gamma = variogram.vario_estimate_axis(z) + gamma = gs.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 0.4917, places=4) def test_ints(self): z = np.array((10, 20, 30, 40), dtype=int) - gamma = variogram.vario_estimate_axis(z) + gamma = gs.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 50.0, places=4) def test_np_int(self): z = np.array((10, 20, 30, 40), dtype=np.int) - gamma = variogram.vario_estimate_axis(z) + gamma = gs.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 50.0, places=4) def test_mixed(self): @@ -35,22 +35,22 @@ def test_mixed(self): (41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3), dtype=np.double, ) - gamma = variogram.vario_estimate_axis(z) + gamma = gs.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 0.4917, places=4) z = np.array((10, 20, 30, 40), dtype=int) - gamma = variogram.vario_estimate_axis(z) + gamma = gs.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 50.0, places=4) def test_list(self): z = [41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3] - gamma = variogram.vario_estimate_axis(z) + gamma = gs.vario_estimate_axis(z) self.assertAlmostEqual(gamma[1], 0.4917, places=4) def test_cressie_1d(self): z = [41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3] - gamma = variogram.vario_estimate_axis(z, estimator="cressie") + gamma = gs.vario_estimate_axis(z, estimator="cressie") self.assertAlmostEqual(gamma[1], 1.546 / 2.0, places=3) def test_1d(self): @@ -59,7 +59,7 @@ def test_1d(self): (41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3), dtype=np.double, ) - gamma = variogram.vario_estimate_axis(z) + gamma = gs.vario_estimate_axis(z) self.assertAlmostEqual(gamma[0], 0.0000, places=4) self.assertAlmostEqual(gamma[1], 0.4917, places=4) self.assertAlmostEqual(gamma[2], 0.7625, places=4) @@ -71,12 +71,12 @@ def test_masked_1d(self): dtype=np.double, ) z_ma = np.ma.masked_array(z, mask=[0, 0, 0, 0, 0, 0, 0, 0, 0, 0]) - gamma = variogram.vario_estimate_axis(z_ma) + gamma = gs.vario_estimate_axis(z_ma) self.assertAlmostEqual(gamma[0], 0.0000, places=4) self.assertAlmostEqual(gamma[1], 0.4917, places=4) self.assertAlmostEqual(gamma[2], 0.7625, places=4) z_ma = np.ma.masked_array(z, mask=[1, 0, 0, 0, 0, 0, 0, 0, 0, 0]) - gamma = variogram.vario_estimate_axis(z_ma) + gamma = gs.vario_estimate_axis(z_ma) self.assertAlmostEqual(gamma[0], 0.0000, places=4) self.assertAlmostEqual(gamma[1], 0.4906, places=4) self.assertAlmostEqual(gamma[2], 0.7107, places=4) @@ -87,8 +87,8 @@ def test_masked_2d(self): mask = np.zeros_like(field) field_ma = np.ma.masked_array(field, mask=mask) - gamma_x = variogram.vario_estimate_axis(field_ma, direction="x") - gamma_y = variogram.vario_estimate_axis(field_ma, direction="y") + gamma_x = gs.vario_estimate_axis(field_ma, direction="x") + gamma_y = gs.vario_estimate_axis(field_ma, direction="y") var = 1.0 / 12.0 self.assertAlmostEqual(gamma_x[0], 0.0, places=2) @@ -101,8 +101,8 @@ def test_masked_2d(self): mask = np.zeros_like(field) mask[0, 0] = 1 field = np.ma.masked_array(field, mask=mask) - gamma_x = variogram.vario_estimate_axis(field_ma, direction="x") - gamma_y = variogram.vario_estimate_axis(field_ma, direction="y") + gamma_x = gs.vario_estimate_axis(field_ma, direction="x") + gamma_y = gs.vario_estimate_axis(field_ma, direction="y") self.assertAlmostEqual(gamma_x[0], 0.0, places=2) self.assertAlmostEqual(gamma_y[0], 0.0, places=2) @@ -112,9 +112,9 @@ def test_masked_3d(self): mask = np.zeros_like(field) field_ma = np.ma.masked_array(field, mask=mask) - gamma_x = variogram.vario_estimate_axis(field_ma, direction="x") - gamma_y = variogram.vario_estimate_axis(field_ma, direction="y") - gamma_z = variogram.vario_estimate_axis(field_ma, direction="z") + gamma_x = gs.vario_estimate_axis(field_ma, direction="x") + gamma_y = gs.vario_estimate_axis(field_ma, direction="y") + gamma_z = gs.vario_estimate_axis(field_ma, direction="z") var = 1.0 / 12.0 self.assertAlmostEqual(gamma_x[0], 0.0, places=2) @@ -130,9 +130,9 @@ def test_masked_3d(self): mask = np.zeros_like(field) mask[0, 0, 0] = 1 field = np.ma.masked_array(field, mask=mask) - gamma_x = variogram.vario_estimate_axis(field_ma, direction="x") - gamma_y = variogram.vario_estimate_axis(field_ma, direction="y") - gamma_z = variogram.vario_estimate_axis(field_ma, direction="z") + gamma_x = gs.vario_estimate_axis(field_ma, direction="x") + gamma_y = gs.vario_estimate_axis(field_ma, direction="y") + gamma_z = gs.vario_estimate_axis(field_ma, direction="z") self.assertAlmostEqual(gamma_x[0], 0.0, places=2) self.assertAlmostEqual(gamma_y[0], 0.0, places=2) self.assertAlmostEqual(gamma_z[0], 0.0, places=2) @@ -144,8 +144,8 @@ def test_uncorrelated_2d(self): rng = np.random.RandomState(1479373475) field = rng.rand(len(x), len(y)) - gamma_x = variogram.vario_estimate_axis(field, direction="x") - gamma_y = variogram.vario_estimate_axis(field, direction="y") + gamma_x = gs.vario_estimate_axis(field, direction="x") + gamma_y = gs.vario_estimate_axis(field, direction="y") var = 1.0 / 12.0 self.assertAlmostEqual(gamma_x[0], 0.0, places=2) @@ -162,10 +162,10 @@ def test_uncorrelated_cressie_2d(self): rng = np.random.RandomState(1479373475) field = rng.rand(len(x), len(y)) - gamma_x = variogram.vario_estimate_axis( + gamma_x = gs.vario_estimate_axis( field, direction="x", estimator="cressie" ) - gamma_y = variogram.vario_estimate_axis( + gamma_y = gs.vario_estimate_axis( field, direction="y", estimator="cressie" ) @@ -183,9 +183,9 @@ def test_uncorrelated_3d(self): rng = np.random.RandomState(1479373475) field = rng.rand(len(x), len(y), len(z)) - gamma = variogram.vario_estimate_axis(field, "x") - gamma = variogram.vario_estimate_axis(field, "y") - gamma = variogram.vario_estimate_axis(field, "z") + gamma = gs.vario_estimate_axis(field, "x") + gamma = gs.vario_estimate_axis(field, "y") + gamma = gs.vario_estimate_axis(field, "z") var = 1.0 / 12.0 self.assertAlmostEqual(gamma[0], 0.0, places=2) @@ -203,11 +203,11 @@ def test_directions_2d(self): # random values repeated along x-axis field_y = np.tile(y_rand, (len(x), 1)) - gamma_x_x = variogram.vario_estimate_axis(field_x, direction="x") - gamma_x_y = variogram.vario_estimate_axis(field_x, direction="y") + gamma_x_x = gs.vario_estimate_axis(field_x, direction="x") + gamma_x_y = gs.vario_estimate_axis(field_x, direction="y") - gamma_y_x = variogram.vario_estimate_axis(field_y, direction="x") - gamma_y_y = variogram.vario_estimate_axis(field_y, direction="y") + gamma_y_x = gs.vario_estimate_axis(field_y, direction="x") + gamma_y_y = gs.vario_estimate_axis(field_y, direction="y") self.assertAlmostEqual(gamma_x_y[1], 0.0) self.assertAlmostEqual(gamma_x_y[len(gamma_x_y) // 2], 0.0) @@ -229,17 +229,17 @@ def test_directions_3d(self): field_y = np.tile(y_rand.reshape((1, len(y), 1)), (len(x), 1, len(z))) field_z = np.tile(z_rand.reshape((1, 1, len(z))), (len(x), len(y), 1)) - gamma_x_x = variogram.vario_estimate_axis(field_x, direction="x") - gamma_x_y = variogram.vario_estimate_axis(field_x, direction="y") - gamma_x_z = variogram.vario_estimate_axis(field_x, direction="z") + gamma_x_x = gs.vario_estimate_axis(field_x, direction="x") + gamma_x_y = gs.vario_estimate_axis(field_x, direction="y") + gamma_x_z = gs.vario_estimate_axis(field_x, direction="z") - gamma_y_x = variogram.vario_estimate_axis(field_y, direction="x") - gamma_y_y = variogram.vario_estimate_axis(field_y, direction="y") - gamma_y_z = variogram.vario_estimate_axis(field_y, direction="z") + gamma_y_x = gs.vario_estimate_axis(field_y, direction="x") + gamma_y_y = gs.vario_estimate_axis(field_y, direction="y") + gamma_y_z = gs.vario_estimate_axis(field_y, direction="z") - gamma_z_x = variogram.vario_estimate_axis(field_z, direction="x") - gamma_z_y = variogram.vario_estimate_axis(field_z, direction="y") - gamma_z_z = variogram.vario_estimate_axis(field_z, direction="z") + gamma_z_x = gs.vario_estimate_axis(field_z, direction="x") + gamma_z_y = gs.vario_estimate_axis(field_z, direction="y") + gamma_z_z = gs.vario_estimate_axis(field_z, direction="z") self.assertAlmostEqual(gamma_x_y[1], 0.0) self.assertAlmostEqual(gamma_x_y[len(gamma_x_y) // 2], 0.0) @@ -264,9 +264,7 @@ def test_exceptions(self): x = np.linspace(0.0, 10.0, 20) rng = np.random.RandomState(1479373475) x_rand = rng.rand(len(x)) - self.assertRaises( - ValueError, variogram.vario_estimate_axis, x, "a" - ) + self.assertRaises(ValueError, gs.vario_estimate_axis, x, "a") if __name__ == "__main__": diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 13c7cff3..60e037a1 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -5,45 +5,12 @@ import unittest import numpy as np -from gstools import vario_estimate, Exponential, SRF import gstools as gs class TestVariogramUnstructured(unittest.TestCase): def setUp(self): - - # this code generates the testdata for the 2D rotation test cases - x = np.random.RandomState(19970221).rand(20) * 10.0 - 5 - y = np.zeros_like(x) - model = gs.Exponential(dim=2, var=2, len_scale=8) - srf = gs.SRF(model, mean=0, seed=19970221) - field = srf((x, y)) - - bins = np.arange(10) - bin_center, gamma = gs.vario_estimate((x,), field, bins) - idx = np.argsort(x) - self.test_data_rotation_1 = { - "gamma": gamma, - "x": x[idx], - "field": field[idx], - "bins": bins, - "bin_center": bin_center, - } - - # CODE ABOVE SHOULD GENERATE THIS DATA - # x = np.array([ - # -4.86210059, -4.1984934 , -3.9246953 , -3.28490663, -2.16332379, - # -1.87553275, -1.74125124, -1.27224687, -1.20931578, -0.2413368 , - # 0.03200921, 1.17099773, 1.53863105, 1.64478688, 2.75252136, - # 3.3556915 , 3.89828775, 4.21485964, 4.5364357 , 4.79236969]), - # field = np.array([ - # -1.10318365, -0.53566629, -0.41789049, -1.06167529, 0.38449961, - # -0.36550477, -0.98905552, -0.19352766, 0.16264266, 0.26920833, - # 0.05379665, 0.71275006, 0.36651935, 0.17366865, 1.20022343, - # 0.79385446, 0.69456069, 1.0733393 , 0.71191592, 0.71969766]) - # gamma_exp = np.array([ - # 0.14260989, 0.18301197, 0.25855841, 0.29990083, 0.67914526, - # 0.60136535, 0.92875492, 1.46910435, 1.10165104]) + pass def test_doubles(self): x = np.arange(1, 11, 1, dtype=np.double) @@ -52,21 +19,21 @@ def test_doubles(self): dtype=np.double, ) bins = np.arange(1, 11, 1, dtype=np.double) - bin_centres, gamma = vario_estimate([x], z, bins) + bin_centres, gamma = gs.vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 0.4917, places=4) def test_ints(self): x = np.arange(1, 5, 1, dtype=int) z = np.array((10, 20, 30, 40), dtype=int) bins = np.arange(1, 11, 1, dtype=int) - bin_centres, gamma = vario_estimate([x], z, bins) + bin_centres, gamma = gs.vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 50.0, places=4) def test_np_int(self): x = np.arange(1, 5, 1, dtype=np.int) z = np.array((10, 20, 30, 40), dtype=np.int) bins = np.arange(1, 11, 1, dtype=np.int) - bin_centres, gamma = vario_estimate([x], z, bins) + bin_centres, gamma = gs.vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 50.0, places=4) def test_mixed(self): @@ -76,26 +43,26 @@ def test_mixed(self): dtype=np.double, ) bins = np.arange(1, 11, 1, dtype=int) - bin_centres, gamma = vario_estimate([x], z, bins) + bin_centres, gamma = gs.vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 0.4917, places=4) x = np.arange(1, 5, 1, dtype=np.double) z = np.array((10, 20, 30, 40), dtype=int) bins = np.arange(1, 11, 1, dtype=int) - bin_centres, gamma = vario_estimate([x], z, bins) + bin_centres, gamma = gs.vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 50.0, places=4) x = np.arange(1, 5, 1, dtype=np.double) z = np.array((10, 20, 30, 40), dtype=int) bins = np.arange(1, 11, 1, dtype=np.double) - bin_centres, gamma = vario_estimate([x], z, bins) + bin_centres, gamma = gs.vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 50.0, places=4) def test_list(self): x = np.arange(1, 11, 1, dtype=np.double) z = [41.2, 40.2, 39.7, 39.2, 40.1, 38.3, 39.1, 40.0, 41.1, 40.3] bins = np.arange(1, 11, 1, dtype=np.double) - bin_centres, gamma = vario_estimate([x], z, bins) + bin_centres, gamma = gs.vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[1], 0.7625, places=4) def test_1d(self): @@ -106,7 +73,7 @@ def test_1d(self): dtype=np.double, ) bins = np.arange(1, 11, 1, dtype=np.double) - bin_centres, gamma = vario_estimate([x], z, bins) + bin_centres, gamma = gs.vario_estimate([x], z, bins) self.assertAlmostEqual(gamma[0], 0.4917, places=4) self.assertAlmostEqual(gamma[1], 0.7625, places=4) @@ -122,7 +89,7 @@ def test_uncorrelated_2d(self): bins = np.arange(0, 100, 10) - bin_centres, gamma = vario_estimate((x, y), field, bins) + bin_centres, gamma = gs.vario_estimate((x, y), field, bins) var = 1.0 / 12.0 self.assertAlmostEqual(gamma[0], var, places=2) @@ -143,9 +110,7 @@ def test_uncorrelated_3d(self): bins = np.arange(0, 100, 10) - bin_centres, gamma = vario_estimate( - (x, y, z), field, bins - ) + bin_centres, gamma = gs.vario_estimate((x, y, z), field, bins) var = 1.0 / 12.0 self.assertAlmostEqual(gamma[0], var, places=2) @@ -160,7 +125,7 @@ def test_sampling_1d(self): bins = np.arange(0, 100, 10) - bin_centres, gamma = vario_estimate( + bin_centres, gamma = gs.vario_estimate( [x], field, bins, sampling_size=5000, sampling_seed=1479373475 ) @@ -181,7 +146,7 @@ def test_sampling_2d(self): bins = np.arange(0, 100, 10) - bin_centres, gamma = vario_estimate( + bin_centres, gamma = gs.vario_estimate( (x, y), field, bins, sampling_size=2000, sampling_seed=1479373475 ) @@ -204,7 +169,7 @@ def test_sampling_3d(self): bins = np.arange(0, 100, 10) - bin_centres, gamma = vario_estimate( + bin_centres, gamma = gs.vario_estimate( (x, y, z), field, bins, @@ -228,27 +193,21 @@ def test_assertions(self): field = np.arange(0, 10) field_e = np.arange(0, 9) + self.assertRaises(ValueError, gs.vario_estimate, [x_e], field, bins) + self.assertRaises(ValueError, gs.vario_estimate, (x, y_e), field, bins) self.assertRaises( - ValueError, vario_estimate, [x_e], field, bins - ) - self.assertRaises( - ValueError, vario_estimate, (x, y_e), field, bins - ) - self.assertRaises( - ValueError, vario_estimate, (x, y_e, z), field, bins + ValueError, gs.vario_estimate, (x, y_e, z), field, bins ) self.assertRaises( - ValueError, vario_estimate, (x, y, z_e), field, bins + ValueError, gs.vario_estimate, (x, y, z_e), field, bins ) self.assertRaises( - ValueError, vario_estimate, (x_e, y, z), field, bins + ValueError, gs.vario_estimate, (x_e, y, z), field, bins ) self.assertRaises( - ValueError, vario_estimate, (x, y, z), field_e, bins - ) - self.assertRaises( - ValueError, vario_estimate, [x], field_e, bins + ValueError, gs.vario_estimate, (x, y, z), field_e, bins ) + self.assertRaises(ValueError, gs.vario_estimate, [x], field_e, bins) def test_multi_field(self): x = np.random.RandomState(19970221).rand(100) * 100.0 @@ -259,9 +218,7 @@ def test_multi_field(self): bins = np.arange(20) * 2 bin_center, gamma1 = gs.vario_estimate(x, field1, bins) bin_center, gamma2 = gs.vario_estimate(x, field2, bins) - bin_center, gamma = gs.vario_estimate( - x, [field1, field2], bins - ) + bin_center, gamma = gs.vario_estimate(x, [field1, field2], bins) gamma_mean = 0.5 * (gamma1 + gamma2) for i in range(len(gamma)): self.assertAlmostEqual(gamma[i], gamma_mean[i], places=2) @@ -278,6 +235,28 @@ def test_no_data(self): for i in range(len(gamma1)): self.assertAlmostEqual(gamma1[i], gamma2[i], places=2) + def test_direction(self): + model = gs.Exponential(dim=3, len_scale=[12, 6, 3]) + x = y = z = range(10) + srf = gs.SRF(model, seed=123456) + field = np.ma.array(srf((x, y, z), mesh_type="structured")) + field.mask = np.abs(field) < 0.1 + + bins = range(10) + __, vario_u = gs.vario_estimate( + *((x, y, z), field, bins), + direction=((1, 0, 0), (0, 1, 0), (0, 0, 1)), # x-, y- and z-axis + bandwidth=0.25, # bandwith small enough to only match lines + mesh_type="structured", + ) + vario_s_x = gs.vario_estimate_axis(field, "x") + vario_s_y = gs.vario_estimate_axis(field, "y") + vario_s_z = gs.vario_estimate_axis(field, "z") + + self.assertTrue(np.all(np.isclose(vario_u[0], vario_s_x[:-1]))) + self.assertTrue(np.all(np.isclose(vario_u[1], vario_s_y[:-1]))) + self.assertTrue(np.all(np.isclose(vario_u[2], vario_s_z[:-1]))) + if __name__ == "__main__": unittest.main() From cb863848c6a221ba1b0c931a95c62df895c22902 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 9 Nov 2020 13:49:03 +0100 Subject: [PATCH 25/44] tests: add test for missing value in vario est. struct; comment out unused variables --- tests/test_variogram_structured.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/tests/test_variogram_structured.py b/tests/test_variogram_structured.py index 73f566c4..e93e499c 100644 --- a/tests/test_variogram_structured.py +++ b/tests/test_variogram_structured.py @@ -203,11 +203,11 @@ def test_directions_2d(self): # random values repeated along x-axis field_y = np.tile(y_rand, (len(x), 1)) - gamma_x_x = gs.vario_estimate_axis(field_x, direction="x") + # gamma_x_x = gs.vario_estimate_axis(field_x, direction="x") gamma_x_y = gs.vario_estimate_axis(field_x, direction="y") gamma_y_x = gs.vario_estimate_axis(field_y, direction="x") - gamma_y_y = gs.vario_estimate_axis(field_y, direction="y") + # gamma_y_y = gs.vario_estimate_axis(field_y, direction="y") self.assertAlmostEqual(gamma_x_y[1], 0.0) self.assertAlmostEqual(gamma_x_y[len(gamma_x_y) // 2], 0.0) @@ -229,17 +229,17 @@ def test_directions_3d(self): field_y = np.tile(y_rand.reshape((1, len(y), 1)), (len(x), 1, len(z))) field_z = np.tile(z_rand.reshape((1, 1, len(z))), (len(x), len(y), 1)) - gamma_x_x = gs.vario_estimate_axis(field_x, direction="x") + # gamma_x_x = gs.vario_estimate_axis(field_x, direction="x") gamma_x_y = gs.vario_estimate_axis(field_x, direction="y") gamma_x_z = gs.vario_estimate_axis(field_x, direction="z") gamma_y_x = gs.vario_estimate_axis(field_y, direction="x") - gamma_y_y = gs.vario_estimate_axis(field_y, direction="y") + # gamma_y_y = gs.vario_estimate_axis(field_y, direction="y") gamma_y_z = gs.vario_estimate_axis(field_y, direction="z") gamma_z_x = gs.vario_estimate_axis(field_z, direction="x") gamma_z_y = gs.vario_estimate_axis(field_z, direction="y") - gamma_z_z = gs.vario_estimate_axis(field_z, direction="z") + # gamma_z_z = gs.vario_estimate_axis(field_z, direction="z") self.assertAlmostEqual(gamma_x_y[1], 0.0) self.assertAlmostEqual(gamma_x_y[len(gamma_x_y) // 2], 0.0) @@ -262,10 +262,21 @@ def test_directions_3d(self): def test_exceptions(self): x = np.linspace(0.0, 10.0, 20) - rng = np.random.RandomState(1479373475) - x_rand = rng.rand(len(x)) + # rng = np.random.RandomState(1479373475) + # x_rand = rng.rand(len(x)) self.assertRaises(ValueError, gs.vario_estimate_axis, x, "a") + def test_missing(self): + x = np.linspace(0.0, 10.0, 10) + x_nan = x.copy() + x_nan[0] = np.nan + x_mask = np.isnan(x_nan) + x = np.ma.array(x, mask=x_mask) + v1 = gs.vario_estimate_axis(x_nan) + v2 = gs.vario_estimate_axis(x) + for i in range(len(v1)): + self.assertAlmostEqual(v1[i], v2[i]) + if __name__ == "__main__": unittest.main() From 04e2cdccda473168c9ff903510112445ac9216c5 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 9 Nov 2020 17:19:10 +0100 Subject: [PATCH 26/44] vario: handle numerical instabilities for vectors in same direction --- gstools/variogram/estimator.pyx | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 3028a3b1..9a3529dd 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -44,25 +44,31 @@ cdef inline bint dir_test( cdef double s_prod = 0.0 cdef double p_norm = 0.0 cdef double b_dist = 0.0 + cdef double tmp cdef int k cdef bint in_band = True - cdef bint in_angle + cdef bint in_angle = True # scalar-product calculation for bandwidth projection and angle calculation for k in range(dim): - s_prod += (pos[k,i] - pos[k,j]) * direction[d,k] - p_norm += (pos[k,i] - pos[k,j])**2 + tmp = (pos[k,i] - pos[k,j]) + s_prod += tmp * direction[d,k] + p_norm += tmp * tmp p_norm = sqrt(p_norm) # calculate band-distance by projection of point-pair-vec to direction line if bandwidth > 0.0: for k in range(dim): - b_dist += ((pos[k,i] - pos[k,j]) - s_prod * direction[d,k]) ** 2 - b_dist = sqrt(b_dist) - in_band = b_dist < bandwidth - - # use smallest angle by taking absolut value for arccos angle formula - in_angle = acos(fabs(s_prod) / p_norm) < angles_tol + tmp = (pos[k,i] - pos[k,j]) - s_prod * direction[d,k] + b_dist += tmp * tmp + in_band = sqrt(b_dist) < bandwidth + + # allow repeating points + if p_norm > 0.0: + # use smallest angle by taking absolut value for arccos angle formula + tmp = fabs(s_prod) / p_norm + if tmp < 1.0: # else same direction (prevent numerical errors) + in_angle = acos(tmp) < angles_tol return in_band and in_angle From 5ac4dd8489e443288e174a041bdcf6b48e78c95f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 9 Nov 2020 17:19:57 +0100 Subject: [PATCH 27/44] tests: testing directional variograms no.1 --- tests/test_variogram_unstructured.py | 42 ++++++++++++++++++++-------- 1 file changed, 30 insertions(+), 12 deletions(-) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 60e037a1..bdd2b8e7 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -10,7 +10,11 @@ class TestVariogramUnstructured(unittest.TestCase): def setUp(self): - pass + model = gs.Exponential(dim=3, len_scale=[12, 6, 3]) + x = y = z = range(10) + self.pos = (x, y, z) + srf = gs.SRF(model, seed=123456) + self.field = srf((x, y, z), mesh_type="structured") def test_doubles(self): x = np.arange(1, 11, 1, dtype=np.double) @@ -235,16 +239,12 @@ def test_no_data(self): for i in range(len(gamma1)): self.assertAlmostEqual(gamma1[i], gamma2[i], places=2) - def test_direction(self): - model = gs.Exponential(dim=3, len_scale=[12, 6, 3]) - x = y = z = range(10) - srf = gs.SRF(model, seed=123456) - field = np.ma.array(srf((x, y, z), mesh_type="structured")) + def test_direction_axis(self): + field = np.ma.array(self.field) field.mask = np.abs(field) < 0.1 - bins = range(10) __, vario_u = gs.vario_estimate( - *((x, y, z), field, bins), + *(self.pos, field, bins), direction=((1, 0, 0), (0, 1, 0), (0, 0, 1)), # x-, y- and z-axis bandwidth=0.25, # bandwith small enough to only match lines mesh_type="structured", @@ -252,10 +252,28 @@ def test_direction(self): vario_s_x = gs.vario_estimate_axis(field, "x") vario_s_y = gs.vario_estimate_axis(field, "y") vario_s_z = gs.vario_estimate_axis(field, "z") - - self.assertTrue(np.all(np.isclose(vario_u[0], vario_s_x[:-1]))) - self.assertTrue(np.all(np.isclose(vario_u[1], vario_s_y[:-1]))) - self.assertTrue(np.all(np.isclose(vario_u[2], vario_s_z[:-1]))) + for i in range(len(bins) - 1): + self.assertAlmostEqual(vario_u[0][i], vario_s_x[i]) + self.assertAlmostEqual(vario_u[1][i], vario_s_y[i]) + self.assertAlmostEqual(vario_u[2][i], vario_s_z[i]) + + def test_direction_angle(self): + bins = range(0, 10, 2) + __, v2, c2 = gs.vario_estimate( + *(self.pos[:2], self.field[0], bins), + angles=np.pi / 4, # 45 deg + mesh_type="structured", + return_counts=True, + ) + __, v1, c1 = gs.vario_estimate( + *(self.pos[:2], self.field[0], bins), + direction=(1, 1), # 45 deg + mesh_type="structured", + return_counts=True, + ) + for i in range(len(bins) - 1): + self.assertAlmostEqual(v1[i], v2[i]) + self.assertAlmostEqual(c1[i], c2[i]) if __name__ == "__main__": From afb18d8aa19228ab7e9333cf507514b15ba2f07c Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 9 Nov 2020 22:30:45 +0100 Subject: [PATCH 28/44] vario: re-use calculated distance in for direction check in cython --- gstools/variogram/estimator.pyx | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 9a3529dd..54ad039b 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -34,6 +34,7 @@ cdef inline double distance( cdef inline bint dir_test( const int dim, const double[:,:] pos, + const double dist, const double[:,:] direction, const double angles_tol, const double bandwidth, @@ -41,20 +42,16 @@ cdef inline bint dir_test( const int j, const int d ) nogil: - cdef double s_prod = 0.0 - cdef double p_norm = 0.0 - cdef double b_dist = 0.0 - cdef double tmp + cdef double s_prod = 0.0 # scalar product + cdef double b_dist = 0.0 # band-distance + cdef double tmp # temporary variable cdef int k cdef bint in_band = True cdef bint in_angle = True # scalar-product calculation for bandwidth projection and angle calculation for k in range(dim): - tmp = (pos[k,i] - pos[k,j]) - s_prod += tmp * direction[d,k] - p_norm += tmp * tmp - p_norm = sqrt(p_norm) + s_prod += (pos[k,i] - pos[k,j]) * direction[d,k] # calculate band-distance by projection of point-pair-vec to direction line if bandwidth > 0.0: @@ -63,10 +60,10 @@ cdef inline bint dir_test( b_dist += tmp * tmp in_band = sqrt(b_dist) < bandwidth - # allow repeating points - if p_norm > 0.0: + # allow repeating points (dist = 0) + if dist > 0.0: # use smallest angle by taking absolut value for arccos angle formula - tmp = fabs(s_prod) / p_norm + tmp = fabs(s_prod) / dist if tmp < 1.0: # else same direction (prevent numerical errors) in_angle = acos(tmp) < angles_tol @@ -208,7 +205,7 @@ def directional( dist = distance(dim, pos, j, k) if dist >= bin_edges[i] and dist < bin_edges[i+1]: for d in range(d_max): - if dir_test(dim, pos, direction, angles_tol, bandwidth, k, j, d): + if dir_test(dim, pos, dist, direction, angles_tol, bandwidth, k, j, d): for m in range(f_max): # skip no data values if not (isnan(f[m,k]) or isnan(f[m,j])): From 01bc57c53de102f87de8341752e57a9dc7ebfaf3 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 9 Nov 2020 23:23:18 +0100 Subject: [PATCH 29/44] vario: better shape handling for angles and mask --- gstools/tools/geometric.py | 3 ++- gstools/variogram/variogram.py | 4 +++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index dfd5dd03..ed737eb4 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -185,11 +185,12 @@ def ang2dir(angles, dtype=np.double, dim=None): :class:`numpy.ndarray` the list of direction vectors """ + pre_dim = np.asanyarray(angles).ndim angles = np.array(angles, ndmin=2, dtype=dtype) if len(angles.shape) > 2: raise ValueError("Can't interpret angles array {}".format(angles)) dim = angles.shape[1] + 1 if dim is None else dim - if dim == 2 and angles.shape[0] == 1: + if dim == 2 and angles.shape[0] == 1 and pre_dim < 2: # fix for 2D where only one angle per direction is given angles = angles.T # can't be interpreted if dim=None is given if dim != angles.shape[1] + 1 or dim == 1: diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index d27ed661..e8881f59 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -209,7 +209,9 @@ def vario_estimate( # selected region is the inverted masked (unmasked values) if np.size(mask) > 1: # not only np.ma.nomask select = np.invert( - np.logical_or(mask.reshape(len(x)), np.all(field.mask, axis=0)) + np.logical_or( + np.reshape(mask, len(x)), np.all(field.mask, axis=0) + ) ) else: select = np.invert(np.all(field.mask, axis=0)) From 3e3a07f2cc5e8abc1f3aa1e1eb1f4ceafd15ce0d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 9 Nov 2020 23:23:51 +0100 Subject: [PATCH 30/44] test: vario assertions and no_data checks added --- tests/test_variogram_unstructured.py | 49 ++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index bdd2b8e7..d75563e5 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -275,6 +275,55 @@ def test_direction_angle(self): self.assertAlmostEqual(v1[i], v2[i]) self.assertAlmostEqual(c1[i], c2[i]) + def test_direction_assertion(self): + pos = [[1, 2, 3], [1, 2, 3]] + bns = [1, 2] + fld = np.ma.array([1, 2, 3]) + self.assertRaises( # degenerated direction + ValueError, gs.vario_estimate, pos, fld, bns, direction=[0, 0] + ) + self.assertRaises( # wrong shape of direction + ValueError, gs.vario_estimate, pos, fld, bns, direction=[[[3, 1]]] + ) + self.assertRaises( # wrong dimension of direction + ValueError, gs.vario_estimate, pos, fld, bns, direction=[[3, 1, 2]] + ) + self.assertRaises( # wrong shape of angles + ValueError, gs.vario_estimate, pos, fld, bns, angles=[[[1]]] + ) + self.assertRaises( # wrong dimension of angles + ValueError, gs.vario_estimate, pos, fld, bns, angles=[[1, 1]] + ) + + def test_mask_no_data(self): + pos = [[1, 2, 3, 4, 5], [1, 2, 3, 4, 5]] + bns = [0, 4] + fld1 = np.ma.array([1, 2, 3, 4, 5]) + fld2 = np.ma.array([np.nan, 2, 3, 4, 5]) + fld3 = np.ma.array([1, 2, 3, 4, 5]) + mask = [False, False, True, False, False] + fld1.mask = [True, False, False, False, False] + fld2.mask = mask + __, v1, c1 = gs.vario_estimate( + *(pos, fld1, bns), + mask=mask, + return_counts=True, + ) + __, v2, c2 = gs.vario_estimate( + *(pos, fld2, bns), + return_counts=True, + ) + __, v3, c3 = gs.vario_estimate( + *(pos, fld3, bns), + no_data=1, + mask=mask, + return_counts=True, + ) + self.assertAlmostEqual(v1[0], v2[0]) + self.assertAlmostEqual(v1[0], v3[0]) + self.assertAlmostEqual(c1[0], c2[0]) + self.assertAlmostEqual(c1[0], c3[0]) + if __name__ == "__main__": unittest.main() From 89e0e93372d52aec5226abaf18583db8e68c9a03 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Mon, 9 Nov 2020 23:56:02 +0100 Subject: [PATCH 31/44] vario: no valueError for full masked arrays; test wrong estimator Error --- gstools/variogram/variogram.py | 12 ++++++++---- tests/test_variogram_unstructured.py | 3 +++ 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index e8881f59..02c2af81 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -181,18 +181,22 @@ def vario_estimate( ----- Internally uses double precision and also returns doubles. """ + bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double) + 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 field = np.ma.array(field, ndmin=2, dtype=np.double) masked = np.ma.is_masked(field) or np.any(mask) + # catch special case if everything is masked if masked and np.all(mask): - raise ValueError("Given mask, masks all values.") + estimates = np.zeros_like(bin_centres) + if return_counts: + return bin_centres, estimates, np.zeros_like(estimates, dtype=int) + return bin_centres, estimates if not masked: field = field.filled() - bin_edges = np.array(bin_edges, ndmin=1, dtype=np.double) - x, y, z, dim = pos2xyz(pos, calc_dim=True, dtype=np.double) - bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2.0 # check_mesh shape + x, y, z, dim = pos2xyz(pos, calc_dim=True, dtype=np.double) if mesh_type != "unstructured": x, y, z, __ = reshape_axis_from_struct_to_unstruct(dim, x, y, z) if len(field.shape) > 2 or field.shape[1] != len(x): diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index d75563e5..8a02076b 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -212,6 +212,9 @@ def test_assertions(self): ValueError, gs.vario_estimate, (x, y, z), field_e, bins ) self.assertRaises(ValueError, gs.vario_estimate, [x], field_e, bins) + self.assertRaises( + ValueError, gs.vario_estimate, [x], field, bins, estimator="bla" + ) def test_multi_field(self): x = np.random.RandomState(19970221).rand(100) * 100.0 From ddbf1480c19531a29f0604cf72de240c758b0d73 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 00:11:23 +0100 Subject: [PATCH 32/44] test: test fully masked vario estimate --- tests/test_variogram_unstructured.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 8a02076b..5af6e709 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -276,7 +276,7 @@ def test_direction_angle(self): ) for i in range(len(bins) - 1): self.assertAlmostEqual(v1[i], v2[i]) - self.assertAlmostEqual(c1[i], c2[i]) + self.assertEqual(c1[i], c2[i]) def test_direction_assertion(self): pos = [[1, 2, 3], [1, 2, 3]] @@ -322,10 +322,17 @@ def test_mask_no_data(self): mask=mask, return_counts=True, ) + __, v4, c4 = gs.vario_estimate( + *(pos, fld3, bns), + mask=True, + return_counts=True, + ) self.assertAlmostEqual(v1[0], v2[0]) self.assertAlmostEqual(v1[0], v3[0]) - self.assertAlmostEqual(c1[0], c2[0]) - self.assertAlmostEqual(c1[0], c3[0]) + self.assertEqual(c1[0], c2[0]) + self.assertEqual(c1[0], c3[0]) + self.assertAlmostEqual(v4[0], 0.0) + self.assertEqual(c4[0], 0) if __name__ == "__main__": From f0e9c610b0e3f594db1c56b488d199c9a55c950c Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 00:55:55 +0100 Subject: [PATCH 33/44] examples: 2d example for directional variogram --- examples/03_variogram/04_directional_2d.py | 48 ++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100755 examples/03_variogram/04_directional_2d.py diff --git a/examples/03_variogram/04_directional_2d.py b/examples/03_variogram/04_directional_2d.py new file mode 100755 index 00000000..fb74ba82 --- /dev/null +++ b/examples/03_variogram/04_directional_2d.py @@ -0,0 +1,48 @@ +""" +Directional variogram estimation in 2D +-------------------------------------- + +In this example, we demonstrate how to estimate a directional variogram by +setting the direction angles in 2D. +""" +import numpy as np +import gstools as gs +from matplotlib import pyplot as plt + + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 5]) + +############################################################################### +# Generating synthetic field with anisotropy and a rotation of 22.5 degree. + +angle = np.pi / 8 +model = gs.Exponential(dim=2, len_scale=[10, 5], angles=angle) +x = y = range(100) +srf = gs.SRF(model, seed=123456) +field = srf((x, y), mesh_type="structured") +srf.plot(ax=ax2) +ax2.set_aspect("equal") + +############################################################################### +# Now we are going to estimate a directional variogram with an angular +# tolerance of 11.25 degree and a bandwith of 1. +# We provide the rotation angle of the covariance model and the orthogonal +# direction by adding 90 degree. + +bins = range(0, 40, 3) +bin_c, vario, cnt = gs.vario_estimate_unstructured( + *((x, y), field, bins), + angles=(angle, angle + np.pi / 2), + angles_tol=np.pi / 16, + bandwidth=1.0, + mesh_type="structured", + return_counts=True, +) +ax1.plot(bin_c, vario[0], label="emp. vario: pi/8") +ax1.plot(bin_c, vario[1], label="emp. vario: pi*5/8") +ax1.legend(loc="lower right") +fig.show() + +############################################################################### +# Without fitting a model we see, that the correlation length in the main +# direction is greater than the transversal one. From d834c0a5f1911742099f42113c945e8d51a34a8f Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 00:57:47 +0100 Subject: [PATCH 34/44] examples: remove lagacy call --- examples/03_variogram/04_directional_2d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/03_variogram/04_directional_2d.py b/examples/03_variogram/04_directional_2d.py index fb74ba82..30b8d210 100755 --- a/examples/03_variogram/04_directional_2d.py +++ b/examples/03_variogram/04_directional_2d.py @@ -30,7 +30,7 @@ # direction by adding 90 degree. bins = range(0, 40, 3) -bin_c, vario, cnt = gs.vario_estimate_unstructured( +bin_c, vario, cnt = gs.vario_estimate( *((x, y), field, bins), angles=(angle, angle + np.pi / 2), angles_tol=np.pi / 16, From 0c323ec6d745ffe554f33c1fbcc74de21c4b0fc5 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 01:14:49 +0100 Subject: [PATCH 35/44] examples: 3d example for directional variogram --- examples/03_variogram/04_directional_2d.py | 2 +- examples/03_variogram/05_directional_3d.py | 73 ++++++++++++++++++++++ 2 files changed, 74 insertions(+), 1 deletion(-) create mode 100755 examples/03_variogram/05_directional_3d.py diff --git a/examples/03_variogram/04_directional_2d.py b/examples/03_variogram/04_directional_2d.py index 30b8d210..f13fc452 100755 --- a/examples/03_variogram/04_directional_2d.py +++ b/examples/03_variogram/04_directional_2d.py @@ -32,7 +32,7 @@ bins = range(0, 40, 3) bin_c, vario, cnt = gs.vario_estimate( *((x, y), field, bins), - angles=(angle, angle + np.pi / 2), + angles=(angle, angle + np.pi / 2), # main dir. and transversal dir. angles_tol=np.pi / 16, bandwidth=1.0, mesh_type="structured", diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py new file mode 100755 index 00000000..b437edab --- /dev/null +++ b/examples/03_variogram/05_directional_3d.py @@ -0,0 +1,73 @@ +""" +Directional variogram estimation in 2D +-------------------------------------- + +In this example, we demonstrate how to estimate a directional variogram by +setting the estimation directions in 3D. +""" +import numpy as np +import gstools as gs +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + + +fig = plt.figure(figsize=[15, 5]) +ax1 = fig.add_subplot(131, projection=Axes3D.name) +ax2 = fig.add_subplot(132, projection=Axes3D.name) +ax3 = fig.add_subplot(133) + +############################################################################### +# Generating synthetic field with anisotropy and rotation by Tait-Bryan angles. + +dim = 3 +# rotation around z, y, x +angles = [np.pi / 2, np.pi / 4, np.pi / 8] +model = gs.Gaussian(dim=3, len_scale=[16, 8, 4], angles=angles) +x = y = z = range(50) +srf = gs.SRF(model, seed=1001) +field = srf.structured((x, y, z)) +srf.plot(ax=ax1) +ax1.set_aspect("equal") + +############################################################################### +# Here we plot the rotated coordinate system to get an impression, what +# the rotation angles do. + +x1, x2, x3 = (1, 0, 0), (0, 1, 0), (0, 0, 1) +ret = np.array(gs.field.tools.rotate_mesh(dim, angles, x1, x2, x3)) +dir0 = ret[:, 0] # main direction +dir1 = ret[:, 1] # first lateral direction +dir2 = ret[:, 2] # second lateral direction +ax2.plot([0, dir0[0]], [0, dir0[1]], [0, dir0[2]], label="1.") +ax2.plot([0, dir1[0]], [0, dir1[1]], [0, dir1[2]], label="2.") +ax2.plot([0, dir2[0]], [0, dir2[1]], [0, dir2[2]], label="3.") +ax2.set_xlim(-1, 1) +ax2.set_ylim(-1, 1) +ax2.set_zlim(-1, 1) +ax2.set_xlabel("X") +ax2.set_ylabel("Y") +ax2.set_zlabel("Z") +ax2.set_aspect("equal") +ax2.set_title("Tait-Bryan main axis") +ax2.legend(loc="lower left") + +############################################################################### +# Now we estimate the variogram along the main axis. When the main axis are +# unknown, one would need to sample multiple directions and look out for +# the one with the longest correlation length (flattest gradient). +# Then check the transversal directions and so on. + +bins = range(0, 40, 3) +bin_c, vario = gs.vario_estimate( + *([x, y, z], field, bins), + direction=(dir0, dir1, dir2), + bandwidth=10, + sampling_size=2000, + sampling_seed=1001, + mesh_type="structured" +) +ax3.plot(bin_c, vario[0], label="1. axis") +ax3.plot(bin_c, vario[1], label="2. axis") +ax3.plot(bin_c, vario[2], label="3. axis") +ax3.legend() +fig.show() From 9e88943c8c8fefc9fd0157503bc68aa17048bb73 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 02:33:52 +0100 Subject: [PATCH 36/44] tests: vario finally at 100% --- tests/test_variogram_unstructured.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_variogram_unstructured.py b/tests/test_variogram_unstructured.py index 5af6e709..8e4e7ee0 100644 --- a/tests/test_variogram_unstructured.py +++ b/tests/test_variogram_unstructured.py @@ -312,10 +312,7 @@ def test_mask_no_data(self): mask=mask, return_counts=True, ) - __, v2, c2 = gs.vario_estimate( - *(pos, fld2, bns), - return_counts=True, - ) + __, v2, c2 = gs.vario_estimate(*(pos, fld2, bns), return_counts=True) __, v3, c3 = gs.vario_estimate( *(pos, fld3, bns), no_data=1, @@ -327,12 +324,15 @@ def test_mask_no_data(self): mask=True, return_counts=True, ) + __, v5 = gs.vario_estimate(*(pos, fld3, bns), mask=True) + self.assertAlmostEqual(v1[0], v2[0]) self.assertAlmostEqual(v1[0], v3[0]) self.assertEqual(c1[0], c2[0]) self.assertEqual(c1[0], c3[0]) self.assertAlmostEqual(v4[0], 0.0) self.assertEqual(c4[0], 0) + self.assertAlmostEqual(v5[0], 0.0) if __name__ == "__main__": From 539aff414386781dfa11e8166fbfa1a30f89a729 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 02:35:08 +0100 Subject: [PATCH 37/44] BUGFIX: 3D contourf plots in gstools not working with mpl 3.3 -> use 2D --- gstools/field/plot.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/gstools/field/plot.py b/gstools/field/plot.py index 23bcf3b8..73f5b7e8 100644 --- a/gstools/field/plot.py +++ b/gstools/field/plot.py @@ -15,7 +15,6 @@ from scipy import interpolate as inter import matplotlib.pyplot as plt from matplotlib.widgets import Slider, RadioButtons -from mpl_toolkits.mplot3d import Axes3D from gstools.tools import pos2xyz from gstools.covmodel.plot import _get_fig_ax @@ -91,7 +90,7 @@ def _plot_2d(pos, field, mesh_type, fig=None, ax=None): # pragma: no cover def _plot_3d(pos, field, mesh_type, fig=None, ax=None): # pragma: no cover """Plot 3D field.""" dir1, dir2 = np.mgrid[0:1:51j, 0:1:51j] - levels = np.linspace(field.min(), field.max(), 256, endpoint=True) + levels = np.linspace(field.min(), field.max(), 100, endpoint=True) x_min = pos[0].min() x_max = pos[0].max() @@ -110,7 +109,7 @@ def _plot_3d(pos, field, mesh_type, fig=None, ax=None): # pragma: no cover "y": [y_min, y_max, y_range, y_step], "z": [z_min, z_max, z_range, z_step], } - fig, ax = _get_fig_ax(fig, ax, Axes3D.name) + fig, ax = _get_fig_ax(fig, ax) title = "Field 3D " + mesh_type + ": " + str(field.shape) fig.subplots_adjust(left=0.2, right=0.8, bottom=0.25) sax = plt.axes([0.15, 0.1, 0.65, 0.03]) @@ -122,7 +121,7 @@ def _plot_3d(pos, field, mesh_type, fig=None, ax=None): # pragma: no cover valinit=z_min + z_range / 2.0, valstep=z_step, ) - rax = plt.axes([0.05, 0.5, 0.1, 0.15]) + rax = plt.axes([0.05, 0.7, 0.1, 0.15]) radio = RadioButtons(rax, ("x slice", "y slice", "z slice"), active=2) z_dir_tmp = "z" # create container @@ -155,13 +154,11 @@ def get_plane(z_val_in, z_dir): plane = inter.griddata( pos, field, (x_io, y_io, z_io), method="linear" ) - if z_dir == "z": - z_io = plane + if z_dir == "x": + return y_io, z_io, plane elif z_dir == "y": - y_io = plane - else: - x_io = plane - return x_io, y_io, z_io + return x_io, z_io, plane + return x_io, y_io, plane def update(__): """Widget update.""" @@ -188,17 +185,20 @@ def update(__): vmin=field.min(), vmax=field.max(), levels=levels, - zdir=z_dir_in, - offset=z_val, ) - cont.cmap.set_under("k", alpha=0.0) - cont.cmap.set_bad("k", alpha=0.0) - ax.set_xlabel("X") - ax.set_ylabel("Y") - ax.set_zlabel("Z") + # cont.cmap.set_under("k", alpha=0.0) + # cont.cmap.set_bad("k", alpha=0.0) + if z_dir_in == "x": + ax.set_xlabel("Y") + ax.set_ylabel("Z") + elif z_dir_in == "y": + ax.set_xlabel("X") + ax.set_ylabel("Z") + else: + ax.set_xlabel("X") + ax.set_ylabel("Y") ax.set_xlim([x_min, x_max]) ax.set_ylim([y_min, y_max]) - ax.set_zlim([z_min, z_max]) ax.set_title(title) fig.canvas.draw_idle() return cont From 849d49c85366c181c9ee9cd3a745f16781e1f49d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 02:35:30 +0100 Subject: [PATCH 38/44] examples: switch to 2D plot for 3D field --- examples/03_variogram/05_directional_3d.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py index b437edab..245a291b 100755 --- a/examples/03_variogram/05_directional_3d.py +++ b/examples/03_variogram/05_directional_3d.py @@ -12,7 +12,7 @@ fig = plt.figure(figsize=[15, 5]) -ax1 = fig.add_subplot(131, projection=Axes3D.name) +ax1 = fig.add_subplot(131) ax2 = fig.add_subplot(132, projection=Axes3D.name) ax3 = fig.add_subplot(133) @@ -47,7 +47,6 @@ ax2.set_xlabel("X") ax2.set_ylabel("Y") ax2.set_zlabel("Z") -ax2.set_aspect("equal") ax2.set_title("Tait-Bryan main axis") ax2.legend(loc="lower left") From f22a7c311ad5b7b4df3e4e6482457ca2dc2e9666 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 03:00:46 +0100 Subject: [PATCH 39/44] examples: use plt.show instead of fig.show --- examples/00_misc/02_check_rand_meth_sampling.py | 2 +- examples/03_variogram/04_directional_2d.py | 2 +- examples/03_variogram/05_directional_3d.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/00_misc/02_check_rand_meth_sampling.py b/examples/00_misc/02_check_rand_meth_sampling.py index c03213c6..e89c1340 100644 --- a/examples/00_misc/02_check_rand_meth_sampling.py +++ b/examples/00_misc/02_check_rand_meth_sampling.py @@ -58,7 +58,7 @@ def plot_rand_meth_samples(generator): ax.set_xlim([0, np.max(x)]) ax.set_title("Radius samples shown {}/{}".format(sample_in, len(rad))) ax.legend() - fig.show() + plt.show() model = gs.Stable(dim=3, alpha=1.5) diff --git a/examples/03_variogram/04_directional_2d.py b/examples/03_variogram/04_directional_2d.py index f13fc452..5373ee60 100755 --- a/examples/03_variogram/04_directional_2d.py +++ b/examples/03_variogram/04_directional_2d.py @@ -41,7 +41,7 @@ ax1.plot(bin_c, vario[0], label="emp. vario: pi/8") ax1.plot(bin_c, vario[1], label="emp. vario: pi*5/8") ax1.legend(loc="lower right") -fig.show() +plt.show() ############################################################################### # Without fitting a model we see, that the correlation length in the main diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py index 245a291b..a8e13a5f 100755 --- a/examples/03_variogram/05_directional_3d.py +++ b/examples/03_variogram/05_directional_3d.py @@ -69,4 +69,4 @@ ax3.plot(bin_c, vario[1], label="2. axis") ax3.plot(bin_c, vario[2], label="3. axis") ax3.legend() -fig.show() +plt.show() From 9302564b8a9450fff533e054359266c4d89f73e2 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 03:14:48 +0100 Subject: [PATCH 40/44] examples: fix plotting order for sphinx gallery --- examples/03_variogram/04_directional_2d.py | 15 +++++--- examples/03_variogram/05_directional_3d.py | 45 ++++++++++++---------- 2 files changed, 35 insertions(+), 25 deletions(-) diff --git a/examples/03_variogram/04_directional_2d.py b/examples/03_variogram/04_directional_2d.py index 5373ee60..f60bfdd4 100755 --- a/examples/03_variogram/04_directional_2d.py +++ b/examples/03_variogram/04_directional_2d.py @@ -9,9 +9,6 @@ import gstools as gs from matplotlib import pyplot as plt - -fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 5]) - ############################################################################### # Generating synthetic field with anisotropy and a rotation of 22.5 degree. @@ -20,8 +17,6 @@ x = y = range(100) srf = gs.SRF(model, seed=123456) field = srf((x, y), mesh_type="structured") -srf.plot(ax=ax2) -ax2.set_aspect("equal") ############################################################################### # Now we are going to estimate a directional variogram with an angular @@ -38,9 +33,19 @@ mesh_type="structured", return_counts=True, ) + +############################################################################### +# Plotting. + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[10, 5]) + ax1.plot(bin_c, vario[0], label="emp. vario: pi/8") ax1.plot(bin_c, vario[1], label="emp. vario: pi*5/8") ax1.legend(loc="lower right") + +srf.plot(ax=ax2) +ax2.set_aspect("equal") + plt.show() ############################################################################### diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py index a8e13a5f..1bd4e2d7 100755 --- a/examples/03_variogram/05_directional_3d.py +++ b/examples/03_variogram/05_directional_3d.py @@ -10,12 +10,6 @@ import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D - -fig = plt.figure(figsize=[15, 5]) -ax1 = fig.add_subplot(131) -ax2 = fig.add_subplot(132, projection=Axes3D.name) -ax3 = fig.add_subplot(133) - ############################################################################### # Generating synthetic field with anisotropy and rotation by Tait-Bryan angles. @@ -26,11 +20,9 @@ x = y = z = range(50) srf = gs.SRF(model, seed=1001) field = srf.structured((x, y, z)) -srf.plot(ax=ax1) -ax1.set_aspect("equal") ############################################################################### -# Here we plot the rotated coordinate system to get an impression, what +# Here we generate the rotated coordinate system to get an impression, what # the rotation angles do. x1, x2, x3 = (1, 0, 0), (0, 1, 0), (0, 0, 1) @@ -38,17 +30,6 @@ dir0 = ret[:, 0] # main direction dir1 = ret[:, 1] # first lateral direction dir2 = ret[:, 2] # second lateral direction -ax2.plot([0, dir0[0]], [0, dir0[1]], [0, dir0[2]], label="1.") -ax2.plot([0, dir1[0]], [0, dir1[1]], [0, dir1[2]], label="2.") -ax2.plot([0, dir2[0]], [0, dir2[1]], [0, dir2[2]], label="3.") -ax2.set_xlim(-1, 1) -ax2.set_ylim(-1, 1) -ax2.set_zlim(-1, 1) -ax2.set_xlabel("X") -ax2.set_ylabel("Y") -ax2.set_zlabel("Z") -ax2.set_title("Tait-Bryan main axis") -ax2.legend(loc="lower left") ############################################################################### # Now we estimate the variogram along the main axis. When the main axis are @@ -65,6 +46,30 @@ sampling_seed=1001, mesh_type="structured" ) + +############################################################################### +# Plotting. + +fig = plt.figure(figsize=[15, 5]) +ax1 = fig.add_subplot(131) +ax2 = fig.add_subplot(132, projection=Axes3D.name) +ax3 = fig.add_subplot(133) + +srf.plot(ax=ax1) +ax1.set_aspect("equal") + +ax2.plot([0, dir0[0]], [0, dir0[1]], [0, dir0[2]], label="1.") +ax2.plot([0, dir1[0]], [0, dir1[1]], [0, dir1[2]], label="2.") +ax2.plot([0, dir2[0]], [0, dir2[1]], [0, dir2[2]], label="3.") +ax2.set_xlim(-1, 1) +ax2.set_ylim(-1, 1) +ax2.set_zlim(-1, 1) +ax2.set_xlabel("X") +ax2.set_ylabel("Y") +ax2.set_zlabel("Z") +ax2.set_title("Tait-Bryan main axis") +ax2.legend(loc="lower left") + ax3.plot(bin_c, vario[0], label="1. axis") ax3.plot(bin_c, vario[1], label="2. axis") ax3.plot(bin_c, vario[2], label="3. axis") From 1fadb670c0e7c2b2ccd5b10a040011a52da1d1fb Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 10 Nov 2020 03:25:43 +0100 Subject: [PATCH 41/44] examples: typo --- examples/03_variogram/05_directional_3d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py index 1bd4e2d7..358aa43a 100755 --- a/examples/03_variogram/05_directional_3d.py +++ b/examples/03_variogram/05_directional_3d.py @@ -1,5 +1,5 @@ """ -Directional variogram estimation in 2D +Directional variogram estimation in 3D -------------------------------------- In this example, we demonstrate how to estimate a directional variogram by From 80f9f3256288df14bad5ce8470c0a05d20deea0c Mon Sep 17 00:00:00 2001 From: LSchueler Date: Tue, 10 Nov 2020 20:13:22 +0100 Subject: [PATCH 42/44] Fix some very minor typos in comments --- examples/03_variogram/03_multi_vario.py | 4 ++-- examples/03_variogram/04_directional_2d.py | 2 +- examples/03_variogram/05_directional_3d.py | 8 ++++---- gstools/tools/geometric.py | 2 +- gstools/variogram/estimator.pyx | 4 ++-- 5 files changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/03_variogram/03_multi_vario.py b/examples/03_variogram/03_multi_vario.py index 08951c08..62eb3d5b 100755 --- a/examples/03_variogram/03_multi_vario.py +++ b/examples/03_variogram/03_multi_vario.py @@ -3,7 +3,7 @@ -------------------------------- In this example, we demonstrate how to estimate a variogram from multiple -fields at the same point-set, that should have the same statistical properties. +fields on the same point-set that should have the same statistical properties. """ import numpy as np import gstools as gs @@ -32,7 +32,7 @@ bin_center, gamma = gs.vario_estimate((x, y), fields, bins) ############################################################################### -# Now we demonstrate, that the mean variogram from both fields coincides +# Now we demonstrate that the mean variogram from both fields coincides # with the joined estimated one. plt.plot(bin_center, gamma1, label="field 1") diff --git a/examples/03_variogram/04_directional_2d.py b/examples/03_variogram/04_directional_2d.py index f60bfdd4..b1a42920 100755 --- a/examples/03_variogram/04_directional_2d.py +++ b/examples/03_variogram/04_directional_2d.py @@ -49,5 +49,5 @@ plt.show() ############################################################################### -# Without fitting a model we see, that the correlation length in the main +# Without fitting a model, we see that the correlation length in the main # direction is greater than the transversal one. diff --git a/examples/03_variogram/05_directional_3d.py b/examples/03_variogram/05_directional_3d.py index 358aa43a..521cacf9 100755 --- a/examples/03_variogram/05_directional_3d.py +++ b/examples/03_variogram/05_directional_3d.py @@ -22,7 +22,7 @@ field = srf.structured((x, y, z)) ############################################################################### -# Here we generate the rotated coordinate system to get an impression, what +# Here we generate the rotated coordinate system to get an impression what # the rotation angles do. x1, x2, x3 = (1, 0, 0), (0, 1, 0), (0, 0, 1) @@ -32,9 +32,9 @@ dir2 = ret[:, 2] # second lateral direction ############################################################################### -# Now we estimate the variogram along the main axis. When the main axis are -# unknown, one would need to sample multiple directions and look out for -# the one with the longest correlation length (flattest gradient). +# Now we estimate the variogram along the main axis. When the main axis is +# unknown, one would need to sample multiple directions and look for the one +# with the longest correlation length (flattest gradient). # Then check the transversal directions and so on. bins = range(0, 40, 3) diff --git a/gstools/tools/geometric.py b/gstools/tools/geometric.py index ed737eb4..dc34edba 100644 --- a/gstools/tools/geometric.py +++ b/gstools/tools/geometric.py @@ -183,7 +183,7 @@ def ang2dir(angles, dtype=np.double, dim=None): Returns ------- :class:`numpy.ndarray` - the list of direction vectors + the array of direction vectors """ pre_dim = np.asanyarray(angles).ndim angles = np.array(angles, ndmin=2, dtype=dtype) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 54ad039b..abe080d7 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -62,7 +62,7 @@ cdef inline bint dir_test( # allow repeating points (dist = 0) if dist > 0.0: - # use smallest angle by taking absolut value for arccos angle formula + # use smallest angle by taking absolute value for arccos angle formula tmp = fabs(s_prod) / dist if tmp < 1.0: # else same direction (prevent numerical errors) in_angle = acos(tmp) < angles_tol @@ -168,7 +168,7 @@ def directional( const double[:,:] pos, const double[:,:] direction, # should be normed const double angles_tol=M_PI/8.0, - const double bandwidth=-1.0, # negative values to turn of bandwidth search + const double bandwidth=-1.0, # neg. values to turn off bandwidth search str estimator_type='m' ): if pos.shape[1] != f.shape[1]: From 2f4b7f1805429b3e757bdcfe5e6dff8b333e2da5 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 14:54:53 +0100 Subject: [PATCH 43/44] Vario: add separate_dirs argument to directional variogram est routines in cython; some code refactoring --- gstools/variogram/estimator.pyx | 36 ++++++++++++++++++++------------- 1 file changed, 22 insertions(+), 14 deletions(-) diff --git a/gstools/variogram/estimator.pyx b/gstools/variogram/estimator.pyx index 54ad039b..0d8a48f3 100644 --- a/gstools/variogram/estimator.pyx +++ b/gstools/variogram/estimator.pyx @@ -169,6 +169,7 @@ def directional( const double[:,:] direction, # should be normed const double angles_tol=M_PI/8.0, const double bandwidth=-1.0, # negative values to turn of bandwidth search + const bint separate_dirs=False, # whether the direction bands don't overlap str estimator_type='m' ): if pos.shape[1] != f.shape[1]: @@ -203,14 +204,20 @@ def directional( for j in range(j_max): for k in range(j+1, k_max): dist = distance(dim, pos, j, k) - if dist >= bin_edges[i] and dist < bin_edges[i+1]: - for d in range(d_max): - if dir_test(dim, pos, dist, direction, angles_tol, bandwidth, k, j, d): - for m in range(f_max): - # skip no data values - if not (isnan(f[m,k]) or isnan(f[m,j])): - counts[d, i] += 1 - variogram[d, i] += estimator_func(f[m,k] - f[m,j]) + if dist < bin_edges[i] or dist >= bin_edges[i+1]: + continue # skip if not in current bin + for d in range(d_max): + if not dir_test(dim, pos, dist, direction, angles_tol, bandwidth, k, j, d): + continue # skip if not in current direction + for m in range(f_max): + # skip no data values + if not (isnan(f[m,k]) or isnan(f[m,j])): + counts[d, i] += 1 + variogram[d, i] += estimator_func(f[m,k] - f[m,j]) + # once we found a fitting direction + # break the search if directions are separated + if separate_dirs: + break normalization_func_vec(variogram, counts) return np.asarray(variogram), np.asarray(counts) @@ -250,12 +257,13 @@ def unstructured( for j in range(j_max): for k in range(j+1, k_max): dist = distance(dim, pos, j, k) - if dist >= bin_edges[i] and dist < bin_edges[i+1]: - for m in range(f_max): - # skip no data values - if not (isnan(f[m,k]) or isnan(f[m,j])): - counts[i] += 1 - variogram[i] += estimator_func(f[m,k] - f[m,j]) + if dist < bin_edges[i] or dist >= bin_edges[i+1]: + continue # skip if not in current bin + for m in range(f_max): + # skip no data values + if not (isnan(f[m,k]) or isnan(f[m,j])): + counts[i] += 1 + variogram[i] += estimator_func(f[m,k] - f[m,j]) normalization_func(variogram, counts) return np.asarray(variogram), np.asarray(counts) From 6ad2229dc19a5bb7ba3ebc67f5ef5f270bc60d95 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 11 Nov 2020 15:17:08 +0100 Subject: [PATCH 44/44] vario: check if direcitons are separated to optimize search --- gstools/variogram/variogram.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/gstools/variogram/variogram.py b/gstools/variogram/variogram.py index 02c2af81..c3b21351 100644 --- a/gstools/variogram/variogram.py +++ b/gstools/variogram/variogram.py @@ -48,6 +48,18 @@ def _set_estimator(estimator): return cython_estimator +def _separate_dirs_test(direction, angles_tol): + """Check if given directions are separated.""" + if direction is None or direction.shape[0] < 2: + return True + separate_dirs = True + for i in range(direction.shape[0] - 1): + for j in range(i + 1, direction.shape[0]): + s_prod = np.minimum(np.abs(np.dot(direction[i], direction[j])), 1) + separate_dirs &= np.arccos(s_prod) >= 2 * angles_tol + return separate_dirs + + def vario_estimate( pos, field, @@ -276,6 +288,7 @@ def vario_estimate( direction, angles_tol, bandwidth, + separate_dirs=_separate_dirs_test(direction, angles_tol), estimator_type=cython_estimator, ) if dir_no == 1: