From ab44056080ca0fbb3cedca84adff5f377c17b7e9 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Sat, 4 Apr 2020 18:45:32 +0200 Subject: [PATCH 01/41] Readme: docs badge pointing to stable release --- README.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.rst b/README.rst index 8673739..9812510 100644 --- a/README.rst +++ b/README.rst @@ -12,7 +12,7 @@ PyKrige .. image:: https://coveralls.io/repos/github/GeoStat-Framework/PyKrige/badge.svg?branch=master :target: https://coveralls.io/github/GeoStat-Framework/PyKrige?branch=master .. image:: https://readthedocs.org/projects/pykrige/badge/?version=stable - :target: http://pykrige.readthedocs.io/en/latest/?badge=stable + :target: http://pykrige.readthedocs.io/en/stable/?badge=stable :alt: Documentation Status .. image:: https://img.shields.io/badge/code%20style-black-000000.svg :target: https://github.com/psf/black From b6cce9555d2d64fb0ed955063b56cf636da07543 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 5 May 2020 14:59:47 +0200 Subject: [PATCH 02/41] Allways use the pseudo-inverse to solve krige-system by least squares --- pykrige/lib/cok.pyx | 2 +- pykrige/ok.py | 4 ++-- pykrige/ok3d.py | 4 ++-- pykrige/uk.py | 4 ++-- pykrige/uk3d.py | 4 ++-- 5 files changed, 9 insertions(+), 9 deletions(-) diff --git a/pykrige/lib/cok.pyx b/pykrige/lib/cok.pyx index be15dbe..7310d0f 100644 --- a/pykrige/lib/cok.pyx +++ b/pykrige/lib/cok.pyx @@ -35,7 +35,7 @@ cpdef _c_exec_loop(double [:, ::1] a_all, cdef double [::1] variogram_model_parameters = np.asarray(pars['variogram_model_parameters']) - cdef double [::1,:] a_inv = scipy.linalg.inv(a_all) + cdef double [::1,:] a_inv = scipy.linalg.pinvh(a_all) for i in range(npt): # same thing as range(npt) if mask is not defined, otherwise take the non masked elements if mask[i]: diff --git a/pykrige/ok.py b/pykrige/ok.py index 6935ac3..84b160a 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -601,7 +601,7 @@ def _exec_vector(self, a, bd, mask): zero_index = None zero_value = False - a_inv = scipy.linalg.inv(a) + a_inv = scipy.linalg.pinvh(a) if np.any(np.absolute(bd) <= self.eps): zero_value = True @@ -632,7 +632,7 @@ def _exec_loop(self, a, bd_all, mask): zvalues = np.zeros(npt) sigmasq = np.zeros(npt) - a_inv = scipy.linalg.inv(a) + a_inv = scipy.linalg.pinvh(a) for j in np.nonzero(~mask)[ 0 diff --git a/pykrige/ok3d.py b/pykrige/ok3d.py index 679d2c9..dd7eba0 100644 --- a/pykrige/ok3d.py +++ b/pykrige/ok3d.py @@ -586,7 +586,7 @@ def _exec_vector(self, a, bd, mask): zero_index = None zero_value = False - a_inv = scipy.linalg.inv(a) + a_inv = scipy.linalg.pinvh(a) if np.any(np.absolute(bd) <= self.eps): zero_value = True @@ -617,7 +617,7 @@ def _exec_loop(self, a, bd_all, mask): kvalues = np.zeros(npt) sigmasq = np.zeros(npt) - a_inv = scipy.linalg.inv(a) + a_inv = scipy.linalg.pinvh(a) for j in np.nonzero(~mask)[ 0 diff --git a/pykrige/uk.py b/pykrige/uk.py index 1e42e00..ab819c4 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -877,7 +877,7 @@ def _exec_vector(self, a, bd, xy, xy_orig, mask, n_withdrifts, spec_drift_grids) zero_index = None zero_value = False - a_inv = scipy.linalg.inv(a) + a_inv = scipy.linalg.pinvh(a) if np.any(np.absolute(bd) <= self.eps): zero_value = True @@ -961,7 +961,7 @@ def _exec_loop(self, a, bd_all, xy, xy_orig, mask, n_withdrifts, spec_drift_grid zvalues = np.zeros(npt) sigmasq = np.zeros(npt) - a_inv = scipy.linalg.inv(a) + a_inv = scipy.linalg.pinvh(a) for j in np.nonzero(~mask)[ 0 diff --git a/pykrige/uk3d.py b/pykrige/uk3d.py index c94664a..1d096dd 100644 --- a/pykrige/uk3d.py +++ b/pykrige/uk3d.py @@ -701,7 +701,7 @@ def _exec_vector(self, a, bd, xyz, mask, n_withdrifts, spec_drift_grids): zero_index = None zero_value = False - a_inv = scipy.linalg.inv(a) + a_inv = scipy.linalg.pinvh(a) if np.any(np.absolute(bd) <= self.eps): zero_value = True @@ -770,7 +770,7 @@ def _exec_loop(self, a, bd_all, xyz, mask, n_withdrifts, spec_drift_grids): kvalues = np.zeros(npt) sigmasq = np.zeros(npt) - a_inv = scipy.linalg.inv(a) + a_inv = scipy.linalg.pinvh(a) for j in np.nonzero(~mask)[ 0 From 5a409d0b99433a45c2412cde29ddf90a32a93d44 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 5 May 2020 15:27:24 +0200 Subject: [PATCH 03/41] cython: inverse needs to be fortran array --- pykrige/lib/cok.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pykrige/lib/cok.pyx b/pykrige/lib/cok.pyx index 7310d0f..c5fb8e3 100644 --- a/pykrige/lib/cok.pyx +++ b/pykrige/lib/cok.pyx @@ -35,7 +35,7 @@ cpdef _c_exec_loop(double [:, ::1] a_all, cdef double [::1] variogram_model_parameters = np.asarray(pars['variogram_model_parameters']) - cdef double [::1,:] a_inv = scipy.linalg.pinvh(a_all) + cdef double [::1,:] a_inv = np.asfortranarray(scipy.linalg.pinvh(a_all)) for i in range(npt): # same thing as range(npt) if mask is not defined, otherwise take the non masked elements if mask[i]: From a94e8b68a97c04fcf88b731e8a1c6673769c716b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 5 May 2020 16:28:18 +0200 Subject: [PATCH 04/41] test: better checking for 0 --- tests/test_core.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/tests/test_core.py b/tests/test_core.py index d17038b..3178771 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -1369,13 +1369,13 @@ def test_force_exact(): z, ss = ok.execute( "grid", np.arange(0.0, 3.1, 0.1), np.arange(2.1, 3.1, 0.1), backend="vectorized" ) - assert np.any(ss <= 1e-15) - assert not np.any(ss[:9, :30] <= 1e-15) + assert np.any(np.isclose(ss, 0)) + assert not np.any(np.isclose(ss[:9, :30], 0)) assert not np.allclose(z[:9, :30], 0.0) z, ss = ok.execute( "grid", np.arange(0.0, 1.9, 0.1), np.arange(2.1, 3.1, 0.1), backend="vectorized" ) - assert not np.any(ss <= 1e-15) + assert not np.any(np.isclose(ss, 0)) z, ss = ok.execute( "masked", np.arange(2.5, 3.5, 0.1), @@ -1385,7 +1385,7 @@ def test_force_exact(): np.meshgrid(np.arange(2.5, 3.5, 0.1), np.arange(2.5, 3.5, 0.25))[0] == 0.0 ), ) - assert ss[2, 5] <= 1e-15 + assert np.isclose(ss[2, 5], 0) assert not np.allclose(ss, 0.0) z, ss = ok.execute("grid", [1.0, 2.0, 3.0], [1.0, 2.0, 3.0], backend="loop") @@ -1425,13 +1425,13 @@ def test_force_exact(): z, ss = ok.execute( "grid", np.arange(0.0, 3.1, 0.1), np.arange(2.1, 3.1, 0.1), backend="loop" ) - assert np.any(ss <= 1e-15) - assert not np.any(ss[:9, :30] <= 1e-15) + assert np.any(np.isclose(ss, 0)) + assert not np.any(np.isclose(ss[:9, :30], 0)) assert not np.allclose(z[:9, :30], 0.0) z, ss = ok.execute( "grid", np.arange(0.0, 1.9, 0.1), np.arange(2.1, 3.1, 0.1), backend="loop" ) - assert not np.any(ss <= 1e-15) + assert not np.any(np.isclose(ss, 0)) z, ss = ok.execute( "masked", np.arange(2.5, 3.5, 0.1), @@ -1441,7 +1441,7 @@ def test_force_exact(): np.meshgrid(np.arange(2.5, 3.5, 0.1), np.arange(2.5, 3.5, 0.25))[0] == 0.0 ), ) - assert ss[2, 5] <= 1e-15 + assert np.isclose(ss[2, 5], 0) assert not np.allclose(ss, 0.0) uk = UniversalKriging(data[:, 0], data[:, 1], data[:, 2]) @@ -1482,13 +1482,13 @@ def test_force_exact(): z, ss = uk.execute( "grid", np.arange(0.0, 3.1, 0.1), np.arange(2.1, 3.1, 0.1), backend="vectorized" ) - assert np.any(ss <= 1e-15) - assert not np.any(ss[:9, :30] <= 1e-15) + assert np.any(np.isclose(ss, 0)) + assert not np.any(np.isclose(ss[:9, :30], 0)) assert not np.allclose(z[:9, :30], 0.0) z, ss = uk.execute( "grid", np.arange(0.0, 1.9, 0.1), np.arange(2.1, 3.1, 0.1), backend="vectorized" ) - assert not (np.any(ss <= 1e-15)) + assert not (np.any(np.isclose(ss, 0))) z, ss = uk.execute( "masked", np.arange(2.5, 3.5, 0.1), @@ -1498,7 +1498,7 @@ def test_force_exact(): np.meshgrid(np.arange(2.5, 3.5, 0.1), np.arange(2.5, 3.5, 0.25))[0] == 0.0 ), ) - assert ss[2, 5] <= 1e-15 + assert np.isclose(ss[2, 5], 0) assert not np.allclose(ss, 0.0) z, ss = uk.execute("grid", [1.0, 2.0, 3.0], [1.0, 2.0, 3.0], backend="loop") assert z[0, 0] == approx(2.0) @@ -1537,13 +1537,13 @@ def test_force_exact(): z, ss = uk.execute( "grid", np.arange(0.0, 3.1, 0.1), np.arange(2.1, 3.1, 0.1), backend="loop" ) - assert np.any(ss <= 1e-15) - assert not np.any(ss[:9, :30] <= 1e-15) + assert np.any(np.isclose(ss, 0)) + assert not np.any(np.isclose(ss[:9, :30], 0)) assert not np.allclose(z[:9, :30], 0.0) z, ss = uk.execute( "grid", np.arange(0.0, 1.9, 0.1), np.arange(2.1, 3.1, 0.1), backend="loop" ) - assert not np.any(ss <= 1e-15) + assert not np.any(np.isclose(ss, 0)) z, ss = uk.execute( "masked", np.arange(2.5, 3.5, 0.1), @@ -1553,7 +1553,7 @@ def test_force_exact(): np.meshgrid(np.arange(2.5, 3.5, 0.1), np.arange(2.5, 3.5, 0.25))[0] == 0.0 ), ) - assert ss[2, 5] <= 1e-15 + assert np.isclose(ss[2, 5], 0) assert not np.allclose(ss, 0.0) z, ss = core._krige( From 89379ba437f37531319f91129fbac0fdad245ea3 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Thu, 21 May 2020 17:01:49 -0700 Subject: [PATCH 05/41] Adds execute() exact_values option --- pykrige/ok.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/pykrige/ok.py b/pykrige/ok.py index 6935ac3..e9afc44 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -568,7 +568,7 @@ def print_statistics(self): print("Q2 =", self.Q2) print("cR =", self.cR) - def _get_kriging_matrix(self, n): + def _get_kriging_matrix(self, n, exact_values): """Assembles the kriging matrix.""" if self.coordinates_type == "euclidean": @@ -585,11 +585,16 @@ def _get_kriging_matrix(self, n): ) a = np.zeros((n + 1, n + 1)) a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d) - np.fill_diagonal(a, 0.0) + + if self.variogram_model != 'linear' and not exact_values: + np.fill_diagonal(a, self.variogram_model_parameters[2]) + elif self.variogram_model == 'linear' and not exact_values: + np.fill_diagonal(a, self.variogram_model_parameters[1]) + else: + np.fill_diagonal(a, 0.0) a[n, :] = 1.0 a[:, n] = 1.0 a[n, n] = 0.0 - return a def _exec_vector(self, a, bd, mask): @@ -702,6 +707,7 @@ def execute( mask=None, backend="vectorized", n_closest_points=None, + exact_values=True ): """Calculates a kriged grid and the associated variance. @@ -794,13 +800,15 @@ def execute( # If this is not checked, nondescriptive errors emerge # later in the code. raise ValueError("n_closest_points has to be at least two!") + if not isinstance(exact_values, bool): + raise ValueError("exact_values has to be boolean True or False") xpts = np.atleast_1d(np.squeeze(np.array(xpoints, copy=True))) ypts = np.atleast_1d(np.squeeze(np.array(ypoints, copy=True))) n = self.X_ADJUSTED.shape[0] nx = xpts.size ny = ypts.size - a = self._get_kriging_matrix(n) + a = self._get_kriging_matrix(n, exact_values) if style in ["grid", "masked"]: if style == "masked": From 32de8ae602c4f72ae9a619fe4ec14f410d573263 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Thu, 21 May 2020 17:01:56 -0700 Subject: [PATCH 06/41] Adds execute() exact_values option --- pykrige/uk.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/pykrige/uk.py b/pykrige/uk.py index 1e42e00..f76f074 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -808,7 +808,7 @@ def print_statistics(self): print("Q2 =", self.Q2) print("cR =", self.cR) - def _get_kriging_matrix(self, n, n_withdrifts): + def _get_kriging_matrix(self, n, n_withdrifts, exact_values): """Assembles the kriging matrix.""" xy = np.concatenate( @@ -820,7 +820,13 @@ def _get_kriging_matrix(self, n, n_withdrifts): else: a = np.zeros((n_withdrifts, n_withdrifts)) a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d) - np.fill_diagonal(a, 0.0) + + if self.variogram_model != 'linear' and not exact_values: + np.fill_diagonal(a, self.variogram_model_parameters[2]) + elif self.variogram_model == 'linear' and not exact_values: + np.fill_diagonal(a, self.variogram_model_parameters[1]) + else: + np.fill_diagonal(a, 0.0) i = n if self.regional_linear_drift: @@ -1035,6 +1041,7 @@ def execute( mask=None, backend="vectorized", specified_drift_arrays=None, + exact_values=True ): """Calculates a kriged grid and the associated variance. Includes drift terms. @@ -1124,6 +1131,9 @@ def execute( if style != "grid" and style != "masked" and style != "points": raise ValueError("style argument must be 'grid', 'points', or 'masked'") + if not isinstance(exact_values, bool): + raise ValueError("exact_values has to be boolean True or False") + n = self.X_ADJUSTED.shape[0] n_withdrifts = n xpts = np.atleast_1d(np.squeeze(np.array(xpoints, copy=True))) @@ -1140,7 +1150,7 @@ def execute( n_withdrifts += len(self.specified_drift_data_arrays) if self.functional_drift: n_withdrifts += len(self.functional_drift_terms) - a = self._get_kriging_matrix(n, n_withdrifts) + a = self._get_kriging_matrix(n, n_withdrifts, exact_values) if style in ["grid", "masked"]: if style == "masked": From 65e21e13f731886b1b0c04f861d26ad86a3809f5 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Thu, 21 May 2020 17:02:17 -0700 Subject: [PATCH 07/41] Adds basic tests for exact_values option --- tests/test_core.py | 64 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/tests/test_core.py b/tests/test_core.py index d17038b..4fe3a1e 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -538,6 +538,21 @@ def test_ok_execute(sample_data_2d): assert np.amax(ss) != np.amin(ss) assert not np.ma.is_masked(z) + with pytest.raises(ValueError): + ok.execute("grid", gridx, gridy, exact_values="blurg") + + z1, ss1 = ok.execute("grid", gridx, gridy, backend="loop", exact_values=True) + assert_allclose(z1, z) + assert_allclose(ss1, ss) + + z, ss = ok.execute("grid", gridx, gridy, backend="loop", exact_values=False) + shape = (gridy.size, gridx.size) + assert z.shape == shape + assert ss.shape == shape + assert np.amax(z) != np.amin(z) + assert np.amax(ss) != np.amin(ss) + assert not np.ma.is_masked(z) + with pytest.raises(IOError): ok.execute("masked", gridx, gridy, backend="vectorized") mask = np.array([True, False]) @@ -570,6 +585,13 @@ def test_ok_execute(sample_data_2d): assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked + z, ss = ok.execute("masked", gridx, gridy, mask=mask_ref.T, backend="loop", exact_values=False) + assert np.ma.is_masked(z) + assert np.ma.is_masked(ss) + assert z[0, 0] is np.ma.masked + assert ss[0, 0] is np.ma.masked + + with pytest.raises(ValueError): ok.execute( "points", @@ -599,6 +621,12 @@ def test_cython_ok(sample_data_2d): assert_allclose(z1, z2) assert_allclose(ss1, ss2) + z1, ss1 = ok.execute("grid", gridx, gridy, backend="loop", exact_values=False) + z2, ss2 = ok.execute("grid", gridx, gridy, backend="C", exact_values=False) + assert_allclose(z1, z2) + assert_allclose(ss1, ss2) + + closest_points = 4 z1, ss1 = ok.execute( @@ -610,6 +638,14 @@ def test_cython_ok(sample_data_2d): assert_allclose(z1, z2) assert_allclose(ss1, ss2) + z1, ss1 = ok.execute( + "grid", gridx, gridy, backend="loop", n_closest_points=closest_points, exact_values=False + ) + z2, ss2 = ok.execute( + "grid", gridx, gridy, backend="C", n_closest_points=closest_points, exact_values=False + ) + assert_allclose(z1, z2) + assert_allclose(ss1, ss2) def test_uk(validation_ref): @@ -822,6 +858,8 @@ def test_uk_execute(sample_data_2d): uk.execute("blurg", gridx, gridy) with pytest.raises(ValueError): uk.execute("grid", gridx, gridy, backend="mrow") + with pytest.raises(ValueError): + uk.execute("grid", gridx, gridy, exact_values="blurg") z, ss = uk.execute("grid", gridx, gridy, backend="vectorized") shape = (gridy.size, gridx.size) @@ -831,6 +869,18 @@ def test_uk_execute(sample_data_2d): assert np.amax(ss) != np.amin(ss) assert not np.ma.is_masked(z) + z1, ss1 = uk.execute("grid", gridx, gridy, backend="vectorized", exact_values=True) + assert_allclose(z1, z) + assert_allclose(ss1, ss) + + z, ss = uk.execute("grid", gridx, gridy, backend="vectorized", exact_values=False) + shape = (gridy.size, gridx.size) + assert z.shape == shape + assert ss.shape == shape + assert np.amax(z) != np.amin(z) + assert np.amax(ss) != np.amin(ss) + assert not np.ma.is_masked(z) + z, ss = uk.execute("grid", gridx, gridy, backend="loop") shape = (gridy.size, gridx.size) assert z.shape == shape @@ -871,6 +921,12 @@ def test_uk_execute(sample_data_2d): assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked + z, ss = uk.execute("masked", gridx, gridy, mask=mask_ref.T, backend="loop", exact_values=False) + assert np.ma.is_masked(z) + assert np.ma.is_masked(ss) + assert z[0, 0] is np.ma.masked + assert ss[0, 0] is np.ma.masked + with pytest.raises(ValueError): uk.execute( "points", @@ -918,11 +974,19 @@ def test_ok_uk_produce_same_result(validation_ref): assert_allclose(z_ok, z_uk) assert_allclose(ss_ok, ss_uk) + z_uk, ss_uk = uk.execute("grid", gridx, gridy, backend="vectorized", exact_values=False) + assert_allclose(z_ok, z_uk) + assert_allclose(ss_ok, ss_uk) + z_ok, ss_ok = ok.execute("grid", gridx, gridy, backend="loop") z_uk, ss_uk = uk.execute("grid", gridx, gridy, backend="loop") assert_allclose(z_ok, z_uk) assert_allclose(ss_ok, ss_uk) + z_uk, ss_uk = uk.execute("grid", gridx, gridy, backend="loop", exact_values=False) + assert_allclose(z_ok, z_uk) + assert_allclose(ss_ok, ss_uk) + def test_ok_backends_produce_same_result(validation_ref): From c553332529bb6839fbe7ccc3b599d8563bf19926 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Thu, 21 May 2020 17:04:45 -0700 Subject: [PATCH 08/41] black files --- pykrige/ok.py | 6 +++--- pykrige/uk.py | 6 +++--- tests/test_core.py | 29 ++++++++++++++++++++++------- 3 files changed, 28 insertions(+), 13 deletions(-) diff --git a/pykrige/ok.py b/pykrige/ok.py index e9afc44..0ff94ea 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -586,9 +586,9 @@ def _get_kriging_matrix(self, n, exact_values): a = np.zeros((n + 1, n + 1)) a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d) - if self.variogram_model != 'linear' and not exact_values: + if self.variogram_model != "linear" and not exact_values: np.fill_diagonal(a, self.variogram_model_parameters[2]) - elif self.variogram_model == 'linear' and not exact_values: + elif self.variogram_model == "linear" and not exact_values: np.fill_diagonal(a, self.variogram_model_parameters[1]) else: np.fill_diagonal(a, 0.0) @@ -707,7 +707,7 @@ def execute( mask=None, backend="vectorized", n_closest_points=None, - exact_values=True + exact_values=True, ): """Calculates a kriged grid and the associated variance. diff --git a/pykrige/uk.py b/pykrige/uk.py index f76f074..e097755 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -821,9 +821,9 @@ def _get_kriging_matrix(self, n, n_withdrifts, exact_values): a = np.zeros((n_withdrifts, n_withdrifts)) a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d) - if self.variogram_model != 'linear' and not exact_values: + if self.variogram_model != "linear" and not exact_values: np.fill_diagonal(a, self.variogram_model_parameters[2]) - elif self.variogram_model == 'linear' and not exact_values: + elif self.variogram_model == "linear" and not exact_values: np.fill_diagonal(a, self.variogram_model_parameters[1]) else: np.fill_diagonal(a, 0.0) @@ -1041,7 +1041,7 @@ def execute( mask=None, backend="vectorized", specified_drift_arrays=None, - exact_values=True + exact_values=True, ): """Calculates a kriged grid and the associated variance. Includes drift terms. diff --git a/tests/test_core.py b/tests/test_core.py index 4fe3a1e..1f46d45 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -585,13 +585,14 @@ def test_ok_execute(sample_data_2d): assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked - z, ss = ok.execute("masked", gridx, gridy, mask=mask_ref.T, backend="loop", exact_values=False) + z, ss = ok.execute( + "masked", gridx, gridy, mask=mask_ref.T, backend="loop", exact_values=False + ) assert np.ma.is_masked(z) assert np.ma.is_masked(ss) assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked - with pytest.raises(ValueError): ok.execute( "points", @@ -626,7 +627,6 @@ def test_cython_ok(sample_data_2d): assert_allclose(z1, z2) assert_allclose(ss1, ss2) - closest_points = 4 z1, ss1 = ok.execute( @@ -639,14 +639,25 @@ def test_cython_ok(sample_data_2d): assert_allclose(ss1, ss2) z1, ss1 = ok.execute( - "grid", gridx, gridy, backend="loop", n_closest_points=closest_points, exact_values=False + "grid", + gridx, + gridy, + backend="loop", + n_closest_points=closest_points, + exact_values=False, ) z2, ss2 = ok.execute( - "grid", gridx, gridy, backend="C", n_closest_points=closest_points, exact_values=False + "grid", + gridx, + gridy, + backend="C", + n_closest_points=closest_points, + exact_values=False, ) assert_allclose(z1, z2) assert_allclose(ss1, ss2) + def test_uk(validation_ref): # Test to compare UK with linear drift to results from KT3D_H2O. @@ -921,7 +932,9 @@ def test_uk_execute(sample_data_2d): assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked - z, ss = uk.execute("masked", gridx, gridy, mask=mask_ref.T, backend="loop", exact_values=False) + z, ss = uk.execute( + "masked", gridx, gridy, mask=mask_ref.T, backend="loop", exact_values=False + ) assert np.ma.is_masked(z) assert np.ma.is_masked(ss) assert z[0, 0] is np.ma.masked @@ -974,7 +987,9 @@ def test_ok_uk_produce_same_result(validation_ref): assert_allclose(z_ok, z_uk) assert_allclose(ss_ok, ss_uk) - z_uk, ss_uk = uk.execute("grid", gridx, gridy, backend="vectorized", exact_values=False) + z_uk, ss_uk = uk.execute( + "grid", gridx, gridy, backend="vectorized", exact_values=False + ) assert_allclose(z_ok, z_uk) assert_allclose(ss_ok, ss_uk) From 832f787c26746feeee17e51d4005624da61cb991 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Thu, 21 May 2020 17:25:37 -0700 Subject: [PATCH 09/41] Updates docstring --- pykrige/ok.py | 30 ++++++++++++------------------ pykrige/uk.py | 27 +++++++++++---------------- 2 files changed, 23 insertions(+), 34 deletions(-) diff --git a/pykrige/ok.py b/pykrige/ok.py index 0ff94ea..3b30e88 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -712,24 +712,7 @@ def execute( """Calculates a kriged grid and the associated variance. This is now the method that performs the main kriging calculation. - Note that currently measurements (i.e., z values) are considered - 'exact'. This means that, when a specified coordinate for interpolation - is exactly the same as one of the data points, the variogram evaluated - at the point is forced to be zero. Also, the diagonal of the kriging - matrix is also always forced to be zero. In forcing the variogram - evaluated at data points to be zero, we are effectively saying that - there is no variance at that point (no uncertainty, - so the value is 'exact'). - - In the future, the code may include an extra 'exact_values' boolean - flag that can be adjusted to specify whether to treat the measurements - as 'exact'. Setting the flag to false would indicate that the - variogram should not be forced to be zero at zero distance - (i.e., when evaluated at data points). Instead, the uncertainty in - the point will be equal to the nugget. This would mean that the - diagonal of the kriging matrix would be set to - the nugget instead of to zero. - + Parameters ---------- style : str @@ -777,6 +760,17 @@ def execute( with caution. As Kitanidis notes, kriging with a moving window can produce unexpected oddities if the variogram model is not carefully chosen. + exact_values : bool, optional + When exact_values is true and a specified coordinate for interpolation + is exactly the same as one of the data points, the variogram evaluated + at the point is forced to be zero. This effectively says that + there is no variance at that point (no uncertainty, + so the value is 'exact'). Setting exact_values to false indicates that the + variogram should not be forced to be zero at zero distance + (i.e., when evaluated at data points). Instead, the uncertainty in + the point will be equal to the nugget. This means that the + diagonal of the kriging matrix would be set to + the nugget instead of to zero. Returns ------- diff --git a/pykrige/uk.py b/pykrige/uk.py index e097755..ba8315e 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -1047,22 +1047,6 @@ def execute( Includes drift terms. This is now the method that performs the main kriging calculation. - Note that currently measurements (i.e., z values) are considered - 'exact'. This means that, when a specified coordinate for interpolation - is exactly the same as one of the data points, the variogram evaluated - at the point is forced to be zero. Also, the diagonal of the kriging - matrix is also always forced to be zero. In forcing the variogram - evaluated at data points to be zero, we are effectively saying that - there is no variance at that point (no uncertainty, - so the value is 'exact'). - - In the future, the code may include an extra 'exact_values' boolean - flag that can be adjusted to specify whether to treat the measurements - as 'exact'. Setting the flag to false would indicate that the variogram - should not be forced to be zero at zero distance (i.e., when evaluated - at data points). Instead, the uncertainty in the point will be equal to - the nugget. This would mean that the diagonal of the kriging matrix - would be set to the nugget instead of to zero. Parameters ---------- @@ -1112,6 +1096,17 @@ def execute( shape (M, N), where M is the number of y grid-points and N is the number of x grid-points, or shape (M, ) or (N, 1), where M is the number of points at which to evaluate the kriging system. + exact_values : bool, optional + When exact_values is true and a specified coordinate for interpolation + is exactly the same as one of the data points, the variogram evaluated + at the point is forced to be zero. This effectively says that + there is no variance at that point (no uncertainty, + so the value is 'exact'). Setting exact_values to false indicates that the + variogram should not be forced to be zero at zero distance + (i.e., when evaluated at data points). Instead, the uncertainty in + the point will be equal to the nugget. This means that the + diagonal of the kriging matrix would be set to + the nugget instead of to zero. Returns ------- From a067eb11eed24b4115a72610641d3859a20ed14f Mon Sep 17 00:00:00 2001 From: Roman Yurchak Date: Tue, 26 May 2020 13:39:47 +0200 Subject: [PATCH 10/41] Fix CI with scikit-learn related fixes (#157) --- .travis.yml | 1 + examples/krige_cv.py | 4 ++-- pykrige/compat.py | 7 ------- tests/test_api.py | 6 +++++- 4 files changed, 8 insertions(+), 10 deletions(-) diff --git a/.travis.yml b/.travis.yml index 64f379f..d1e4a87 100644 --- a/.travis.yml +++ b/.travis.yml @@ -11,6 +11,7 @@ env: - CIBW_BEFORE_BUILD="pip install numpy==1.17.3 scipy==1.3.2 cython==0.29.14 setuptools" - CIBW_TEST_REQUIRES="pytest scikit-learn gstools" - CIBW_TEST_COMMAND="pytest -v {project}/tests" + - OMP_NUM_THREADS=2 before_install: - | diff --git a/examples/krige_cv.py b/examples/krige_cv.py index b39f02d..0b967c5 100644 --- a/examples/krige_cv.py +++ b/examples/krige_cv.py @@ -8,7 +8,7 @@ import numpy as np from pykrige.rk import Krige -from pykrige.compat import GridSearchCV +from sklearn.model_selection import GridSearchCV # 2D Kring param opt @@ -20,7 +20,7 @@ # "weight": [True, False] } -estimator = GridSearchCV(Krige(), param_dict, verbose=True) +estimator = GridSearchCV(Krige(), param_dict, verbose=True, return_train_score=True) # dummy data X = np.random.randint(0, 400, size=(100, 2)).astype(float) diff --git a/pykrige/compat.py b/pykrige/compat.py index 4fbe3a8..73434aa 100644 --- a/pykrige/compat.py +++ b/pykrige/compat.py @@ -1,7 +1,6 @@ # coding: utf-8 # pylint: disable= invalid-name, unused-import """For compatibility""" -import inspect from functools import partial @@ -11,12 +10,6 @@ from sklearn.model_selection import train_test_split SKLEARN_INSTALLED = True - arg_spec = inspect.getfullargspec(GridSearchCV)[0] - # https://stackoverflow.com/a/56618067/6696397 - if "return_train_score" in arg_spec: - GridSearchCV = partial(GridSearchCV, return_train_score=True) - if "iid" in arg_spec: - GridSearchCV = partial(GridSearchCV, iid=False) except ImportError: SKLEARN_INSTALLED = False diff --git a/tests/test_api.py b/tests/test_api.py index 7905cf8..98d2de0 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -1,10 +1,10 @@ from itertools import product import numpy as np +import pytest from pykrige.rk import Krige from pykrige.rk import threed_krige -from pykrige.compat import GridSearchCV def _method_and_vergiogram(): @@ -15,6 +15,9 @@ def _method_and_vergiogram(): def test_krige(): # dummy data + pytest.importorskip("sklearn") + from sklearn.model_selection import GridSearchCV + np.random.seed(1) X = np.random.randint(0, 400, size=(20, 3)).astype(float) y = 5 * np.random.rand(20) @@ -28,6 +31,7 @@ def test_krige(): n_jobs=-1, pre_dispatch="2*n_jobs", verbose=False, + return_train_score=True, cv=5, ) # run the gridsearch From 5b5420abceeb90054b7a5b1b9b3f81b7feea8c74 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Mon, 22 Jun 2020 22:23:48 -0700 Subject: [PATCH 11/41] Adds exact_values to zero_value conditions --- pykrige/ok.py | 57 +++++++++++++++++++++++++++------------------------ 1 file changed, 30 insertions(+), 27 deletions(-) diff --git a/pykrige/ok.py b/pykrige/ok.py index 3b30e88..bda88ce 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -165,6 +165,7 @@ def __init__( variogram_model="linear", variogram_parameters=None, variogram_function=None, + exact_values=True, nlags=6, weight=False, anisotropy_scaling=1.0, @@ -178,6 +179,9 @@ def __init__( # set up variogram model and parameters... self.variogram_model = variogram_model self.model = None + if not isinstance(exact_values, bool): + raise ValueError("exact_values has to be boolean True or False") + self.exact_values = exact_values # check if a GSTools covariance model is given if hasattr(self.variogram_model, "pykrige_kwargs"): # save the model in the class @@ -568,7 +572,7 @@ def print_statistics(self): print("Q2 =", self.Q2) print("cR =", self.cR) - def _get_kriging_matrix(self, n, exact_values): + def _get_kriging_matrix(self, n): """Assembles the kriging matrix.""" if self.coordinates_type == "euclidean": @@ -586,12 +590,7 @@ def _get_kriging_matrix(self, n, exact_values): a = np.zeros((n + 1, n + 1)) a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d) - if self.variogram_model != "linear" and not exact_values: - np.fill_diagonal(a, self.variogram_model_parameters[2]) - elif self.variogram_model == "linear" and not exact_values: - np.fill_diagonal(a, self.variogram_model_parameters[1]) - else: - np.fill_diagonal(a, 0.0) + np.fill_diagonal(a, 0.0) a[n, :] = 1.0 a[:, n] = 1.0 a[n, n] = 0.0 @@ -614,7 +613,7 @@ def _exec_vector(self, a, bd, mask): b = np.zeros((npt, n + 1, 1)) b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd) - if zero_value: + if zero_value and self.exact_values: b[zero_index[0], zero_index[1], 0] = 0.0 b[:, n, 0] = 1.0 @@ -652,7 +651,7 @@ def _exec_loop(self, a, bd_all, mask): b = np.zeros((n + 1, 1)) b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd) - if zero_value: + if zero_value and self.exact_values: b[zero_index[0], 0] = 0.0 b[n, 0] = 1.0 x = np.dot(a_inv, b) @@ -688,7 +687,7 @@ def _exec_loop_moving_window(self, a_all, bd_all, mask, bd_idx): zero_value = False b = np.zeros((n + 1, 1)) b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd) - if zero_value: + if zero_value and self.exact_values: b[zero_index[0], 0] = 0.0 b[n, 0] = 1.0 @@ -706,13 +705,29 @@ def execute( ypoints, mask=None, backend="vectorized", - n_closest_points=None, - exact_values=True, + n_closest_points=None ): """Calculates a kriged grid and the associated variance. This is now the method that performs the main kriging calculation. - + Note that currently measurements (i.e., z values) are considered + 'exact'. This means that, when a specified coordinate for interpolation + is exactly the same as one of the data points, the variogram evaluated + at the point is forced to be zero. Also, the diagonal of the kriging + matrix is also always forced to be zero. In forcing the variogram + evaluated at data points to be zero, we are effectively saying that + there is no variance at that point (no uncertainty, + so the value is 'exact'). + + In the future, the code may include an extra 'exact_values' boolean + flag that can be adjusted to specify whether to treat the measurements + as 'exact'. Setting the flag to false would indicate that the + variogram should not be forced to be zero at zero distance + (i.e., when evaluated at data points). Instead, the uncertainty in + the point will be equal to the nugget. This would mean that the + diagonal of the kriging matrix would be set to + the nugget instead of to zero. + Parameters ---------- style : str @@ -760,17 +775,6 @@ def execute( with caution. As Kitanidis notes, kriging with a moving window can produce unexpected oddities if the variogram model is not carefully chosen. - exact_values : bool, optional - When exact_values is true and a specified coordinate for interpolation - is exactly the same as one of the data points, the variogram evaluated - at the point is forced to be zero. This effectively says that - there is no variance at that point (no uncertainty, - so the value is 'exact'). Setting exact_values to false indicates that the - variogram should not be forced to be zero at zero distance - (i.e., when evaluated at data points). Instead, the uncertainty in - the point will be equal to the nugget. This means that the - diagonal of the kriging matrix would be set to - the nugget instead of to zero. Returns ------- @@ -794,15 +798,13 @@ def execute( # If this is not checked, nondescriptive errors emerge # later in the code. raise ValueError("n_closest_points has to be at least two!") - if not isinstance(exact_values, bool): - raise ValueError("exact_values has to be boolean True or False") xpts = np.atleast_1d(np.squeeze(np.array(xpoints, copy=True))) ypts = np.atleast_1d(np.squeeze(np.array(ypoints, copy=True))) n = self.X_ADJUSTED.shape[0] nx = xpts.size ny = ypts.size - a = self._get_kriging_matrix(n, exact_values) + a = self._get_kriging_matrix(n) if style in ["grid", "masked"]: if style == "masked": @@ -878,6 +880,7 @@ def execute( "eps", "variogram_model_parameters", "variogram_function", + "exact_values" ] } From ffbd472996a9c1da10e5d78423f848555338efaf Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Mon, 22 Jun 2020 22:24:03 -0700 Subject: [PATCH 12/41] Skip check_vect_b call if not exact_values --- pykrige/lib/cok.pyx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pykrige/lib/cok.pyx b/pykrige/lib/cok.pyx index be15dbe..056c544 100644 --- a/pykrige/lib/cok.pyx +++ b/pykrige/lib/cok.pyx @@ -48,7 +48,8 @@ cpdef _c_exec_loop(double [:, ::1] a_all, b[k] *= -1 b[n] = 1.0 - check_b_vect(n, bd, b, eps) + if not pars['exact_values']: + check_b_vect(n, bd, b, eps) # Do the BLAS matrix-vector multiplication call From 20641a948dd64b3daf22965ee8d3da525ff801a6 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Tue, 23 Jun 2020 11:41:00 -0700 Subject: [PATCH 13/41] Updates documentation, adds ref for non-exact --- pykrige/ok.py | 29 +++++++++-------------------- pykrige/uk.py | 46 ++++++++++++++++++---------------------------- 2 files changed, 27 insertions(+), 48 deletions(-) diff --git a/pykrige/ok.py b/pykrige/ok.py index bda88ce..53c7e63 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -140,11 +140,17 @@ class OrdinaryKriging: coordinates, x is interpreted as longitude and y as latitude coordinates, both given in degree. Longitudes are expected in [0, 360] and latitudes in [-90, 90]. Default is 'euclidean'. + exact_values : bool, optional + If True, interpolation provides input values at input locations. + If False interpolation accounts for variance within input values + at input locations and does not behave as an exact-interpolator [2]. References ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. + .. [2] N. Cressie, Statistics for spatial data, + (Wiley Series in Probability and Statistics, 1993) 137 p. """ eps = 1.0e-10 # Cutoff for comparison to zero @@ -165,7 +171,6 @@ def __init__( variogram_model="linear", variogram_parameters=None, variogram_function=None, - exact_values=True, nlags=6, weight=False, anisotropy_scaling=1.0, @@ -174,14 +179,17 @@ def __init__( enable_plotting=False, enable_statistics=False, coordinates_type="euclidean", + exact_values=True ): # set up variogram model and parameters... self.variogram_model = variogram_model self.model = None + if not isinstance(exact_values, bool): raise ValueError("exact_values has to be boolean True or False") self.exact_values = exact_values + # check if a GSTools covariance model is given if hasattr(self.variogram_model, "pykrige_kwargs"): # save the model in the class @@ -709,25 +717,6 @@ def execute( ): """Calculates a kriged grid and the associated variance. - This is now the method that performs the main kriging calculation. - Note that currently measurements (i.e., z values) are considered - 'exact'. This means that, when a specified coordinate for interpolation - is exactly the same as one of the data points, the variogram evaluated - at the point is forced to be zero. Also, the diagonal of the kriging - matrix is also always forced to be zero. In forcing the variogram - evaluated at data points to be zero, we are effectively saying that - there is no variance at that point (no uncertainty, - so the value is 'exact'). - - In the future, the code may include an extra 'exact_values' boolean - flag that can be adjusted to specify whether to treat the measurements - as 'exact'. Setting the flag to false would indicate that the - variogram should not be forced to be zero at zero distance - (i.e., when evaluated at data points). Instead, the uncertainty in - the point will be equal to the nugget. This would mean that the - diagonal of the kriging matrix would be set to - the nugget instead of to zero. - Parameters ---------- style : str diff --git a/pykrige/uk.py b/pykrige/uk.py index ba8315e..a908882 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -172,11 +172,17 @@ class UniversalKriging: Default is False (off). enable_plotting : boolean, optional Enables plotting to display variogram. Default is False (off). + exact_values : bool, optional + If True, interpolation provides input values at input locations. + If False interpolation accounts for variance within input values + at input locations and does not behave as an exact-interpolator [2]. References ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. + .. [2] N. Cressie, Statistics for spatial data, + (Wiley Series in Probability and Statistics, 1993) 137 p. """ UNBIAS = True # This can be changed to remove the unbiasedness condition @@ -212,6 +218,7 @@ def __init__( functional_drift=None, verbose=False, enable_plotting=False, + exact_values=True ): # Deal with mutable default argument @@ -225,6 +232,11 @@ def __init__( # set up variogram model and parameters... self.variogram_model = variogram_model self.model = None + + if not isinstance(exact_values, bool): + raise ValueError("exact_values has to be boolean True or False") + self.exact_values = exact_values + # check if a GSTools covariance model is given if hasattr(self.variogram_model, "pykrige_kwargs"): # save the model in the class @@ -808,7 +820,7 @@ def print_statistics(self): print("Q2 =", self.Q2) print("cR =", self.cR) - def _get_kriging_matrix(self, n, n_withdrifts, exact_values): + def _get_kriging_matrix(self, n, n_withdrifts): """Assembles the kriging matrix.""" xy = np.concatenate( @@ -821,12 +833,7 @@ def _get_kriging_matrix(self, n, n_withdrifts, exact_values): a = np.zeros((n_withdrifts, n_withdrifts)) a[:n, :n] = -self.variogram_function(self.variogram_model_parameters, d) - if self.variogram_model != "linear" and not exact_values: - np.fill_diagonal(a, self.variogram_model_parameters[2]) - elif self.variogram_model == "linear" and not exact_values: - np.fill_diagonal(a, self.variogram_model_parameters[1]) - else: - np.fill_diagonal(a, 0.0) + np.fill_diagonal(a, 0.0) i = n if self.regional_linear_drift: @@ -894,7 +901,7 @@ def _exec_vector(self, a, bd, xy, xy_orig, mask, n_withdrifts, spec_drift_grids) else: b = np.zeros((npt, n_withdrifts, 1)) b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd) - if zero_value: + if zero_value and self.exact_values: b[zero_index[0], zero_index[1], 0] = 0.0 i = n @@ -985,7 +992,7 @@ def _exec_loop(self, a, bd_all, xy, xy_orig, mask, n_withdrifts, spec_drift_grid else: b = np.zeros((n_withdrifts, 1)) b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd) - if zero_value: + if zero_value and self.exact_values: b[zero_index[0], 0] = 0.0 i = n @@ -1040,14 +1047,11 @@ def execute( ypoints, mask=None, backend="vectorized", - specified_drift_arrays=None, - exact_values=True, + specified_drift_arrays=None ): """Calculates a kriged grid and the associated variance. Includes drift terms. - This is now the method that performs the main kriging calculation. - Parameters ---------- style : str @@ -1096,17 +1100,6 @@ def execute( shape (M, N), where M is the number of y grid-points and N is the number of x grid-points, or shape (M, ) or (N, 1), where M is the number of points at which to evaluate the kriging system. - exact_values : bool, optional - When exact_values is true and a specified coordinate for interpolation - is exactly the same as one of the data points, the variogram evaluated - at the point is forced to be zero. This effectively says that - there is no variance at that point (no uncertainty, - so the value is 'exact'). Setting exact_values to false indicates that the - variogram should not be forced to be zero at zero distance - (i.e., when evaluated at data points). Instead, the uncertainty in - the point will be equal to the nugget. This means that the - diagonal of the kriging matrix would be set to - the nugget instead of to zero. Returns ------- @@ -1126,9 +1119,6 @@ def execute( if style != "grid" and style != "masked" and style != "points": raise ValueError("style argument must be 'grid', 'points', or 'masked'") - if not isinstance(exact_values, bool): - raise ValueError("exact_values has to be boolean True or False") - n = self.X_ADJUSTED.shape[0] n_withdrifts = n xpts = np.atleast_1d(np.squeeze(np.array(xpoints, copy=True))) @@ -1145,7 +1135,7 @@ def execute( n_withdrifts += len(self.specified_drift_data_arrays) if self.functional_drift: n_withdrifts += len(self.functional_drift_terms) - a = self._get_kriging_matrix(n, n_withdrifts, exact_values) + a = self._get_kriging_matrix(n, n_withdrifts) if style in ["grid", "masked"]: if style == "masked": From 4905decc25b28ba851ca12441727730161f8e786 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Tue, 23 Jun 2020 11:41:33 -0700 Subject: [PATCH 14/41] Updates preliminary testing --- tests/test_core.py | 87 +++++++++++++++++++++++++++------------------- 1 file changed, 51 insertions(+), 36 deletions(-) diff --git a/tests/test_core.py b/tests/test_core.py index 1f46d45..42ee968 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -519,6 +519,11 @@ def test_ok_execute(sample_data_2d): ok = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2]) + with pytest.raises(ValueError): + OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], exact_values="blurg") + + ok_non_exact = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], exact_values=False) + with pytest.raises(ValueError): ok.execute("blurg", gridx, gridy) @@ -538,14 +543,11 @@ def test_ok_execute(sample_data_2d): assert np.amax(ss) != np.amin(ss) assert not np.ma.is_masked(z) - with pytest.raises(ValueError): - ok.execute("grid", gridx, gridy, exact_values="blurg") - - z1, ss1 = ok.execute("grid", gridx, gridy, backend="loop", exact_values=True) + z1, ss1 = ok_non_exact.execute("grid", gridx, gridy, backend="loop") assert_allclose(z1, z) assert_allclose(ss1, ss) - z, ss = ok.execute("grid", gridx, gridy, backend="loop", exact_values=False) + z, ss = ok_non_exact.execute("grid", gridx, gridy, backend="loop") shape = (gridy.size, gridx.size) assert z.shape == shape assert ss.shape == shape @@ -585,14 +587,13 @@ def test_ok_execute(sample_data_2d): assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked - z, ss = ok.execute( - "masked", gridx, gridy, mask=mask_ref.T, backend="loop", exact_values=False - ) + z, ss = ok_non_exact.execute("masked", gridx, gridy, mask=mask_ref.T, backend="loop") assert np.ma.is_masked(z) assert np.ma.is_masked(ss) assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked + with pytest.raises(ValueError): ok.execute( "points", @@ -617,16 +618,19 @@ def test_cython_ok(sample_data_2d): data, (gridx, gridy, _), mask_ref = sample_data_2d ok = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2]) + ok_non_exact = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], exact_values=False) + z1, ss1 = ok.execute("grid", gridx, gridy, backend="loop") z2, ss2 = ok.execute("grid", gridx, gridy, backend="C") assert_allclose(z1, z2) assert_allclose(ss1, ss2) - z1, ss1 = ok.execute("grid", gridx, gridy, backend="loop", exact_values=False) - z2, ss2 = ok.execute("grid", gridx, gridy, backend="C", exact_values=False) + z1, ss1 = ok_non_exact.execute("grid", gridx, gridy, backend="loop") + z2, ss2 = ok_non_exact.execute("grid", gridx, gridy, backend="C") assert_allclose(z1, z2) assert_allclose(ss1, ss2) + closest_points = 4 z1, ss1 = ok.execute( @@ -638,26 +642,15 @@ def test_cython_ok(sample_data_2d): assert_allclose(z1, z2) assert_allclose(ss1, ss2) - z1, ss1 = ok.execute( - "grid", - gridx, - gridy, - backend="loop", - n_closest_points=closest_points, - exact_values=False, + z1, ss1 = ok_non_exact.execute( + "grid", gridx, gridy, backend="loop", n_closest_points=closest_points ) - z2, ss2 = ok.execute( - "grid", - gridx, - gridy, - backend="C", - n_closest_points=closest_points, - exact_values=False, + z2, ss2 = ok_non_exact.execute( + "grid", gridx, gridy, backend="C", n_closest_points=closest_points ) assert_allclose(z1, z2) assert_allclose(ss1, ss2) - def test_uk(validation_ref): # Test to compare UK with linear drift to results from KT3D_H2O. @@ -865,12 +858,29 @@ def test_uk_execute(sample_data_2d): drift_terms=["regional_linear"], ) + with pytest.raises(ValueError): + UniversalKriging( + data[:, 0], + data[:, 1], + data[:, 2], + variogram_model="linear", + drift_terms=["regional_linear"], + exact_values="blurg" + ) + + uk_non_exact = UniversalKriging( + data[:, 0], + data[:, 1], + data[:, 2], + variogram_model="linear", + drift_terms=["regional_linear"] + ) + + with pytest.raises(ValueError): uk.execute("blurg", gridx, gridy) with pytest.raises(ValueError): uk.execute("grid", gridx, gridy, backend="mrow") - with pytest.raises(ValueError): - uk.execute("grid", gridx, gridy, exact_values="blurg") z, ss = uk.execute("grid", gridx, gridy, backend="vectorized") shape = (gridy.size, gridx.size) @@ -880,11 +890,11 @@ def test_uk_execute(sample_data_2d): assert np.amax(ss) != np.amin(ss) assert not np.ma.is_masked(z) - z1, ss1 = uk.execute("grid", gridx, gridy, backend="vectorized", exact_values=True) + z1, ss1 = uk_non_exact.execute("grid", gridx, gridy, backend="vectorized") assert_allclose(z1, z) assert_allclose(ss1, ss) - z, ss = uk.execute("grid", gridx, gridy, backend="vectorized", exact_values=False) + z, ss = uk_non_exact.execute("grid", gridx, gridy, backend="vectorized") shape = (gridy.size, gridx.size) assert z.shape == shape assert ss.shape == shape @@ -932,9 +942,7 @@ def test_uk_execute(sample_data_2d): assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked - z, ss = uk.execute( - "masked", gridx, gridy, mask=mask_ref.T, backend="loop", exact_values=False - ) + z, ss = uk_non_exact.execute("masked", gridx, gridy, mask=mask_ref.T, backend="loop") assert np.ma.is_masked(z) assert np.ma.is_masked(ss) assert z[0, 0] is np.ma.masked @@ -983,13 +991,20 @@ def test_ok_uk_produce_same_result(validation_ref): verbose=False, enable_plotting=False, ) + uk_non_exact = UniversalKriging( + data[:, 0], + data[:, 1], + data[:, 2], + variogram_model="linear", + verbose=False, + enable_plotting=False, + exact_values=False + ) z_uk, ss_uk = uk.execute("grid", gridx, gridy, backend="vectorized") assert_allclose(z_ok, z_uk) assert_allclose(ss_ok, ss_uk) - z_uk, ss_uk = uk.execute( - "grid", gridx, gridy, backend="vectorized", exact_values=False - ) + z_uk, ss_uk = uk_non_exact.execute("grid", gridx, gridy, backend="vectorized") assert_allclose(z_ok, z_uk) assert_allclose(ss_ok, ss_uk) @@ -998,7 +1013,7 @@ def test_ok_uk_produce_same_result(validation_ref): assert_allclose(z_ok, z_uk) assert_allclose(ss_ok, ss_uk) - z_uk, ss_uk = uk.execute("grid", gridx, gridy, backend="loop", exact_values=False) + z_uk, ss_uk = uk_non_exact.execute("grid", gridx, gridy, backend="loop") assert_allclose(z_ok, z_uk) assert_allclose(ss_ok, ss_uk) From 32e5ac5cbcb2a197e2fa06f28390140e0b302cda Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Tue, 23 Jun 2020 13:24:49 -0700 Subject: [PATCH 15/41] Adds exact_values --- pykrige/ok3d.py | 18 +++++++++++++++--- pykrige/uk3d.py | 16 ++++++++++++++-- 2 files changed, 29 insertions(+), 5 deletions(-) diff --git a/pykrige/ok3d.py b/pykrige/ok3d.py index 679d2c9..d488acb 100644 --- a/pykrige/ok3d.py +++ b/pykrige/ok3d.py @@ -151,11 +151,17 @@ class OrdinaryKriging3D: Default is False (off). enable_plotting : bool, optional Enables plotting to display variogram. Default is False (off). + exact_values : bool, optional + If True, interpolation provides input values at input locations. + If False interpolation accounts for variance within input values + at input locations and does not behave as an exact-interpolator [2]. References ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. + .. [2] N. Cressie, Statistics for spatial data, + (Wiley Series in Probability and Statistics, 1993) 137 p. """ eps = 1.0e-10 # Cutoff for comparison to zero @@ -186,11 +192,17 @@ def __init__( anisotropy_angle_z=0.0, verbose=False, enable_plotting=False, + exact_values=True ): # set up variogram model and parameters... self.variogram_model = variogram_model self.model = None + + if not isinstance(exact_values, bool): + raise ValueError("exact_values has to be boolean True or False") + self.exact_values = exact_values + # check if a GSTools covariance model is given if hasattr(self.variogram_model, "pykrige_kwargs"): # save the model in the class @@ -594,7 +606,7 @@ def _exec_vector(self, a, bd, mask): b = np.zeros((npt, n + 1, 1)) b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd) - if zero_value: + if zero_value and self.exact_values: b[zero_index[0], zero_index[1], 0] = 0.0 b[:, n, 0] = 1.0 @@ -632,7 +644,7 @@ def _exec_loop(self, a, bd_all, mask): b = np.zeros((n + 1, 1)) b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd) - if zero_value: + if zero_value and self.exact_values: b[zero_index[0], 0] = 0.0 b[n, 0] = 1.0 @@ -669,7 +681,7 @@ def _exec_loop_moving_window(self, a_all, bd_all, mask, bd_idx): zero_index = None b = np.zeros((n + 1, 1)) b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd) - if zero_value: + if zero_value and self.exact_values: b[zero_index[0], 0] = 0.0 b[n, 0] = 1.0 diff --git a/pykrige/uk3d.py b/pykrige/uk3d.py index c94664a..e50e47c 100644 --- a/pykrige/uk3d.py +++ b/pykrige/uk3d.py @@ -166,11 +166,17 @@ class UniversalKriging3D: Default is False (off). enable_plotting : boolean, optional Enables plotting to display variogram. Default is False (off). + exact_values : bool, optional + If True, interpolation provides input values at input locations. + If False interpolation accounts for variance within input values + at input locations and does not behave as an exact-interpolator [2]. References ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. + .. [2] N. Cressie, Statistics for spatial data, + (Wiley Series in Probability and Statistics, 1993) 137 p. """ UNBIAS = True # This can be changed to remove the unbiasedness condition @@ -206,6 +212,7 @@ def __init__( functional_drift=None, verbose=False, enable_plotting=False, + exact_values=True ): # Deal with mutable default argument @@ -219,6 +226,11 @@ def __init__( # set up variogram model and parameters... self.variogram_model = variogram_model self.model = None + + if not isinstance(exact_values, bool): + raise ValueError("exact_values has to be boolean True or False") + self.exact_values = exact_values + # check if a GSTools covariance model is given if hasattr(self.variogram_model, "pykrige_kwargs"): # save the model in the class @@ -712,7 +724,7 @@ def _exec_vector(self, a, bd, xyz, mask, n_withdrifts, spec_drift_grids): else: b = np.zeros((npt, n_withdrifts, 1)) b[:, :n, 0] = -self.variogram_function(self.variogram_model_parameters, bd) - if zero_value: + if zero_value and self.exact_values: b[zero_index[0], zero_index[1], 0] = 0.0 i = n @@ -788,7 +800,7 @@ def _exec_loop(self, a, bd_all, xyz, mask, n_withdrifts, spec_drift_grids): else: b = np.zeros((n_withdrifts, 1)) b[:n, 0] = -self.variogram_function(self.variogram_model_parameters, bd) - if zero_value: + if zero_value and self.exact_values: b[zero_index[0], 0] = 0.0 i = n From 71b8dac9b04c51ca574872c513217b03fbc48fc6 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Tue, 23 Jun 2020 13:25:24 -0700 Subject: [PATCH 16/41] Adds test for comparison of backends --- tests/test_core.py | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/tests/test_core.py b/tests/test_core.py index 42ee968..6e92f70 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -1836,6 +1836,28 @@ def test_ok3d(validation_ref): variogram_model="exponential", variogram_parameters=[500.0, 3000.0, 0.0], ) + + with pytest.raises(ValueError): + OrdinaryKriging3D( + data[:, 0], + data[:, 1], + np.zeros(data[:, 1].shape), + data[:, 2], + variogram_model="exponential", + variogram_parameters=[500.0, 3000.0, 0.0], + exact_values="blurg" + ) + + ok3d_non_exact = OrdinaryKriging3D( + data[:, 0], + data[:, 1], + np.zeros(data[:, 1].shape), + data[:, 2], + variogram_model="exponential", + variogram_parameters=[500.0, 3000.0, 0.0], + exact_values=False + ) + k, ss = k3d.execute( "grid", gridx_ref, gridy_ref, np.array([0.0]), backend="vectorized" ) @@ -2066,6 +2088,17 @@ def test_ok3d_backends_produce_same_result(sample_data_3d): k3d = OrdinaryKriging3D( data[:, 0], data[:, 1], data[:, 2], data[:, 3], variogram_model="linear" ) + + ok3d_non_exact = OrdinaryKriging3D( + data[:, 0], + data[:, 1], + np.zeros(data[:, 1].shape), + data[:, 2], + variogram_model="exponential", + variogram_parameters=[500.0, 3000.0, 0.0], + exact_values=False + ) + k_k3d_v, ss_k3d_v = k3d.execute( "grid", gridx_ref, gridy_ref, gridz_ref, backend="vectorized" ) @@ -2075,6 +2108,11 @@ def test_ok3d_backends_produce_same_result(sample_data_3d): assert_allclose(k_k3d_v, k_k3d_l, rtol=1e-05, atol=1e-8) assert_allclose(ss_k3d_v, ss_k3d_l, rtol=1e-05, atol=1e-8) + k, ss = ok3d_non_exact.execute("grid", np.arange(10.0), np.arange(10.0), np.arange(10.0), backend="loop") + k1, ss1 = ok3d_non_exact.execute("grid", np.arange(10.0), np.arange(10.0), np.arange(10.0), backend="vectorized") + assert_allclose(k1, k) + assert_allclose(ss1, ss) + def test_ok3d_execute(sample_data_3d): @@ -2085,6 +2123,7 @@ def test_ok3d_execute(sample_data_3d): with pytest.raises(ValueError): k3d.execute("blurg", gridx_ref, gridy_ref, gridz_ref) + k, ss = k3d.execute("grid", gridx_ref, gridy_ref, gridz_ref, backend="vectorized") shape = (gridz_ref.size, gridy_ref.size, gridx_ref.size) assert k.shape == shape From 8445ba561cf808a748caf46e8359309169d1d2ed Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Tue, 23 Jun 2020 14:13:15 -0700 Subject: [PATCH 17/41] Black --- pykrige/ok.py | 8 ++++---- pykrige/ok3d.py | 2 +- pykrige/uk.py | 4 ++-- pykrige/uk3d.py | 2 +- tests/test_core.py | 49 +++++++++++++++++++++++++++------------------- 5 files changed, 37 insertions(+), 28 deletions(-) diff --git a/pykrige/ok.py b/pykrige/ok.py index 53c7e63..fb0b417 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -179,7 +179,7 @@ def __init__( enable_plotting=False, enable_statistics=False, coordinates_type="euclidean", - exact_values=True + exact_values=True, ): # set up variogram model and parameters... @@ -189,7 +189,7 @@ def __init__( if not isinstance(exact_values, bool): raise ValueError("exact_values has to be boolean True or False") self.exact_values = exact_values - + # check if a GSTools covariance model is given if hasattr(self.variogram_model, "pykrige_kwargs"): # save the model in the class @@ -713,7 +713,7 @@ def execute( ypoints, mask=None, backend="vectorized", - n_closest_points=None + n_closest_points=None, ): """Calculates a kriged grid and the associated variance. @@ -869,7 +869,7 @@ def execute( "eps", "variogram_model_parameters", "variogram_function", - "exact_values" + "exact_values", ] } diff --git a/pykrige/ok3d.py b/pykrige/ok3d.py index d488acb..ea366a2 100644 --- a/pykrige/ok3d.py +++ b/pykrige/ok3d.py @@ -192,7 +192,7 @@ def __init__( anisotropy_angle_z=0.0, verbose=False, enable_plotting=False, - exact_values=True + exact_values=True, ): # set up variogram model and parameters... diff --git a/pykrige/uk.py b/pykrige/uk.py index a908882..44c3760 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -218,7 +218,7 @@ def __init__( functional_drift=None, verbose=False, enable_plotting=False, - exact_values=True + exact_values=True, ): # Deal with mutable default argument @@ -1047,7 +1047,7 @@ def execute( ypoints, mask=None, backend="vectorized", - specified_drift_arrays=None + specified_drift_arrays=None, ): """Calculates a kriged grid and the associated variance. Includes drift terms. diff --git a/pykrige/uk3d.py b/pykrige/uk3d.py index e50e47c..4b69bd6 100644 --- a/pykrige/uk3d.py +++ b/pykrige/uk3d.py @@ -212,7 +212,7 @@ def __init__( functional_drift=None, verbose=False, enable_plotting=False, - exact_values=True + exact_values=True, ): # Deal with mutable default argument diff --git a/tests/test_core.py b/tests/test_core.py index 6e92f70..3edea16 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -522,7 +522,9 @@ def test_ok_execute(sample_data_2d): with pytest.raises(ValueError): OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], exact_values="blurg") - ok_non_exact = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], exact_values=False) + ok_non_exact = OrdinaryKriging( + data[:, 0], data[:, 1], data[:, 2], exact_values=False + ) with pytest.raises(ValueError): ok.execute("blurg", gridx, gridy) @@ -587,13 +589,14 @@ def test_ok_execute(sample_data_2d): assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked - z, ss = ok_non_exact.execute("masked", gridx, gridy, mask=mask_ref.T, backend="loop") + z, ss = ok_non_exact.execute( + "masked", gridx, gridy, mask=mask_ref.T, backend="loop" + ) assert np.ma.is_masked(z) assert np.ma.is_masked(ss) assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked - with pytest.raises(ValueError): ok.execute( "points", @@ -618,7 +621,9 @@ def test_cython_ok(sample_data_2d): data, (gridx, gridy, _), mask_ref = sample_data_2d ok = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2]) - ok_non_exact = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], exact_values=False) + ok_non_exact = OrdinaryKriging( + data[:, 0], data[:, 1], data[:, 2], exact_values=False + ) z1, ss1 = ok.execute("grid", gridx, gridy, backend="loop") z2, ss2 = ok.execute("grid", gridx, gridy, backend="C") @@ -630,7 +635,6 @@ def test_cython_ok(sample_data_2d): assert_allclose(z1, z2) assert_allclose(ss1, ss2) - closest_points = 4 z1, ss1 = ok.execute( @@ -651,6 +655,7 @@ def test_cython_ok(sample_data_2d): assert_allclose(z1, z2) assert_allclose(ss1, ss2) + def test_uk(validation_ref): # Test to compare UK with linear drift to results from KT3D_H2O. @@ -865,7 +870,7 @@ def test_uk_execute(sample_data_2d): data[:, 2], variogram_model="linear", drift_terms=["regional_linear"], - exact_values="blurg" + exact_values="blurg", ) uk_non_exact = UniversalKriging( @@ -873,10 +878,9 @@ def test_uk_execute(sample_data_2d): data[:, 1], data[:, 2], variogram_model="linear", - drift_terms=["regional_linear"] + drift_terms=["regional_linear"], ) - with pytest.raises(ValueError): uk.execute("blurg", gridx, gridy) with pytest.raises(ValueError): @@ -942,7 +946,9 @@ def test_uk_execute(sample_data_2d): assert z[0, 0] is np.ma.masked assert ss[0, 0] is np.ma.masked - z, ss = uk_non_exact.execute("masked", gridx, gridy, mask=mask_ref.T, backend="loop") + z, ss = uk_non_exact.execute( + "masked", gridx, gridy, mask=mask_ref.T, backend="loop" + ) assert np.ma.is_masked(z) assert np.ma.is_masked(ss) assert z[0, 0] is np.ma.masked @@ -998,7 +1004,7 @@ def test_ok_uk_produce_same_result(validation_ref): variogram_model="linear", verbose=False, enable_plotting=False, - exact_values=False + exact_values=False, ) z_uk, ss_uk = uk.execute("grid", gridx, gridy, backend="vectorized") assert_allclose(z_ok, z_uk) @@ -1844,8 +1850,8 @@ def test_ok3d(validation_ref): np.zeros(data[:, 1].shape), data[:, 2], variogram_model="exponential", - variogram_parameters=[500.0, 3000.0, 0.0], - exact_values="blurg" + variogram_parameters=[500.0, 3000.0, 0.0], + exact_values="blurg", ) ok3d_non_exact = OrdinaryKriging3D( @@ -1854,8 +1860,8 @@ def test_ok3d(validation_ref): np.zeros(data[:, 1].shape), data[:, 2], variogram_model="exponential", - variogram_parameters=[500.0, 3000.0, 0.0], - exact_values=False + variogram_parameters=[500.0, 3000.0, 0.0], + exact_values=False, ) k, ss = k3d.execute( @@ -2095,10 +2101,10 @@ def test_ok3d_backends_produce_same_result(sample_data_3d): np.zeros(data[:, 1].shape), data[:, 2], variogram_model="exponential", - variogram_parameters=[500.0, 3000.0, 0.0], - exact_values=False + variogram_parameters=[500.0, 3000.0, 0.0], + exact_values=False, ) - + k_k3d_v, ss_k3d_v = k3d.execute( "grid", gridx_ref, gridy_ref, gridz_ref, backend="vectorized" ) @@ -2108,8 +2114,12 @@ def test_ok3d_backends_produce_same_result(sample_data_3d): assert_allclose(k_k3d_v, k_k3d_l, rtol=1e-05, atol=1e-8) assert_allclose(ss_k3d_v, ss_k3d_l, rtol=1e-05, atol=1e-8) - k, ss = ok3d_non_exact.execute("grid", np.arange(10.0), np.arange(10.0), np.arange(10.0), backend="loop") - k1, ss1 = ok3d_non_exact.execute("grid", np.arange(10.0), np.arange(10.0), np.arange(10.0), backend="vectorized") + k, ss = ok3d_non_exact.execute( + "grid", np.arange(10.0), np.arange(10.0), np.arange(10.0), backend="loop" + ) + k1, ss1 = ok3d_non_exact.execute( + "grid", np.arange(10.0), np.arange(10.0), np.arange(10.0), backend="vectorized" + ) assert_allclose(k1, k) assert_allclose(ss1, ss) @@ -2123,7 +2133,6 @@ def test_ok3d_execute(sample_data_3d): with pytest.raises(ValueError): k3d.execute("blurg", gridx_ref, gridy_ref, gridz_ref) - k, ss = k3d.execute("grid", gridx_ref, gridy_ref, gridz_ref, backend="vectorized") shape = (gridz_ref.size, gridy_ref.size, gridx_ref.size) assert k.shape == shape From 47630266141812cc5e34842bd4c429b43a5ed335 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Wed, 24 Jun 2020 09:45:57 -0700 Subject: [PATCH 18/41] Updates doc with nugget hint --- pykrige/ok.py | 10 ++++++++-- pykrige/ok3d.py | 10 ++++++++-- pykrige/uk.py | 11 ++++++++--- pykrige/uk3d.py | 10 ++++++++-- 4 files changed, 32 insertions(+), 9 deletions(-) diff --git a/pykrige/ok.py b/pykrige/ok.py index fb0b417..b34e4c5 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -15,6 +15,8 @@ ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. +.. [2] N. Cressie, Statistics for spatial data, + (Wiley Series in Probability and Statistics, 1993) 137 p. Copyright (c) 2015-2020, PyKrige Developers """ @@ -142,8 +144,12 @@ class OrdinaryKriging: [0, 360] and latitudes in [-90, 90]. Default is 'euclidean'. exact_values : bool, optional If True, interpolation provides input values at input locations. - If False interpolation accounts for variance within input values - at input locations and does not behave as an exact-interpolator [2]. + If False, interpolation accounts for variance/nugget within input + values at input locations and does not behave as an + exact-interpolator [2]. Note that this only has an effect if + there is variance/nugget present within the input data since it is + interpreted as measurement error. If the nugget is zero, the kriged + field will behave as an exact interpolator. References ---------- diff --git a/pykrige/ok3d.py b/pykrige/ok3d.py index ea366a2..396f0a3 100644 --- a/pykrige/ok3d.py +++ b/pykrige/ok3d.py @@ -14,6 +14,8 @@ ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. +.. [2] N. Cressie, Statistics for spatial data, +(Wiley Series in Probability and Statistics, 1993) 137 p. Copyright (c) 2015-2020, PyKrige Developers """ @@ -153,8 +155,12 @@ class OrdinaryKriging3D: Enables plotting to display variogram. Default is False (off). exact_values : bool, optional If True, interpolation provides input values at input locations. - If False interpolation accounts for variance within input values - at input locations and does not behave as an exact-interpolator [2]. + If False, interpolation accounts for variance/nugget within input + values at input locations and does not behave as an + exact-interpolator [2]. Note that this only has an effect if + there is variance/nugget present within the input data since it is + interpreted as measurement error. If the nugget is zero, the kriged + field will behave as an exact interpolator. References ---------- diff --git a/pykrige/uk.py b/pykrige/uk.py index 44c3760..7858e17 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -15,7 +15,8 @@ ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. - +.. [2] N. Cressie, Statistics for spatial data, + (Wiley Series in Probability and Statistics, 1993) 137 p. Copyright (c) 2015-2020, PyKrige Developers """ import numpy as np @@ -174,8 +175,12 @@ class UniversalKriging: Enables plotting to display variogram. Default is False (off). exact_values : bool, optional If True, interpolation provides input values at input locations. - If False interpolation accounts for variance within input values - at input locations and does not behave as an exact-interpolator [2]. + If False, interpolation accounts for variance/nugget within input + values at input locations and does not behave as an + exact-interpolator [2]. Note that this only has an effect if + there is variance/nugget present within the input data since it is + interpreted as measurement error. If the nugget is zero, the kriged + field will behave as an exact interpolator. References ---------- diff --git a/pykrige/uk3d.py b/pykrige/uk3d.py index 4b69bd6..b78d6ca 100644 --- a/pykrige/uk3d.py +++ b/pykrige/uk3d.py @@ -14,6 +14,8 @@ ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. +.. [2] N. Cressie, Statistics for spatial data, + (Wiley Series in Probability and Statistics, 1993) 137 p. Copyright (c) 2015-2020, PyKrige Developers """ @@ -168,8 +170,12 @@ class UniversalKriging3D: Enables plotting to display variogram. Default is False (off). exact_values : bool, optional If True, interpolation provides input values at input locations. - If False interpolation accounts for variance within input values - at input locations and does not behave as an exact-interpolator [2]. + If False, interpolation accounts for variance/nugget within input + values at input locations and does not behave as an + exact-interpolator [2]. Note that this only has an effect if + there is variance/nugget present within the input data since it is + interpreted as measurement error. If the nugget is zero, the kriged + field will behave as an exact interpolator. References ---------- From f23ef041955730f107e24df5a7f868bc63ba4158 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Wed, 24 Jun 2020 09:48:52 -0700 Subject: [PATCH 19/41] Adds exact_values example --- examples/exact_values_example_1D.py | 50 +++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 examples/exact_values_example_1D.py diff --git a/examples/exact_values_example_1D.py b/examples/exact_values_example_1D.py new file mode 100644 index 0000000..be41d75 --- /dev/null +++ b/examples/exact_values_example_1D.py @@ -0,0 +1,50 @@ +from pykrige.ok import OrdinaryKriging +import matplotlib.pyplot as plt +import numpy as np + +plt.style.use("ggplot") + +np.random.seed(42) + +x = np.linspace(0, 12.5, 50) +xpred = np.linspace(0, 12.5, 393) +y = np.sin(x) * np.exp(-0.25 * x) + np.random.normal(-0.25, 0.25, 50) + +# compare OrdinaryKriging as an exact and non exact interpolator +uk = OrdinaryKriging( + x, np.zeros(x.shape), y, variogram_model="linear", exact_values=False +) +uk_exact = OrdinaryKriging(x, np.zeros(x.shape), y, variogram_model="linear") + +y_pred, y_std = uk.execute("grid", xpred, np.array([0.0]), backend="loop") +y_pred_exact, y_std_exact = uk_exact.execute( + "grid", xpred, np.array([0.0]), backend="loop" +) + + +y_pred = np.squeeze(y_pred) +y_std = np.squeeze(y_std) + +y_pred_exact = np.squeeze(y_pred_exact) +y_std_exact = np.squeeze(y_std_exact) + + +fig, ax = plt.subplots(1, 1, figsize=(10, 4)) + +ax.scatter(x, y, label="Input Data") +ax.plot(xpred, y_pred_exact, label="Exact Prediction") +ax.plot(xpred, y_pred, label="Non Exact Prediction") + +ax.fill_between( + xpred, + y_pred - 3 * y_std, + y_pred + 3 * y_std, + alpha=0.3, + label="Confidence interval", +) +ax.legend(loc=9) +ax.set_ylim(-1.8, 1.3) +ax.legend(loc=9) +plt.xlabel("X") +plt.ylabel("Field") +plt.show() From 4444ec6c558872a6fcc6455d52c9e44d6b186c88 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Fri, 26 Jun 2020 11:28:15 -0700 Subject: [PATCH 20/41] Adds non_exact testing --- tests/test_core.py | 55 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/tests/test_core.py b/tests/test_core.py index 3edea16..5a13b04 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -436,6 +436,61 @@ def test_core_krige_3d(): assert ss == approx(0.0, rel=1e-3) +def test_non_exact(): + # custom data for this test + data = np.array( + [[0.0, 0.0, 0.47], [1.5, 1.5, 0.56], [3, 3, 0.74], [4.5, 4.5, 1.47],] + ) + + # construct grid points so diagonal + # is identical to input points + gridx = np.arange(0.0, 4.51, 1.5) + gridy = np.arange(0.0, 4.51, 1.5) + + ok = OrdinaryKriging( + data[:, 0], + data[:, 1], + data[:, 2], + variogram_model="exponential", + variogram_parameters=[500.0, 3000.0, 5.0], + ) + z, ss = ok.execute("grid", gridx, gridy, backend="vectorized") + + ok_non_exact = OrdinaryKriging( + data[:, 0], + data[:, 1], + data[:, 2], + variogram_model="exponential", + variogram_parameters=[500.0, 3000.0, 5.0], + exact_values=False, + ) + z_non_exact, ss_non_exact = ok_non_exact.execute( + "grid", gridx, gridy, backend="vectorized" + ) + + in_values = np.diag(z) + + # test that krig field + # at input location are identical + # to the inputs themselves with + # exact_values == True + assert_allclose(in_values, data[:, 2]) + + # test that krig field + # at input location are different + # than the inputs themselves + # with exact_values == False + assert ~np.allclose(in_values, data[:, 2]) + + # test that off diagonal values are the same + # by filling with dummy value and comparing + # each entry in array + np.fill_diagonal(z, 0.0) + np.fill_diagonal(z_non_exact, 0.0) + + assert_allclose(z, z_non_exact) + + def test_ok(validation_ref): # Test to compare OK results to those obtained using KT3D_H2O. From f1dda1e12e0f2b4b525a8f474066d8a077c78f48 Mon Sep 17 00:00:00 2001 From: Nic Annau Date: Mon, 29 Jun 2020 07:50:46 -0700 Subject: [PATCH 21/41] Adds headline and description --- examples/exact_values_example_1D.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/examples/exact_values_example_1D.py b/examples/exact_values_example_1D.py index be41d75..5ceeadc 100644 --- a/examples/exact_values_example_1D.py +++ b/examples/exact_values_example_1D.py @@ -1,3 +1,12 @@ +# -*- coding: utf-8 -*- +""" +Exact Values +================= + +PyKrige demonstration and usage +as a non-exact interpolator in 1D. +""" + from pykrige.ok import OrdinaryKriging import matplotlib.pyplot as plt import numpy as np From 74a59835cbc427937418ee3032e1170b001af112 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 15 Jul 2020 19:32:09 +0200 Subject: [PATCH 22/41] travis: don't build PRs of branches twice --- .travis.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.travis.yml b/.travis.yml index d1e4a87..4e41a65 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,3 +1,6 @@ +# do not build pull request from base repo twice (only build PRs of forks) +if: type != pull_request OR fork = true + language: python python: 3.8 From ecff11fcc8164e45330c758d80cd46476aa81f0b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 15 Jul 2020 19:32:50 +0200 Subject: [PATCH 23/41] core: add dict with p-inv routines --- pykrige/core.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pykrige/core.py b/pykrige/core.py index 781c55c..1e02b02 100755 --- a/pykrige/core.py +++ b/pykrige/core.py @@ -25,10 +25,15 @@ import numpy as np from scipy.spatial.distance import pdist, squareform, cdist from scipy.optimize import least_squares +import scipy.linalg as spl + eps = 1.0e-10 # Cutoff for comparison to zero +P_INV = {1: spl.pinv, 2: spl.pinv2, 3: spl.pinvh} + + def great_circle_distance(lon1, lat1, lon2, lat2): """Calculate the great circle distance between one or multiple pairs of points given in spherical coordinates. Spherical coordinates are expected From 94b777a069ed0cdbe7e3a4519a530e38fc74bc4b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 15 Jul 2020 19:33:41 +0200 Subject: [PATCH 24/41] pseudo inv: add switch and selector for usage of pseudo-inverse matrix --- pykrige/lib/cok.pyx | 16 +++++++++++++++- pykrige/ok.py | 45 ++++++++++++++++++++++++++++++++++++++------- pykrige/ok3d.py | 43 ++++++++++++++++++++++++++++++++++++------- pykrige/uk.py | 43 ++++++++++++++++++++++++++++++++++++------- pykrige/uk3d.py | 43 ++++++++++++++++++++++++++++++++++++------- 5 files changed, 161 insertions(+), 29 deletions(-) diff --git a/pykrige/lib/cok.pyx b/pykrige/lib/cok.pyx index 14b3fc7..80e0b21 100644 --- a/pykrige/lib/cok.pyx +++ b/pykrige/lib/cok.pyx @@ -35,7 +35,21 @@ cpdef _c_exec_loop(double [:, ::1] a_all, cdef double [::1] variogram_model_parameters = np.asarray(pars['variogram_model_parameters']) - cdef double [::1,:] a_inv = np.asfortranarray(scipy.linalg.pinvh(a_all)) + cdef double [::1,:] a_inv = np.empty_like(a_all) + + + if pars['pseudo_inv']: + if pars['pseudo_inv_type'] == 1: + a_inv = np.asfortranarray(scipy.linalg.pinv(a_all)) + elif pars['pseudo_inv_type'] == 2: + a_inv = np.asfortranarray(scipy.linalg.pinv2(a_all)) + elif pars['pseudo_inv_type'] == 3: + a_inv = np.asfortranarray(scipy.linalg.pinvh(a_all)) + else: + raise ValueError('Unknown pseudo inverse method selected.') + else: + a_inv = scipy.linalg.inv(a_all) + for i in range(npt): # same thing as range(npt) if mask is not defined, otherwise take the non masked elements if mask[i]: diff --git a/pykrige/ok.py b/pykrige/ok.py index 92c152d..b51911d 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -15,7 +15,7 @@ ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. -.. [2] N. Cressie, Statistics for spatial data, +.. [2] N. Cressie, Statistics for spatial data, (Wiley Series in Probability and Statistics, 1993) 137 p. Copyright (c) 2015-2020, PyKrige Developers @@ -31,6 +31,7 @@ _initialize_variogram_model, _make_variogram_parameter_list, _find_statistics, + P_INV, ) import warnings @@ -143,19 +144,32 @@ class OrdinaryKriging: coordinates, both given in degree. Longitudes are expected in [0, 360] and latitudes in [-90, 90]. Default is 'euclidean'. exact_values : bool, optional - If True, interpolation provides input values at input locations. - If False, interpolation accounts for variance/nugget within input - values at input locations and does not behave as an + If True, interpolation provides input values at input locations. + If False, interpolation accounts for variance/nugget within input + values at input locations and does not behave as an exact-interpolator [2]. Note that this only has an effect if there is variance/nugget present within the input data since it is interpreted as measurement error. If the nugget is zero, the kriged field will behave as an exact interpolator. + pseudo_inv : :class:`bool`, optional + Whether the kriging system is solved with the pseudo inverted + kriging matrix. If `True`, this leads to more numerical stability + and redundant points are averaged. But it can take more time. + Default: False + pseudo_inv_type : :class:`int`, optional + Here you can select the algorithm to compute the pseudo-inverse matrix: + + * `1`: use `pinv` from `scipy` which uses `lstsq` + * `2`: use `pinv2` from `scipy` which uses `SVD` + * `3`: use `pinvh` from `scipy` which uses eigen-values + + Default: `1` References ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. - .. [2] N. Cressie, Statistics for spatial data, + .. [2] N. Cressie, Statistics for spatial data, (Wiley Series in Probability and Statistics, 1993) 137 p. """ @@ -186,7 +200,14 @@ def __init__( enable_statistics=False, coordinates_type="euclidean", exact_values=True, + pseudo_inv=False, + pseudo_inv_type=1, ): + # config the pseudo inverse + self.pseudo_inv = bool(pseudo_inv) + self.pseudo_inv_type = int(pseudo_inv_type) + if self.pseudo_inv_type not in [1, 2, 3]: + raise ValueError("pseudo inv type needs to be in [1,2,3]") # set up variogram model and parameters... self.variogram_model = variogram_model @@ -619,7 +640,11 @@ def _exec_vector(self, a, bd, mask): zero_index = None zero_value = False - a_inv = scipy.linalg.pinvh(a) + # use the desired method to invert the kriging matrix + if self.pseudo_inv: + a_inv = P_INV[self.pseudo_inv_type](a) + else: + a_inv = scipy.linalg.inv(a) if np.any(np.absolute(bd) <= self.eps): zero_value = True @@ -650,7 +675,11 @@ def _exec_loop(self, a, bd_all, mask): zvalues = np.zeros(npt) sigmasq = np.zeros(npt) - a_inv = scipy.linalg.pinvh(a) + # use the desired method to invert the kriging matrix + if self.pseudo_inv: + a_inv = P_INV[self.pseudo_inv_type](a) + else: + a_inv = scipy.linalg.inv(a) for j in np.nonzero(~mask)[ 0 @@ -876,6 +905,8 @@ def execute( "variogram_model_parameters", "variogram_function", "exact_values", + "pseudo_inv", + "pseudo_inv_type", ] } diff --git a/pykrige/ok3d.py b/pykrige/ok3d.py index a0404c8..4007f62 100644 --- a/pykrige/ok3d.py +++ b/pykrige/ok3d.py @@ -14,7 +14,7 @@ ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. -.. [2] N. Cressie, Statistics for spatial data, +.. [2] N. Cressie, Statistics for spatial data, (Wiley Series in Probability and Statistics, 1993) 137 p. Copyright (c) 2015-2020, PyKrige Developers @@ -30,6 +30,7 @@ _initialize_variogram_model, _make_variogram_parameter_list, _find_statistics, + P_INV, ) import warnings @@ -154,19 +155,32 @@ class OrdinaryKriging3D: enable_plotting : bool, optional Enables plotting to display variogram. Default is False (off). exact_values : bool, optional - If True, interpolation provides input values at input locations. - If False, interpolation accounts for variance/nugget within input - values at input locations and does not behave as an + If True, interpolation provides input values at input locations. + If False, interpolation accounts for variance/nugget within input + values at input locations and does not behave as an exact-interpolator [2]. Note that this only has an effect if there is variance/nugget present within the input data since it is interpreted as measurement error. If the nugget is zero, the kriged field will behave as an exact interpolator. + pseudo_inv : :class:`bool`, optional + Whether the kriging system is solved with the pseudo inverted + kriging matrix. If `True`, this leads to more numerical stability + and redundant points are averaged. But it can take more time. + Default: False + pseudo_inv_type : :class:`int`, optional + Here you can select the algorithm to compute the pseudo-inverse matrix: + + * `1`: use `pinv` from `scipy` which uses `lstsq` + * `2`: use `pinv2` from `scipy` which uses `SVD` + * `3`: use `pinvh` from `scipy` which uses eigen-values + + Default: `1` References ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. - .. [2] N. Cressie, Statistics for spatial data, + .. [2] N. Cressie, Statistics for spatial data, (Wiley Series in Probability and Statistics, 1993) 137 p. """ @@ -199,7 +213,14 @@ def __init__( verbose=False, enable_plotting=False, exact_values=True, + pseudo_inv=False, + pseudo_inv_type=1, ): + # config the pseudo inverse + self.pseudo_inv = bool(pseudo_inv) + self.pseudo_inv_type = int(pseudo_inv_type) + if self.pseudo_inv_type not in [1, 2, 3]: + raise ValueError("pseudo inv type needs to be in [1,2,3]") # set up variogram model and parameters... self.variogram_model = variogram_model @@ -604,7 +625,11 @@ def _exec_vector(self, a, bd, mask): zero_index = None zero_value = False - a_inv = scipy.linalg.pinvh(a) + # use the desired method to invert the kriging matrix + if self.pseudo_inv: + a_inv = P_INV[self.pseudo_inv_type](a) + else: + a_inv = scipy.linalg.inv(a) if np.any(np.absolute(bd) <= self.eps): zero_value = True @@ -635,7 +660,11 @@ def _exec_loop(self, a, bd_all, mask): kvalues = np.zeros(npt) sigmasq = np.zeros(npt) - a_inv = scipy.linalg.pinvh(a) + # use the desired method to invert the kriging matrix + if self.pseudo_inv: + a_inv = P_INV[self.pseudo_inv_type](a) + else: + a_inv = scipy.linalg.inv(a) for j in np.nonzero(~mask)[ 0 diff --git a/pykrige/uk.py b/pykrige/uk.py index 5eaba1f..59bdcb0 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -15,7 +15,7 @@ ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. -.. [2] N. Cressie, Statistics for spatial data, +.. [2] N. Cressie, Statistics for spatial data, (Wiley Series in Probability and Statistics, 1993) 137 p. Copyright (c) 2015-2020, PyKrige Developers """ @@ -29,6 +29,7 @@ _initialize_variogram_model, _make_variogram_parameter_list, _find_statistics, + P_INV, ) import warnings @@ -174,19 +175,32 @@ class UniversalKriging: enable_plotting : boolean, optional Enables plotting to display variogram. Default is False (off). exact_values : bool, optional - If True, interpolation provides input values at input locations. - If False, interpolation accounts for variance/nugget within input - values at input locations and does not behave as an + If True, interpolation provides input values at input locations. + If False, interpolation accounts for variance/nugget within input + values at input locations and does not behave as an exact-interpolator [2]. Note that this only has an effect if there is variance/nugget present within the input data since it is interpreted as measurement error. If the nugget is zero, the kriged field will behave as an exact interpolator. + pseudo_inv : :class:`bool`, optional + Whether the kriging system is solved with the pseudo inverted + kriging matrix. If `True`, this leads to more numerical stability + and redundant points are averaged. But it can take more time. + Default: False + pseudo_inv_type : :class:`int`, optional + Here you can select the algorithm to compute the pseudo-inverse matrix: + + * `1`: use `pinv` from `scipy` which uses `lstsq` + * `2`: use `pinv2` from `scipy` which uses `SVD` + * `3`: use `pinvh` from `scipy` which uses eigen-values + + Default: `1` References ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. - .. [2] N. Cressie, Statistics for spatial data, + .. [2] N. Cressie, Statistics for spatial data, (Wiley Series in Probability and Statistics, 1993) 137 p. """ @@ -224,7 +238,14 @@ def __init__( verbose=False, enable_plotting=False, exact_values=True, + pseudo_inv=False, + pseudo_inv_type=1, ): + # config the pseudo inverse + self.pseudo_inv = bool(pseudo_inv) + self.pseudo_inv_type = int(pseudo_inv_type) + if self.pseudo_inv_type not in [1, 2, 3]: + raise ValueError("pseudo inv type needs to be in [1,2,3]") # Deal with mutable default argument if drift_terms is None: @@ -895,7 +916,11 @@ def _exec_vector(self, a, bd, xy, xy_orig, mask, n_withdrifts, spec_drift_grids) zero_index = None zero_value = False - a_inv = scipy.linalg.pinvh(a) + # use the desired method to invert the kriging matrix + if self.pseudo_inv: + a_inv = P_INV[self.pseudo_inv_type](a) + else: + a_inv = scipy.linalg.inv(a) if np.any(np.absolute(bd) <= self.eps): zero_value = True @@ -979,7 +1004,11 @@ def _exec_loop(self, a, bd_all, xy, xy_orig, mask, n_withdrifts, spec_drift_grid zvalues = np.zeros(npt) sigmasq = np.zeros(npt) - a_inv = scipy.linalg.pinvh(a) + # use the desired method to invert the kriging matrix + if self.pseudo_inv: + a_inv = P_INV[self.pseudo_inv_type](a) + else: + a_inv = scipy.linalg.inv(a) for j in np.nonzero(~mask)[ 0 diff --git a/pykrige/uk3d.py b/pykrige/uk3d.py index 1a64c2b..cc05f5e 100644 --- a/pykrige/uk3d.py +++ b/pykrige/uk3d.py @@ -14,7 +14,7 @@ ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. -.. [2] N. Cressie, Statistics for spatial data, +.. [2] N. Cressie, Statistics for spatial data, (Wiley Series in Probability and Statistics, 1993) 137 p. Copyright (c) 2015-2020, PyKrige Developers @@ -29,6 +29,7 @@ _initialize_variogram_model, _make_variogram_parameter_list, _find_statistics, + P_INV, ) import warnings @@ -169,19 +170,32 @@ class UniversalKriging3D: enable_plotting : boolean, optional Enables plotting to display variogram. Default is False (off). exact_values : bool, optional - If True, interpolation provides input values at input locations. - If False, interpolation accounts for variance/nugget within input - values at input locations and does not behave as an + If True, interpolation provides input values at input locations. + If False, interpolation accounts for variance/nugget within input + values at input locations and does not behave as an exact-interpolator [2]. Note that this only has an effect if there is variance/nugget present within the input data since it is interpreted as measurement error. If the nugget is zero, the kriged field will behave as an exact interpolator. + pseudo_inv : :class:`bool`, optional + Whether the kriging system is solved with the pseudo inverted + kriging matrix. If `True`, this leads to more numerical stability + and redundant points are averaged. But it can take more time. + Default: False + pseudo_inv_type : :class:`int`, optional + Here you can select the algorithm to compute the pseudo-inverse matrix: + + * `1`: use `pinv` from `scipy` which uses `lstsq` + * `2`: use `pinv2` from `scipy` which uses `SVD` + * `3`: use `pinvh` from `scipy` which uses eigen-values + + Default: `1` References ---------- .. [1] P.K. Kitanidis, Introduction to Geostatistcs: Applications in Hydrogeology, (Cambridge University Press, 1997) 272 p. - .. [2] N. Cressie, Statistics for spatial data, + .. [2] N. Cressie, Statistics for spatial data, (Wiley Series in Probability and Statistics, 1993) 137 p. """ @@ -219,7 +233,14 @@ def __init__( verbose=False, enable_plotting=False, exact_values=True, + pseudo_inv=False, + pseudo_inv_type=1, ): + # config the pseudo inverse + self.pseudo_inv = bool(pseudo_inv) + self.pseudo_inv_type = int(pseudo_inv_type) + if self.pseudo_inv_type not in [1, 2, 3]: + raise ValueError("pseudo inv type needs to be in [1,2,3]") # Deal with mutable default argument if drift_terms is None: @@ -719,7 +740,11 @@ def _exec_vector(self, a, bd, xyz, mask, n_withdrifts, spec_drift_grids): zero_index = None zero_value = False - a_inv = scipy.linalg.pinvh(a) + # use the desired method to invert the kriging matrix + if self.pseudo_inv: + a_inv = P_INV[self.pseudo_inv_type](a) + else: + a_inv = scipy.linalg.inv(a) if np.any(np.absolute(bd) <= self.eps): zero_value = True @@ -788,7 +813,11 @@ def _exec_loop(self, a, bd_all, xyz, mask, n_withdrifts, spec_drift_grids): kvalues = np.zeros(npt) sigmasq = np.zeros(npt) - a_inv = scipy.linalg.pinvh(a) + # use the desired method to invert the kriging matrix + if self.pseudo_inv: + a_inv = P_INV[self.pseudo_inv_type](a) + else: + a_inv = scipy.linalg.inv(a) for j in np.nonzero(~mask)[ 0 From ac641c7413fe94cb59c0a3d67e212f67e4bdf267 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 15 Jul 2020 19:39:53 +0200 Subject: [PATCH 25/41] cython: assure inv matrix to be fortran contiguous --- pykrige/lib/cok.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pykrige/lib/cok.pyx b/pykrige/lib/cok.pyx index 80e0b21..b3b187e 100644 --- a/pykrige/lib/cok.pyx +++ b/pykrige/lib/cok.pyx @@ -35,7 +35,7 @@ cpdef _c_exec_loop(double [:, ::1] a_all, cdef double [::1] variogram_model_parameters = np.asarray(pars['variogram_model_parameters']) - cdef double [::1,:] a_inv = np.empty_like(a_all) + cdef double [::1,:] a_inv = np.asfortranarray(np.empty_like(a_all)) if pars['pseudo_inv']: @@ -48,7 +48,7 @@ cpdef _c_exec_loop(double [:, ::1] a_all, else: raise ValueError('Unknown pseudo inverse method selected.') else: - a_inv = scipy.linalg.inv(a_all) + a_inv = np.asfortranarray(scipy.linalg.inv(a_all)) for i in range(npt): # same thing as range(npt) if mask is not defined, otherwise take the non masked elements From 198394913a376692a834558ff67e85118273c93b Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 15 Jul 2020 20:05:46 +0200 Subject: [PATCH 26/41] Test: check if usage of pseudo inv results in mean values for redundant points --- tests/test_core.py | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/tests/test_core.py b/tests/test_core.py index 5def740..4ea1af4 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -2891,3 +2891,42 @@ def test_gstools_variogram(model): y_id = int(data[i, 1] * 10) x_id = int(data[i, 0] * 10) assert np.isclose(z1[y_id, x_id], data[i, 2]) + + +@pytest.mark.parametrize("model", [OrdinaryKriging, UniversalKriging]) +def test_pseudo_2d(model): + # test data + data = np.array([[0.0, 0.0, 1.0], [0.0, 0.0, 3.0], [1.0, 0.0, 6.0],]) + for p_type in [1, 2, 3]: + # create the krige field + krige = model( + data[:, 0], + data[:, 1], + data[:, 2], + variogram_parameters=[1.0, 0.0], + pseudo_inv=True, + pseudo_inv_type=p_type, + ) + z1, ss1 = krige.execute("points", 0.0, 0.0) + # check if the field coincides with the mean of the redundant data + assert np.isclose(z1.item(), 2.0) + + +@pytest.mark.parametrize("model", [OrdinaryKriging3D, UniversalKriging3D]) +def test_pseudo_3d(model): + # test data + data = np.array([[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 3.0], [1.0, 0.0, 0.0, 6.0],]) + for p_type in [1, 2, 3]: + # create the krige field + krige = model( + data[:, 0], + data[:, 1], + data[:, 2], + data[:, 3], + variogram_parameters=[1.0, 0.0], + pseudo_inv=True, + pseudo_inv_type=p_type, + ) + z1, ss1 = krige.execute("points", 0.0, 0.0, 0.0) + # check if the field coincides with the mean of the redundant data + assert np.isclose(z1.item(), 2.0) From a32bc2b45300e7f49179ee7e0d0358068b5eaf9d Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 16 Jul 2020 00:28:39 +0200 Subject: [PATCH 27/41] add lstsq for statistics when pseudo_inv in use --- pykrige/core.py | 31 ++++++++++++++++++++++++++++--- pykrige/ok.py | 2 ++ pykrige/ok3d.py | 2 ++ pykrige/uk.py | 2 ++ pykrige/uk3d.py | 2 ++ tests/test_core.py | 4 ++-- 6 files changed, 38 insertions(+), 5 deletions(-) diff --git a/pykrige/core.py b/pykrige/core.py index 1e02b02..67d7ce5 100755 --- a/pykrige/core.py +++ b/pykrige/core.py @@ -679,7 +679,13 @@ def _calculate_variogram_model( def _krige( - X, y, coords, variogram_function, variogram_model_parameters, coordinates_type + X, + y, + coords, + variogram_function, + variogram_model_parameters, + coordinates_type, + pseudo_inv=False, ): """Sets up and solves the ordinary kriging system for the given coordinate pair. This function is only used for the statistics calculations. @@ -699,6 +705,11 @@ def _krige( coordinates_type: str type of coordinates in X array, can be 'euclidean' for standard rectangular coordinates or 'geographic' if the coordinates are lat/lon + pseudo_inv : :class:`bool`, optional + Whether the kriging system is solved with the pseudo inverted + kriging matrix. If `True`, this leads to more numerical stability + and redundant points are averaged. But it can take more time. + Default: False Returns ------- @@ -760,7 +771,10 @@ def _krige( b[n, 0] = 1.0 # solve - res = np.linalg.solve(a, b) + if pseudo_inv: + res = np.linalg.lstsq(a, b) + else: + res = np.linalg.solve(a, b) zinterp = np.sum(res[:n, 0] * y) sigmasq = np.sum(res[:, 0] * -b[:, 0]) @@ -768,7 +782,12 @@ def _krige( def _find_statistics( - X, y, variogram_function, variogram_model_parameters, coordinates_type + X, + y, + variogram_function, + variogram_model_parameters, + coordinates_type, + pseudo_inv=False, ): """Calculates variogram fit statistics. Returns the delta, sigma, and epsilon values for the variogram fit. @@ -787,6 +806,11 @@ def _find_statistics( coordinates_type: str type of coordinates in X array, can be 'euclidean' for standard rectangular coordinates or 'geographic' if the coordinates are lat/lon + pseudo_inv : :class:`bool`, optional + Whether the kriging system is solved with the pseudo inverted + kriging matrix. If `True`, this leads to more numerical stability + and redundant points are averaged. But it can take more time. + Default: False Returns ------- @@ -815,6 +839,7 @@ def _find_statistics( variogram_function, variogram_model_parameters, coordinates_type, + pseudo_inv, ) # if the estimation error is zero, it's probably because diff --git a/pykrige/ok.py b/pykrige/ok.py index b51911d..9ef63fa 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -356,6 +356,7 @@ def __init__( self.variogram_function, self.variogram_model_parameters, self.coordinates_type, + self.pseudo_inv, ) self.Q1 = core.calcQ1(self.epsilon) self.Q2 = core.calcQ2(self.epsilon) @@ -527,6 +528,7 @@ def update_variogram_model( self.variogram_function, self.variogram_model_parameters, self.coordinates_type, + self.pseudo_inv, ) self.Q1 = core.calcQ1(self.epsilon) self.Q2 = core.calcQ2(self.epsilon) diff --git a/pykrige/ok3d.py b/pykrige/ok3d.py index 4007f62..1a92ea2 100644 --- a/pykrige/ok3d.py +++ b/pykrige/ok3d.py @@ -353,6 +353,7 @@ def __init__( self.variogram_function, self.variogram_model_parameters, "euclidean", + self.pseudo_inv, ) self.Q1 = core.calcQ1(self.epsilon) self.Q2 = core.calcQ2(self.epsilon) @@ -535,6 +536,7 @@ def update_variogram_model( self.variogram_function, self.variogram_model_parameters, "euclidean", + self.pseudo_inv, ) self.Q1 = core.calcQ1(self.epsilon) self.Q2 = core.calcQ2(self.epsilon) diff --git a/pykrige/uk.py b/pykrige/uk.py index 59bdcb0..6d12998 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -377,6 +377,7 @@ def __init__( self.variogram_function, self.variogram_model_parameters, "euclidean", + self.pseudo_inv, ) self.Q1 = core.calcQ1(self.epsilon) self.Q2 = core.calcQ2(self.epsilon) @@ -766,6 +767,7 @@ def update_variogram_model( self.variogram_function, self.variogram_model_parameters, "euclidean", + self.pseudo_inv, ) self.Q1 = core.calcQ1(self.epsilon) self.Q2 = core.calcQ2(self.epsilon) diff --git a/pykrige/uk3d.py b/pykrige/uk3d.py index cc05f5e..56894e8 100644 --- a/pykrige/uk3d.py +++ b/pykrige/uk3d.py @@ -381,6 +381,7 @@ def __init__( self.variogram_function, self.variogram_model_parameters, "euclidean", + self.pseudo_inv, ) self.Q1 = core.calcQ1(self.epsilon) self.Q2 = core.calcQ2(self.epsilon) @@ -620,6 +621,7 @@ def update_variogram_model( self.variogram_function, self.variogram_model_parameters, "euclidean", + self.pseudo_inv, ) self.Q1 = core.calcQ1(self.epsilon) self.Q2 = core.calcQ2(self.epsilon) diff --git a/tests/test_core.py b/tests/test_core.py index 4ea1af4..d251c32 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -2896,7 +2896,7 @@ def test_gstools_variogram(model): @pytest.mark.parametrize("model", [OrdinaryKriging, UniversalKriging]) def test_pseudo_2d(model): # test data - data = np.array([[0.0, 0.0, 1.0], [0.0, 0.0, 3.0], [1.0, 0.0, 6.0],]) + data = np.array([[0.0, 0.0, 1.0], [0.0, 0.0, 3.0], [1.0, 0.0, 6.0]]) for p_type in [1, 2, 3]: # create the krige field krige = model( @@ -2915,7 +2915,7 @@ def test_pseudo_2d(model): @pytest.mark.parametrize("model", [OrdinaryKriging3D, UniversalKriging3D]) def test_pseudo_3d(model): # test data - data = np.array([[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 3.0], [1.0, 0.0, 0.0, 6.0],]) + data = np.array([[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 3.0], [1.0, 0.0, 0.0, 6.0]]) for p_type in [1, 2, 3]: # create the krige field krige = model( From 23b8e40996631d45e8bb140e94e45cf4b58327bb Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Thu, 16 Jul 2020 11:05:30 +0200 Subject: [PATCH 28/41] core: correctly use lstsq --- pykrige/core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pykrige/core.py b/pykrige/core.py index 67d7ce5..ba8ecfc 100755 --- a/pykrige/core.py +++ b/pykrige/core.py @@ -772,7 +772,7 @@ def _krige( # solve if pseudo_inv: - res = np.linalg.lstsq(a, b) + res = np.linalg.lstsq(a, b, rcond=None)[0] else: res = np.linalg.solve(a, b) zinterp = np.sum(res[:n, 0] * y) From a00297a376f437ae9b17afc85c86b53d1e893c90 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Sebastian=20M=C3=BCller?= Date: Fri, 7 Aug 2020 20:53:36 +0200 Subject: [PATCH 29/41] Grid search update (#158) --- README.rst | 2 +- examples/krige_cv.py | 2 +- pykrige/rk.py | 295 +++++++++++++++++++++++++++++++++---------- 3 files changed, 227 insertions(+), 72 deletions(-) diff --git a/README.rst b/README.rst index 8673739..9812510 100644 --- a/README.rst +++ b/README.rst @@ -12,7 +12,7 @@ PyKrige .. image:: https://coveralls.io/repos/github/GeoStat-Framework/PyKrige/badge.svg?branch=master :target: https://coveralls.io/github/GeoStat-Framework/PyKrige?branch=master .. image:: https://readthedocs.org/projects/pykrige/badge/?version=stable - :target: http://pykrige.readthedocs.io/en/latest/?badge=stable + :target: http://pykrige.readthedocs.io/en/stable/?badge=stable :alt: Documentation Status .. image:: https://img.shields.io/badge/code%20style-black-000000.svg :target: https://github.com/psf/black diff --git a/examples/krige_cv.py b/examples/krige_cv.py index 0b967c5..ece8fd8 100644 --- a/examples/krige_cv.py +++ b/examples/krige_cv.py @@ -53,7 +53,7 @@ # "weight": [True, False] } -estimator = GridSearchCV(Krige(), param_dict3d, verbose=True) +estimator = GridSearchCV(Krige(), param_dict3d, verbose=True, return_train_score=True) # dummy data X3 = np.random.randint(0, 400, size=(100, 3)).astype(float) diff --git a/pykrige/rk.py b/pykrige/rk.py index 1f08ff2..2b3b40f 100644 --- a/pykrige/rk.py +++ b/pykrige/rk.py @@ -19,8 +19,44 @@ threed_krige = ("ordinary3d", "universal3d") +krige_methods_kws = { + "ordinary": [ + "anisotropy_scaling", + "anisotropy_angle", + "enable_statistics", + "coordinates_type", + ], + "universal": [ + "anisotropy_scaling", + "anisotropy_angle", + "drift_terms", + "point_drift", + "external_drift", + "external_drift_x", + "external_drift_y", + "functional_drift", + ], + "ordinary3d": [ + "anisotropy_scaling_y", + "anisotropy_scaling_z", + "anisotropy_angle_x", + "anisotropy_angle_y", + "anisotropy_angle_z", + ], + "universal3d": [ + "anisotropy_scaling_y", + "anisotropy_scaling_z", + "anisotropy_angle_x", + "anisotropy_angle_y", + "anisotropy_angle_z", + "drift_terms", + "functional_drift", + ], +} + def validate_method(method): + """Validate the kriging method in use.""" if method not in krige_methods.keys(): raise ValueError( "Kriging method must be one of {}".format(krige_methods.keys()) @@ -30,9 +66,47 @@ def validate_method(method): class Krige(RegressorMixin, BaseEstimator): """ A scikit-learn wrapper class for Ordinary and Universal Kriging. + This works with both Grid/RandomSearchCv for finding the best Krige parameters combination for a problem. + Parameters + ---------- + method: str, optional + type of kriging to be performed + variogram_model: str, optional + variogram model to be used during Kriging + nlags: int + see OK/UK class description + weight: bool + see OK/UK class description + n_closest_points: int + number of closest points to be used during Ordinary Kriging + verbose: bool + see OK/UK class description + exact_values : bool + see OK/UK class description + variogram_parameters : list or dict + see OK/UK class description + variogram_function : callable + see OK/UK class description + anisotropy_scaling : tuple + single value for 2D (UK/OK) and two values in 3D (UK3D/OK3D) + anisotropy_angle : tuple + single value for 2D (UK/OK) and three values in 3D (UK3D/OK3D) + enable_statistics : bool + see OK class description + coordinates_type : str + see OK/UK class description + drift_terms : list of strings + see UK/UK3D class description + point_drift : array_like + see UK class description + ext_drift_grid : tuple + Holding the three values external_drift, external_drift_x and + external_drift_z for the UK class + functional_drift : list of callable + see UK/UK3D class description """ def __init__( @@ -43,19 +117,42 @@ def __init__( weight=False, n_closest_points=10, verbose=False, + exact_values=True, + variogram_parameters=None, + variogram_function=None, + anisotropy_scaling=(1.0, 1.0), + anisotropy_angle=(0.0, 0.0, 0.0), + enable_statistics=False, + coordinates_type="euclidean", + drift_terms=None, + point_drift=None, + ext_drift_grid=(None, None, None), + functional_drift=None, ): - validate_method(method) self.variogram_model = variogram_model - self.verbose = verbose + self.variogram_parameters = variogram_parameters + self.variogram_function = variogram_function self.nlags = nlags self.weight = weight + self.verbose = verbose + self.exact_values = exact_values + self.anisotropy_scaling = anisotropy_scaling + self.anisotropy_angle = anisotropy_angle + self.enable_statistics = enable_statistics + self.coordinates_type = coordinates_type + self.drift_terms = drift_terms + self.point_drift = point_drift + self.ext_drift_grid = ext_drift_grid + self.functional_drift = functional_drift self.model = None # not trained self.n_closest_points = n_closest_points self.method = method def fit(self, x, y, *args, **kwargs): """ + Fit the current model. + Parameters ---------- x: ndarray @@ -64,28 +161,39 @@ def fit(self, x, y, *args, **kwargs): y: ndarray array of targets (N, ) """ - - points = self._dimensionality_check(x) - - # if condition required to address backward compatibility - if self.method in threed_krige: - self.model = krige_methods[self.method]( - val=y, - variogram_model=self.variogram_model, - nlags=self.nlags, - weight=self.weight, - verbose=self.verbose, - **points - ) - else: - self.model = krige_methods[self.method]( - z=y, - variogram_model=self.variogram_model, - nlags=self.nlags, - weight=self.weight, - verbose=self.verbose, - **points - ) + val_kw = "val" if self.method in threed_krige else "z" + setup = dict( + variogram_model=self.variogram_model, + variogram_parameters=self.variogram_parameters, + variogram_function=self.variogram_function, + nlags=self.nlags, + weight=self.weight, + verbose=self.verbose, + exact_values=self.exact_values, + ) + add_setup = dict( + anisotropy_scaling=self.anisotropy_scaling[0], + anisotropy_angle=self.anisotropy_angle[0], + enable_statistics=self.enable_statistics, + coordinates_type=self.coordinates_type, + anisotropy_scaling_y=self.anisotropy_scaling[0], + anisotropy_scaling_z=self.anisotropy_scaling[1], + anisotropy_angle_x=self.anisotropy_angle[0], + anisotropy_angle_y=self.anisotropy_angle[1], + anisotropy_angle_z=self.anisotropy_angle[2], + drift_terms=self.drift_terms, + point_drift=self.point_drift, + external_drift=self.ext_drift_grid[0], + external_drift_x=self.ext_drift_grid[1], + external_drift_y=self.ext_drift_grid[2], + functional_drift=self.functional_drift, + ) + for kw in krige_methods_kws[self.method]: + setup[kw] = add_setup[kw] + input_kw = self._dimensionality_check(x) + input_kw.update(setup) + input_kw[val_kw] = y + self.model = krige_methods[self.method](**input_kw) def _dimensionality_check(self, x, ext=""): if self.method in ("ordinary", "universal"): @@ -97,58 +205,57 @@ def _dimensionality_check(self, x, ext=""): if x.shape[1] != 3: raise ValueError("3d krige can use only 3d points") else: - return {"x" + ext: x[:, 0], "y" + ext: x[:, 1], "z" + ext: x[:, 2]} + return { + "x" + ext: x[:, 0], + "y" + ext: x[:, 1], + "z" + ext: x[:, 2], + } def predict(self, x, *args, **kwargs): """ + Predict. + Parameters ---------- x: ndarray array of Points, (x, y) pairs of shape (N, 2) for 2d kriging array of Points, (x, y, z) pairs of shape (N, 3) for 3d kriging - Returns ------- Prediction array """ if not self.model: raise Exception("Not trained. Train first") - points = self._dimensionality_check(x, ext="points") - return self.execute(points, *args, **kwargs)[0] def execute(self, points, *args, **kwargs): # TODO array of Points, (x, y) pairs of shape (N, 2) """ + Execute. + Parameters ---------- points: dict - Returns: + Returns ------- Prediction array Variance array """ - if isinstance(self.model, OrdinaryKriging) or isinstance( - self.model, OrdinaryKriging3D - ): - prediction, variance = self.model.execute( - "points", - n_closest_points=self.n_closest_points, - backend="loop", - **points - ) + default_kw = dict(style="points", backend="loop") + default_kw.update(kwargs) + points.update(default_kw) + if isinstance(self.model, (OrdinaryKriging, OrdinaryKriging3D)): + points.update(dict(n_closest_points=self.n_closest_points)) else: print("n_closest_points will be ignored for UniversalKriging") - prediction, variance = self.model.execute( - "points", backend="loop", **points - ) - + prediction, variance = self.model.execute(**points) return prediction, variance def check_sklearn_model(model): + """Check the sklearn method in use.""" if not (isinstance(model, BaseEstimator) and isinstance(model, RegressorMixin)): raise RuntimeError( "Needs to supply an instance of a scikit-learn regression class." @@ -157,8 +264,49 @@ def check_sklearn_model(model): class RegressionKriging: """ - This is an implementation of Regression-Kriging as described here: + An implementation of Regression-Kriging. + + As described here: https://en.wikipedia.org/wiki/Regression-Kriging + + Parameters + ---------- + regression_model: machine learning model instance from sklearn + method: str, optional + type of kriging to be performed + variogram_model: str, optional + variogram model to be used during Kriging + n_closest_points: int + number of closest points to be used during Ordinary Kriging + nlags: int + see OK/UK class description + weight: bool + see OK/UK class description + verbose: bool + see OK/UK class description + exact_values : bool + see OK/UK class description + variogram_parameters : list or dict + see OK/UK class description + variogram_function : callable + see OK/UK class description + anisotropy_scaling : tuple + single value for 2D (UK/OK) and two values in 3D (UK3D/OK3D) + anisotropy_angle : tuple + single value for 2D (UK/OK) and three values in 3D (UK3D/OK3D) + enable_statistics : bool + see OK class description + coordinates_type : str + see OK/UK class description + drift_terms : list of strings + see UK/UK3D class description + point_drift : array_like + see UK class description + ext_drift_grid : tuple + Holding the three values external_drift, external_drift_x and + external_drift_z for the UK class + functional_drift : list of callable + see UK/UK3D class description """ def __init__( @@ -170,24 +318,18 @@ def __init__( nlags=6, weight=False, verbose=False, + exact_values=True, + variogram_parameters=None, + variogram_function=None, + anisotropy_scaling=(1.0, 1.0), + anisotropy_angle=(0.0, 0.0, 0.0), + enable_statistics=False, + coordinates_type="euclidean", + drift_terms=None, + point_drift=None, + ext_drift_grid=(None, None, None), + functional_drift=None, ): - """ - Parameters - ---------- - regression_model: machine learning model instance from sklearn - method: str, optional - type of kriging to be performed - variogram_model: str, optional - variogram model to be used during Kriging - n_closest_points: int - number of closest points to be used during Ordinary Kriging - nlags: int - see OK/UK class description - weight: bool - see OK/UK class description - verbose: bool - see OK/UK class description - """ check_sklearn_model(regression_model) self.regression_model = regression_model self.n_closest_points = n_closest_points @@ -198,11 +340,22 @@ def __init__( weight=weight, n_closest_points=n_closest_points, verbose=verbose, + exact_values=exact_values, + variogram_parameters=variogram_parameters, + variogram_function=variogram_function, + anisotropy_scaling=anisotropy_scaling, + anisotropy_angle=anisotropy_angle, + enable_statistics=enable_statistics, + coordinates_type=coordinates_type, + drift_terms=drift_terms, + point_drift=point_drift, + ext_drift_grid=ext_drift_grid, + functional_drift=functional_drift, ) def fit(self, p, x, y): """ - fit the regression method and also Krige the residual + Fit the regression method and also Krige the residual. Parameters ---------- @@ -223,8 +376,10 @@ def fit(self, p, x, y): self.krige.fit(x=x, y=y - ml_pred) print("Finished kriging residuals") - def predict(self, p, x): + def predict(self, p, x, **kwargs): """ + Predict. + Parameters ---------- p: ndarray @@ -241,11 +396,12 @@ def predict(self, p, x): The expected value of ys for the query inputs, of shape (Ns,). """ + return self.krige_residual(x, **kwargs) + self.regression_model.predict(p) - return self.krige_residual(x) + self.regression_model.predict(p) - - def krige_residual(self, x): + def krige_residual(self, x, **kwargs): """ + Calculate the residuals. + Parameters ---------- x: ndarray @@ -257,11 +413,11 @@ def krige_residual(self, x): residual: ndarray kriged residual values """ - return self.krige.predict(x) + return self.krige.predict(x, **kwargs) - def score(self, p, x, y, sample_weight=None): + def score(self, p, x, y, sample_weight=None, **kwargs): """ - Overloading default regression score method + Overloading default regression score method. Parameters ---------- @@ -275,7 +431,6 @@ def score(self, p, x, y, sample_weight=None): y: ndarray array of targets (Ns, ) """ - return r2_score( - y_pred=self.predict(p, x), y_true=y, sample_weight=sample_weight + y_pred=self.predict(p, x, **kwargs), y_true=y, sample_weight=sample_weight ) From 48a6d2c919ac67e6d3733209fb666d852039adf2 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 18 Aug 2020 14:15:55 +0200 Subject: [PATCH 30/41] pseudo-inv: change p-inv selector to strings --- pykrige/core.py | 2 +- pykrige/lib/cok.pyx | 6 +++--- pykrige/ok.py | 18 +++++++++--------- pykrige/ok3d.py | 18 +++++++++--------- pykrige/uk.py | 18 +++++++++--------- pykrige/uk3d.py | 18 +++++++++--------- 6 files changed, 40 insertions(+), 40 deletions(-) diff --git a/pykrige/core.py b/pykrige/core.py index ba8ecfc..02d8b4f 100755 --- a/pykrige/core.py +++ b/pykrige/core.py @@ -31,7 +31,7 @@ eps = 1.0e-10 # Cutoff for comparison to zero -P_INV = {1: spl.pinv, 2: spl.pinv2, 3: spl.pinvh} +P_INV = {"pinv": spl.pinv, "pinv2": spl.pinv2, "pinvh": spl.pinvh} def great_circle_distance(lon1, lat1, lon2, lat2): diff --git a/pykrige/lib/cok.pyx b/pykrige/lib/cok.pyx index b3b187e..9044c01 100644 --- a/pykrige/lib/cok.pyx +++ b/pykrige/lib/cok.pyx @@ -39,11 +39,11 @@ cpdef _c_exec_loop(double [:, ::1] a_all, if pars['pseudo_inv']: - if pars['pseudo_inv_type'] == 1: + if pars['pseudo_inv_type'] == "pinv": a_inv = np.asfortranarray(scipy.linalg.pinv(a_all)) - elif pars['pseudo_inv_type'] == 2: + elif pars['pseudo_inv_type'] == "pinv2": a_inv = np.asfortranarray(scipy.linalg.pinv2(a_all)) - elif pars['pseudo_inv_type'] == 3: + elif pars['pseudo_inv_type'] == "pinvh": a_inv = np.asfortranarray(scipy.linalg.pinvh(a_all)) else: raise ValueError('Unknown pseudo inverse method selected.') diff --git a/pykrige/ok.py b/pykrige/ok.py index 9ef63fa..358520d 100644 --- a/pykrige/ok.py +++ b/pykrige/ok.py @@ -156,14 +156,14 @@ class OrdinaryKriging: kriging matrix. If `True`, this leads to more numerical stability and redundant points are averaged. But it can take more time. Default: False - pseudo_inv_type : :class:`int`, optional + pseudo_inv_type : :class:`str`, optional Here you can select the algorithm to compute the pseudo-inverse matrix: - * `1`: use `pinv` from `scipy` which uses `lstsq` - * `2`: use `pinv2` from `scipy` which uses `SVD` - * `3`: use `pinvh` from `scipy` which uses eigen-values + * `"pinv"`: use `pinv` from `scipy` which uses `lstsq` + * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD` + * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values - Default: `1` + Default: `"pinv"` References ---------- @@ -201,13 +201,13 @@ def __init__( coordinates_type="euclidean", exact_values=True, pseudo_inv=False, - pseudo_inv_type=1, + pseudo_inv_type="pinv", ): # config the pseudo inverse self.pseudo_inv = bool(pseudo_inv) - self.pseudo_inv_type = int(pseudo_inv_type) - if self.pseudo_inv_type not in [1, 2, 3]: - raise ValueError("pseudo inv type needs to be in [1,2,3]") + self.pseudo_inv_type = str(pseudo_inv_type) + if self.pseudo_inv_type not in P_INV: + raise ValueError("pseudo inv type not valid: " + str(pseudo_inv_type)) # set up variogram model and parameters... self.variogram_model = variogram_model diff --git a/pykrige/ok3d.py b/pykrige/ok3d.py index 1a92ea2..5f52746 100644 --- a/pykrige/ok3d.py +++ b/pykrige/ok3d.py @@ -167,14 +167,14 @@ class OrdinaryKriging3D: kriging matrix. If `True`, this leads to more numerical stability and redundant points are averaged. But it can take more time. Default: False - pseudo_inv_type : :class:`int`, optional + pseudo_inv_type : :class:`str`, optional Here you can select the algorithm to compute the pseudo-inverse matrix: - * `1`: use `pinv` from `scipy` which uses `lstsq` - * `2`: use `pinv2` from `scipy` which uses `SVD` - * `3`: use `pinvh` from `scipy` which uses eigen-values + * `"pinv"`: use `pinv` from `scipy` which uses `lstsq` + * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD` + * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values - Default: `1` + Default: `"pinv"` References ---------- @@ -214,13 +214,13 @@ def __init__( enable_plotting=False, exact_values=True, pseudo_inv=False, - pseudo_inv_type=1, + pseudo_inv_type="pinv", ): # config the pseudo inverse self.pseudo_inv = bool(pseudo_inv) - self.pseudo_inv_type = int(pseudo_inv_type) - if self.pseudo_inv_type not in [1, 2, 3]: - raise ValueError("pseudo inv type needs to be in [1,2,3]") + self.pseudo_inv_type = str(pseudo_inv_type) + if self.pseudo_inv_type not in P_INV: + raise ValueError("pseudo inv type not valid: " + str(pseudo_inv_type)) # set up variogram model and parameters... self.variogram_model = variogram_model diff --git a/pykrige/uk.py b/pykrige/uk.py index 6d12998..707cf34 100644 --- a/pykrige/uk.py +++ b/pykrige/uk.py @@ -187,14 +187,14 @@ class UniversalKriging: kriging matrix. If `True`, this leads to more numerical stability and redundant points are averaged. But it can take more time. Default: False - pseudo_inv_type : :class:`int`, optional + pseudo_inv_type : :class:`str`, optional Here you can select the algorithm to compute the pseudo-inverse matrix: - * `1`: use `pinv` from `scipy` which uses `lstsq` - * `2`: use `pinv2` from `scipy` which uses `SVD` - * `3`: use `pinvh` from `scipy` which uses eigen-values + * `"pinv"`: use `pinv` from `scipy` which uses `lstsq` + * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD` + * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values - Default: `1` + Default: `"pinv"` References ---------- @@ -239,13 +239,13 @@ def __init__( enable_plotting=False, exact_values=True, pseudo_inv=False, - pseudo_inv_type=1, + pseudo_inv_type="pinv", ): # config the pseudo inverse self.pseudo_inv = bool(pseudo_inv) - self.pseudo_inv_type = int(pseudo_inv_type) - if self.pseudo_inv_type not in [1, 2, 3]: - raise ValueError("pseudo inv type needs to be in [1,2,3]") + self.pseudo_inv_type = str(pseudo_inv_type) + if self.pseudo_inv_type not in P_INV: + raise ValueError("pseudo inv type not valid: " + str(pseudo_inv_type)) # Deal with mutable default argument if drift_terms is None: diff --git a/pykrige/uk3d.py b/pykrige/uk3d.py index 56894e8..3fe8813 100644 --- a/pykrige/uk3d.py +++ b/pykrige/uk3d.py @@ -182,14 +182,14 @@ class UniversalKriging3D: kriging matrix. If `True`, this leads to more numerical stability and redundant points are averaged. But it can take more time. Default: False - pseudo_inv_type : :class:`int`, optional + pseudo_inv_type : :class:`str`, optional Here you can select the algorithm to compute the pseudo-inverse matrix: - * `1`: use `pinv` from `scipy` which uses `lstsq` - * `2`: use `pinv2` from `scipy` which uses `SVD` - * `3`: use `pinvh` from `scipy` which uses eigen-values + * `"pinv"`: use `pinv` from `scipy` which uses `lstsq` + * `"pinv2"`: use `pinv2` from `scipy` which uses `SVD` + * `"pinvh"`: use `pinvh` from `scipy` which uses eigen-values - Default: `1` + Default: `"pinv"` References ---------- @@ -234,13 +234,13 @@ def __init__( enable_plotting=False, exact_values=True, pseudo_inv=False, - pseudo_inv_type=1, + pseudo_inv_type="pinv", ): # config the pseudo inverse self.pseudo_inv = bool(pseudo_inv) - self.pseudo_inv_type = int(pseudo_inv_type) - if self.pseudo_inv_type not in [1, 2, 3]: - raise ValueError("pseudo inv type needs to be in [1,2,3]") + self.pseudo_inv_type = str(pseudo_inv_type) + if self.pseudo_inv_type not in P_INV: + raise ValueError("pseudo inv type not valid: " + str(pseudo_inv_type)) # Deal with mutable default argument if drift_terms is None: From c3217cf19d1e2926e7968f4e1b68776edaf2e649 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 18 Aug 2020 14:17:04 +0200 Subject: [PATCH 31/41] pseudo-inv: adopt tests to new selector type --- tests/test_core.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_core.py b/tests/test_core.py index d251c32..cba4e00 100755 --- a/tests/test_core.py +++ b/tests/test_core.py @@ -2897,7 +2897,7 @@ def test_gstools_variogram(model): def test_pseudo_2d(model): # test data data = np.array([[0.0, 0.0, 1.0], [0.0, 0.0, 3.0], [1.0, 0.0, 6.0]]) - for p_type in [1, 2, 3]: + for p_type in ["pinv", "pinv2", "pinvh"]: # create the krige field krige = model( data[:, 0], @@ -2916,7 +2916,7 @@ def test_pseudo_2d(model): def test_pseudo_3d(model): # test data data = np.array([[0.0, 0.0, 0.0, 1.0], [0.0, 0.0, 0.0, 3.0], [1.0, 0.0, 0.0, 6.0]]) - for p_type in [1, 2, 3]: + for p_type in ["pinv", "pinv2", "pinvh"]: # create the krige field krige = model( data[:, 0], From 492fa9db43c97c663721581a335ea30359231b38 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Tue, 18 Aug 2020 15:34:16 +0200 Subject: [PATCH 32/41] RK: add option for pseudo-inv --- pykrige/rk.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/pykrige/rk.py b/pykrige/rk.py index 2b3b40f..e7b4760 100644 --- a/pykrige/rk.py +++ b/pykrige/rk.py @@ -118,6 +118,8 @@ def __init__( n_closest_points=10, verbose=False, exact_values=True, + pseudo_inv=False, + pseudo_inv_type="pinv", variogram_parameters=None, variogram_function=None, anisotropy_scaling=(1.0, 1.0), @@ -137,6 +139,8 @@ def __init__( self.weight = weight self.verbose = verbose self.exact_values = exact_values + self.pseudo_inv = pseudo_inv + self.pseudo_inv_type = pseudo_inv_type self.anisotropy_scaling = anisotropy_scaling self.anisotropy_angle = anisotropy_angle self.enable_statistics = enable_statistics @@ -170,6 +174,8 @@ def fit(self, x, y, *args, **kwargs): weight=self.weight, verbose=self.verbose, exact_values=self.exact_values, + pseudo_inv=self.pseudo_inv, + pseudo_inv_type=self.pseudo_inv_type, ) add_setup = dict( anisotropy_scaling=self.anisotropy_scaling[0], @@ -319,6 +325,8 @@ def __init__( weight=False, verbose=False, exact_values=True, + pseudo_inv=False, + pseudo_inv_type="pinv", variogram_parameters=None, variogram_function=None, anisotropy_scaling=(1.0, 1.0), @@ -341,6 +349,8 @@ def __init__( n_closest_points=n_closest_points, verbose=verbose, exact_values=exact_values, + pseudo_inv=pseudo_inv, + pseudo_inv_type=pseudo_inv_type, variogram_parameters=variogram_parameters, variogram_function=variogram_function, anisotropy_scaling=anisotropy_scaling, From 464c203ea46469b73f31ec0e9c0bb7d6b0247e4e Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 19 Aug 2020 13:13:45 +0200 Subject: [PATCH 33/41] update examples --- examples/01_ordinary.py | 48 +++++++++++++++++++ examples/02_universal.py | 41 ++++++++++++++++ ...ols_covmodel.py => 03_gstools_covmodel.py} | 0 ...ige_geometric.py => 04_krige_geometric.py} | 31 +++++------- examples/{kriging_1D.py => 05_kriging_1D.py} | 0 ...le_1D.py => 06_exact_values_example_1D.py} | 2 +- ...riging2d.py => 07_regression_kriging2d.py} | 0 examples/{krige_cv.py => 08_krige_cv.py} | 0 ...ing_meuse.ipynb => 09_kriging_meuse.ipynb} | 0 examples/output.asc | 18 +++++++ 10 files changed, 121 insertions(+), 19 deletions(-) create mode 100755 examples/01_ordinary.py create mode 100755 examples/02_universal.py rename examples/{gstools_covmodel.py => 03_gstools_covmodel.py} (100%) rename examples/{krige_geometric.py => 04_krige_geometric.py} (75%) rename examples/{kriging_1D.py => 05_kriging_1D.py} (100%) rename examples/{exact_values_example_1D.py => 06_exact_values_example_1D.py} (98%) rename examples/{regression_kriging2d.py => 07_regression_kriging2d.py} (100%) rename examples/{krige_cv.py => 08_krige_cv.py} (100%) rename examples/{kriging_meuse.ipynb => 09_kriging_meuse.ipynb} (100%) create mode 100644 examples/output.asc diff --git a/examples/01_ordinary.py b/examples/01_ordinary.py new file mode 100755 index 0000000..f0cd162 --- /dev/null +++ b/examples/01_ordinary.py @@ -0,0 +1,48 @@ +""" +Ordinary Kriging Example +======================== + +First we will create a 2D dataset together with the associated x, y grids. + +""" + +import numpy as np +import pykrige.kriging_tools as kt +from pykrige.ok import OrdinaryKriging +import matplotlib.pyplot as plt + + +data = np.array([[0.3, 1.2, 0.47], + [1.9, 0.6, 0.56], + [1.1, 3.2, 0.74], + [3.3, 4.4, 1.47], + [4.7, 3.8, 1.74]]) + +gridx = np.arange(0.0, 5.5, 0.5) +gridy = np.arange(0.0, 5.5, 0.5) + +############################################################################### +# Create the ordinary kriging object. Required inputs are the X-coordinates of +# the data points, the Y-coordinates of the data points, and the Z-values of the +# data points. If no variogram model is specified, defaults to a linear variogram +# model. If no variogram model parameters are specified, then the code automatically +# calculates the parameters by fitting the variogram model to the binned +# experimental semivariogram. The verbose kwarg controls code talk-back, and +# the enable_plotting kwarg controls the display of the semivariogram. + +OK = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear', + verbose=False, enable_plotting=False) + +############################################################################### +# Creates the kriged grid and the variance grid. Allows for kriging on a rectangular +# grid of points, on a masked rectangular grid of points, or with arbitrary points. +# (See OrdinaryKriging.__doc__ for more information.) + +z, ss = OK.execute('grid', gridx, gridy) + +############################################################################### +# Writes the kriged grid to an ASCII grid file and plot it. + +kt.write_asc_grid(gridx, gridy, z, filename="output.asc") +plt.imshow(z) +plt.show() diff --git a/examples/02_universal.py b/examples/02_universal.py new file mode 100755 index 0000000..ea59232 --- /dev/null +++ b/examples/02_universal.py @@ -0,0 +1,41 @@ +""" +Universal Kriging Example +========================= + +In this example we apply a regional linear trend to the kriging system. + +""" + +from pykrige.uk import UniversalKriging +import numpy as np +import matplotlib.pyplot as plt + + +data = np.array([[0.3, 1.2, 0.47], + [1.9, 0.6, 0.56], + [1.1, 3.2, 0.74], + [3.3, 4.4, 1.47], + [4.7, 3.8, 1.74]]) + +gridx = np.arange(0.0, 5.5, 0.5) +gridy = np.arange(0.0, 5.5, 0.5) + +############################################################################### +# Create the universal kriging object. Required inputs are the X-coordinates of +# the data points, the Y-coordinates of the data points, and the Z-values of the +# data points. Variogram is handled as in the ordinary kriging case. +# drift_terms is a list of the drift terms to include; currently supported terms +# are 'regional_linear', 'point_log', and 'external_Z'. Refer to +# UniversalKriging.__doc__ for more information. + +UK = UniversalKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear', + drift_terms=['regional_linear']) + +############################################################################### +# Creates the kriged grid and the variance grid. Allows for kriging on a rectangular +# grid of points, on a masked rectangular grid of points, or with arbitrary points. +# (See UniversalKriging.__doc__ for more information.) + +z, ss = UK.execute('grid', gridx, gridy) +plt.imshow(z) +plt.show() diff --git a/examples/gstools_covmodel.py b/examples/03_gstools_covmodel.py similarity index 100% rename from examples/gstools_covmodel.py rename to examples/03_gstools_covmodel.py diff --git a/examples/krige_geometric.py b/examples/04_krige_geometric.py similarity index 75% rename from examples/krige_geometric.py rename to examples/04_krige_geometric.py index 41732b1..0756b8b 100644 --- a/examples/krige_geometric.py +++ b/examples/04_krige_geometric.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- """ Geometric example ------------------- +================= A small example script showing the usage of the 'geographic' coordinates type for ordinary kriging on a sphere. @@ -9,6 +9,8 @@ from pykrige.ok import OrdinaryKriging import numpy as np +from matplotlib import pyplot as plt + # Make this example reproducible: np.random.seed(89239413) @@ -47,7 +49,9 @@ # Execute on grid: z2, ss2 = OK.execute("grid", grid_lon, grid_lat) +############################################################################### # Print data at equator (last longitude index will show periodicity): + print("Original data:") print("Longitude:", lon.astype(int)) print("Latitude: ", lat.astype(int)) @@ -60,26 +64,17 @@ print("Value: ", np.array_str(z2[5, :], precision=2)) print("Sigma²: ", np.array_str(ss2[5, :], precision=2)) -# ====================================OUTPUT================================== -# >>> Original data: -# >>> Longitude: [122 166 92 138 86 122 136] -# >>> Latitude: [-46 -36 -25 -73 -25 50 -29] -# >>> z: [ 2.75 3.36 2.24 3.07 3.37 5.25 2.82] -# >>> -# >>> Krige at 60° latitude: -# >>> ====================== -# >>> Longitude: [ 0. 60. 120. 180. 240. 300. 360.] -# >>> Value: [ 5.32 5.14 5.3 5.18 5.35 5.61 5.32] -# >>> Sigma²: [ 2.19 1.31 0.41 1.22 2.1 2.46 2.19] -# >>> -# >>> Ignoring curvature: -# >>> ===================== -# >>> Value: [ 4.55 4.72 5.25 4.82 4.61 4.53 4.48] -# >>> Sigma²: [ 3.77 1.99 0.39 1.84 3.52 5.43 7.5 ] -# +############################################################################### # We can see that the data point at longitude 122, latitude 50 correctly # dominates the kriged results, since it is the closest node in spherical # distance metric, as longitude differences scale with cos(latitude). # When kriging using longitude / latitude linearly, the value for grid points # with longitude values further away as longitude is now incorrectly # weighted equally as latitude. + +fig, (ax1, ax2) = plt.subplots(1, 2) +ax1.imshow(z1, extent=[0, 360, -90, 90]) +ax1.set_title("geo-coordinates") +ax2.imshow(z2, extent=[0, 360, -90, 90]) +ax2.set_title("non geo-coordinates") +plt.show() diff --git a/examples/kriging_1D.py b/examples/05_kriging_1D.py similarity index 100% rename from examples/kriging_1D.py rename to examples/05_kriging_1D.py diff --git a/examples/exact_values_example_1D.py b/examples/06_exact_values_example_1D.py similarity index 98% rename from examples/exact_values_example_1D.py rename to examples/06_exact_values_example_1D.py index 5ceeadc..715ac85 100644 --- a/examples/exact_values_example_1D.py +++ b/examples/06_exact_values_example_1D.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- """ Exact Values -================= +============ PyKrige demonstration and usage as a non-exact interpolator in 1D. diff --git a/examples/regression_kriging2d.py b/examples/07_regression_kriging2d.py similarity index 100% rename from examples/regression_kriging2d.py rename to examples/07_regression_kriging2d.py diff --git a/examples/krige_cv.py b/examples/08_krige_cv.py similarity index 100% rename from examples/krige_cv.py rename to examples/08_krige_cv.py diff --git a/examples/kriging_meuse.ipynb b/examples/09_kriging_meuse.ipynb similarity index 100% rename from examples/kriging_meuse.ipynb rename to examples/09_kriging_meuse.ipynb diff --git a/examples/output.asc b/examples/output.asc new file mode 100644 index 0000000..0f3d43b --- /dev/null +++ b/examples/output.asc @@ -0,0 +1,18 @@ +NCOLS 11 +NROWS 11 +XLLCENTER 0.00 +YLLCENTER 0.00 +DX 0.50 +DY 0.50 +NODATA_VALUE -999.00 +0.90 0.95 1.02 1.11 1.20 1.31 1.42 1.52 1.61 1.67 1.72 +0.83 0.89 0.96 1.05 1.15 1.27 1.40 1.51 1.61 1.69 1.73 +0.77 0.81 0.88 0.97 1.08 1.21 1.34 1.47 1.59 1.70 1.74 +0.70 0.74 0.79 0.88 1.00 1.13 1.27 1.40 1.52 1.64 1.69 +0.64 0.67 0.71 0.81 0.93 1.05 1.18 1.31 1.43 1.53 1.59 +0.58 0.61 0.66 0.75 0.85 0.97 1.09 1.21 1.32 1.42 1.48 +0.53 0.56 0.61 0.69 0.78 0.89 1.01 1.12 1.23 1.32 1.39 +0.49 0.51 0.56 0.62 0.71 0.81 0.92 1.04 1.14 1.23 1.30 +0.46 0.48 0.52 0.56 0.63 0.74 0.85 0.96 1.06 1.15 1.22 +0.45 0.46 0.49 0.53 0.58 0.69 0.80 0.90 1.00 1.08 1.15 +0.44 0.46 0.49 0.52 0.58 0.67 0.76 0.86 0.94 1.02 1.09 \ No newline at end of file From d821cc689881d5c6b479225192611b07f51dfc4a Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 19 Aug 2020 13:16:08 +0200 Subject: [PATCH 34/41] blacken examples --- examples/01_ordinary.py | 26 ++++++++++++++++++-------- examples/02_universal.py | 25 +++++++++++++++++-------- examples/07_regression_kriging2d.py | 22 +--------------------- 3 files changed, 36 insertions(+), 37 deletions(-) diff --git a/examples/01_ordinary.py b/examples/01_ordinary.py index f0cd162..de15d1d 100755 --- a/examples/01_ordinary.py +++ b/examples/01_ordinary.py @@ -12,11 +12,15 @@ import matplotlib.pyplot as plt -data = np.array([[0.3, 1.2, 0.47], - [1.9, 0.6, 0.56], - [1.1, 3.2, 0.74], - [3.3, 4.4, 1.47], - [4.7, 3.8, 1.74]]) +data = np.array( + [ + [0.3, 1.2, 0.47], + [1.9, 0.6, 0.56], + [1.1, 3.2, 0.74], + [3.3, 4.4, 1.47], + [4.7, 3.8, 1.74], + ] +) gridx = np.arange(0.0, 5.5, 0.5) gridy = np.arange(0.0, 5.5, 0.5) @@ -30,15 +34,21 @@ # experimental semivariogram. The verbose kwarg controls code talk-back, and # the enable_plotting kwarg controls the display of the semivariogram. -OK = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear', - verbose=False, enable_plotting=False) +OK = OrdinaryKriging( + data[:, 0], + data[:, 1], + data[:, 2], + variogram_model="linear", + verbose=False, + enable_plotting=False, +) ############################################################################### # Creates the kriged grid and the variance grid. Allows for kriging on a rectangular # grid of points, on a masked rectangular grid of points, or with arbitrary points. # (See OrdinaryKriging.__doc__ for more information.) -z, ss = OK.execute('grid', gridx, gridy) +z, ss = OK.execute("grid", gridx, gridy) ############################################################################### # Writes the kriged grid to an ASCII grid file and plot it. diff --git a/examples/02_universal.py b/examples/02_universal.py index ea59232..3b5ba01 100755 --- a/examples/02_universal.py +++ b/examples/02_universal.py @@ -11,11 +11,15 @@ import matplotlib.pyplot as plt -data = np.array([[0.3, 1.2, 0.47], - [1.9, 0.6, 0.56], - [1.1, 3.2, 0.74], - [3.3, 4.4, 1.47], - [4.7, 3.8, 1.74]]) +data = np.array( + [ + [0.3, 1.2, 0.47], + [1.9, 0.6, 0.56], + [1.1, 3.2, 0.74], + [3.3, 4.4, 1.47], + [4.7, 3.8, 1.74], + ] +) gridx = np.arange(0.0, 5.5, 0.5) gridy = np.arange(0.0, 5.5, 0.5) @@ -28,14 +32,19 @@ # are 'regional_linear', 'point_log', and 'external_Z'. Refer to # UniversalKriging.__doc__ for more information. -UK = UniversalKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear', - drift_terms=['regional_linear']) +UK = UniversalKriging( + data[:, 0], + data[:, 1], + data[:, 2], + variogram_model="linear", + drift_terms=["regional_linear"], +) ############################################################################### # Creates the kriged grid and the variance grid. Allows for kriging on a rectangular # grid of points, on a masked rectangular grid of points, or with arbitrary points. # (See UniversalKriging.__doc__ for more information.) -z, ss = UK.execute('grid', gridx, gridy) +z, ss = UK.execute("grid", gridx, gridy) plt.imshow(z) plt.show() diff --git a/examples/07_regression_kriging2d.py b/examples/07_regression_kriging2d.py index 46d25fd..19b927d 100644 --- a/examples/07_regression_kriging2d.py +++ b/examples/07_regression_kriging2d.py @@ -4,6 +4,7 @@ An example of regression kriging """ + import sys from sklearn.svm import SVR @@ -42,24 +43,3 @@ m_rk.fit(p_train, x_train, target_train) print("Regression Score: ", m_rk.regression_model.score(p_test, target_test)) print("RK score: ", m_rk.score(p_test, x_test, target_test)) - -# ====================================OUTPUT================================== - -# ======================================== -# regression model: -# Finished learning regression model -# Finished kriging residuals -# Regression Score: -0.034053855457 -# RK score: 0.66195576665 -# ======================================== -# regression model: -# Finished learning regression model -# Finished kriging residuals -# Regression Score: 0.699771164651 -# RK score: 0.737574040386 -# ======================================== -# regression model: -# Finished learning regression model -# Finished kriging residuals -# Regression Score: 0.527796839838 -# RK score: 0.604908933617 From 4132431f51afbb1e541a3d1da1c09f9bbba2b049 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 19 Aug 2020 14:24:31 +0200 Subject: [PATCH 35/41] move examples from readme to gallery --- README.rst | 201 +++--------------- examples/{01_ordinary.py => 00_ordinary.py} | 0 examples/{02_universal.py => 01_universal.py} | 0 examples/02_kriging3D.py | 103 +++++++++ 4 files changed, 138 insertions(+), 166 deletions(-) rename examples/{01_ordinary.py => 00_ordinary.py} (100%) rename examples/{02_universal.py => 01_universal.py} (100%) create mode 100755 examples/02_kriging3D.py diff --git a/README.rst b/README.rst index 9812510..f1c3e20 100644 --- a/README.rst +++ b/README.rst @@ -60,176 +60,45 @@ If you use conda, PyKrige can be installed from the `conda-forge` channel with, conda install -c conda-forge pykrige -Ordinary Kriging Example -^^^^^^^^^^^^^^^^^^^^^^^^ +Features +^^^^^^^^ -First we will create a 2D dataset together with the associated x, y grids, - -.. code:: python - - import numpy as np - import pykrige.kriging_tools as kt - from pykrige.ok import OrdinaryKriging - - data = np.array([[0.3, 1.2, 0.47], - [1.9, 0.6, 0.56], - [1.1, 3.2, 0.74], - [3.3, 4.4, 1.47], - [4.7, 3.8, 1.74]]) - - gridx = np.arange(0.0, 5.5, 0.5) - gridy = np.arange(0.0, 5.5, 0.5) - - # Create the ordinary kriging object. Required inputs are the X-coordinates of - # the data points, the Y-coordinates of the data points, and the Z-values of the - # data points. If no variogram model is specified, defaults to a linear variogram - # model. If no variogram model parameters are specified, then the code automatically - # calculates the parameters by fitting the variogram model to the binned - # experimental semivariogram. The verbose kwarg controls code talk-back, and - # the enable_plotting kwarg controls the display of the semivariogram. - OK = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear', - verbose=False, enable_plotting=False) - - # Creates the kriged grid and the variance grid. Allows for kriging on a rectangular - # grid of points, on a masked rectangular grid of points, or with arbitrary points. - # (See OrdinaryKriging.__doc__ for more information.) - z, ss = OK.execute('grid', gridx, gridy) - - # Writes the kriged grid to an ASCII grid file. - kt.write_asc_grid(gridx, gridy, z, filename="output.asc") - - -Universal Kriging Example -^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. code:: python - - from pykrige.uk import UniversalKriging - import numpy as np - - data = np.array([[0.3, 1.2, 0.47], - [1.9, 0.6, 0.56], - [1.1, 3.2, 0.74], - [3.3, 4.4, 1.47], - [4.7, 3.8, 1.74]]) - - gridx = np.arange(0.0, 5.5, 0.5) - gridy = np.arange(0.0, 5.5, 0.5) - - # Create the ordinary kriging object. Required inputs are the X-coordinates of - # the data points, the Y-coordinates of the data points, and the Z-values of the - # data points. Variogram is handled as in the ordinary kriging case. - # drift_terms is a list of the drift terms to include; currently supported terms - # are 'regional_linear', 'point_log', and 'external_Z'. Refer to - # UniversalKriging.__doc__ for more information. - UK = UniversalKriging(data[:, 0], data[:, 1], data[:, 2], variogram_model='linear', - drift_terms=['regional_linear']) - - # Creates the kriged grid and the variance grid. Allows for kriging on a rectangular - # grid of points, on a masked rectangular grid of points, or with arbitrary points. - # (See UniversalKriging.__doc__ for more information.) - z, ss = UK.execute('grid', gridx, gridy) - - -Three-Dimensional Kriging Example -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. code:: python - - from pykrige.ok3d import OrdinaryKriging3D - from pykrige.uk3d import UniversalKriging3D - import numpy as np - - data = np.array([[0.1, 0.1, 0.3, 0.9], - [0.2, 0.1, 0.4, 0.8], - [0.1, 0.3, 0.1, 0.9], - [0.5, 0.4, 0.4, 0.5], - [0.3, 0.3, 0.2, 0.7]]) - - gridx = np.arange(0.0, 0.6, 0.05) - gridy = np.arange(0.0, 0.6, 0.01) - gridz = np.arange(0.0, 0.6, 0.1) - - # Create the 3D ordinary kriging object and solves for the three-dimension kriged - # volume and variance. Refer to OrdinaryKriging3D.__doc__ for more information. - ok3d = OrdinaryKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3], - variogram_model='linear') - k3d, ss3d = ok3d.execute('grid', gridx, gridy, gridz) - - # Create the 3D universal kriging object and solves for the three-dimension kriged - # volume and variance. Refer to UniversalKriging3D.__doc__ for more information. - uk3d = UniversalKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3], - variogram_model='linear', drift_terms=['regional_linear']) - k3d, ss3d = uk3d.execute('grid', gridx, gridy, gridz) - - # To use the generic 'specified' drift term, the user must provide the drift values - # at each data point and at every grid point. The following example is equivalent to - # using a linear drift in all three spatial dimensions. Refer to - # UniversalKriging3D.__doc__ for more information. - zg, yg, xg = np.meshgrid(gridz, gridy, gridx, indexing='ij') - uk3d = UniversalKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3], - variogram_model='linear', drift_terms=['specified'], - specified_drift=[data[:, 0], data[:, 1]]) - k3d, ss3d = uk3d.execute('grid', gridx, gridy, gridz, specified_drift_arrays=[xg, yg, zg]) - - # To use the generic 'functional' drift term, the user must provide a callable - # function that takes only the spatial dimensions as arguments. The following example - # is equivalent to using a linear drift only in the x-direction. Refer to - # UniversalKriging3D.__doc__ for more information. - func = lambda x, y, z: x - uk3d = UniversalKriging3D(data[:, 0], data[:, 1], data[:, 2], data[:, 3], - variogram_model='linear', drift_terms=['functional'], - functional_drift=[func]) - k3d, ss3d = uk3d.execute('grid', gridx, gridy, gridz) - - # Note that the use of the 'specified' and 'functional' generic drift capabilities is - # essentially identical in the two-dimensional universal kriging class (except for a - # difference in the number of spatial coordinates for the passed drift functions). - # See UniversalKriging.__doc__ for more information. - - -GSTools covariance models -^^^^^^^^^^^^^^^^^^^^^^^^^ - -You can also use `GSTools `_ -covariance models as input for the ``variogram_model`` in the -PyKrige routines: - -.. code:: python - - import numpy as np - from gstools import Gaussian - from pykrige.ok import OrdinaryKriging - from matplotlib import pyplot as plt - - # conditioning data - data = np.array([[0.3, 1.2, 0.47], - [1.9, 0.6, 0.56], - [1.1, 3.2, 0.74], - [3.3, 4.4, 1.47], - [4.7, 3.8, 1.74]]) - # grid definition for output field - gridx = np.arange(0.0, 5.5, 0.1) - gridy = np.arange(0.0, 6.5, 0.1) - # a GSTools based covariance model - cov_model = Gaussian(dim=2, len_scale=1, anis=0.2, angles=-0.5, var=0.5, nugget=0.1) - # ordinary kriging with pykrige - OK1 = OrdinaryKriging(data[:, 0], data[:, 1], data[:, 2], cov_model) - z1, ss1 = OK1.execute('grid', gridx, gridy) - plt.imshow(z1, origin="lower") - plt.show() - -Which gives: - -.. image:: https://raw.githubusercontent.com/GeoStat-Framework/GSTools/master/docs/source/pics/20_pykrige.png - :width: 400px - :align: center -Have a look at the `documentation about the Covariance Model of GSTools `_. +Krigging algorithms +------------------- + + +.. autosummary:: + :toctree: ./generated/ + + pykrige.ok.OrdinaryKriging + pykrige.uk.UniversalKriging + pykrige.ok3d.OrdinaryKriging3D + pykrige.uk3d.UniversalKriging3D + pykrige.rk.RegressionKriging + + +Wrappers +-------- + +.. autosummary:: + :toctree: ./generated/ + + pykrige.rk.Krige + + +Tools +----- + +.. autosummary:: + :toctree: ./generated/ + + pykrige.kriging_tools.write_asc_grid + pykrige.kriging_tools.read_asc_grid Kriging Parameters Tuning -^^^^^^^^^^^^^^^^^^^^^^^^^ +------------------------- A scikit-learn compatible API for parameter tuning by cross-validation is exposed in `sklearn.model_selection.GridSearchCV `_. @@ -238,7 +107,7 @@ example for a more practical illustration. Regression Kriging -^^^^^^^^^^^^^^^^^^ +------------------ `Regression kriging `_ can be performed with `pykrige.rk.RegressionKriging `_. diff --git a/examples/01_ordinary.py b/examples/00_ordinary.py similarity index 100% rename from examples/01_ordinary.py rename to examples/00_ordinary.py diff --git a/examples/02_universal.py b/examples/01_universal.py similarity index 100% rename from examples/02_universal.py rename to examples/01_universal.py diff --git a/examples/02_kriging3D.py b/examples/02_kriging3D.py new file mode 100755 index 0000000..590ed2c --- /dev/null +++ b/examples/02_kriging3D.py @@ -0,0 +1,103 @@ +""" +Three-Dimensional Kriging Example +================================= + +""" + +from pykrige.ok3d import OrdinaryKriging3D +from pykrige.uk3d import UniversalKriging3D +import numpy as np +from matplotlib import pyplot as plt + + +data = np.array( + [ + [0.1, 0.1, 0.3, 0.9], + [0.2, 0.1, 0.4, 0.8], + [0.1, 0.3, 0.1, 0.9], + [0.5, 0.4, 0.4, 0.5], + [0.3, 0.3, 0.2, 0.7], + ] +) + +gridx = np.arange(0.0, 0.6, 0.05) +gridy = np.arange(0.0, 0.6, 0.01) +gridz = np.arange(0.0, 0.6, 0.1) + +############################################################################### +# Create the 3D ordinary kriging object and solves for the three-dimension kriged +# volume and variance. Refer to OrdinaryKriging3D.__doc__ for more information. + +ok3d = OrdinaryKriging3D( + data[:, 0], data[:, 1], data[:, 2], data[:, 3], variogram_model="linear" +) +k3d1, ss3d = ok3d.execute("grid", gridx, gridy, gridz) + +############################################################################### +# Create the 3D universal kriging object and solves for the three-dimension kriged +# volume and variance. Refer to UniversalKriging3D.__doc__ for more information. + +uk3d = UniversalKriging3D( + data[:, 0], + data[:, 1], + data[:, 2], + data[:, 3], + variogram_model="linear", + drift_terms=["regional_linear"], +) +k3d2, ss3d = uk3d.execute("grid", gridx, gridy, gridz) + +############################################################################### +# To use the generic 'specified' drift term, the user must provide the drift values +# at each data point and at every grid point. The following example is equivalent to +# using a linear drift in all three spatial dimensions. Refer to +# UniversalKriging3D.__doc__ for more information. + +zg, yg, xg = np.meshgrid(gridz, gridy, gridx, indexing="ij") +uk3d = UniversalKriging3D( + data[:, 0], + data[:, 1], + data[:, 2], + data[:, 3], + variogram_model="linear", + drift_terms=["specified"], + specified_drift=[data[:, 0], data[:, 1], data[:, 2]], +) +k3d3, ss3d = uk3d.execute( + "grid", gridx, gridy, gridz, specified_drift_arrays=[xg, yg, zg] +) + +############################################################################### +# To use the generic 'functional' drift term, the user must provide a callable +# function that takes only the spatial dimensions as arguments. The following example +# is equivalent to using a linear drift only in the x-direction. Refer to +# UniversalKriging3D.__doc__ for more information. + +func = lambda x, y, z: x +uk3d = UniversalKriging3D( + data[:, 0], + data[:, 1], + data[:, 2], + data[:, 3], + variogram_model="linear", + drift_terms=["functional"], + functional_drift=[func], +) +k3d4, ss3d = uk3d.execute("grid", gridx, gridy, gridz) + +############################################################################### +# Note that the use of the 'specified' and 'functional' generic drift capabilities is +# essentially identical in the two-dimensional universal kriging class (except for a +# difference in the number of spatial coordinates for the passed drift functions). +# See UniversalKriging.__doc__ for more information. + +fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2) +ax1.imshow(k3d1[:, :, 0]) +ax1.set_title("ordinary kriging") +ax2.imshow(k3d2[:, :, 0]) +ax2.set_title("regional lin. drift") +ax3.imshow(k3d3[:, :, 0]) +ax3.set_title("specified drift") +ax4.imshow(k3d3[:, :, 0]) +ax4.set_title("functional drift") +plt.show() From 8fc2f54dd286848521977e14416a518127a705b0 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 19 Aug 2020 14:25:09 +0200 Subject: [PATCH 36/41] update changelog for 1.5.1 --- CHANGELOG.md | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index a9930b3..76f7ba1 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,22 @@ Changelog ========= +Version 1.5.1 +------------- +*August 20, 2020* + +**New features** + +* update Regression Kriging class to be compatible with all kriging features (#158) +* added option to enable/disable "exact values" to all kriging routines (#153) +* added option to use the pseudo-inverse in all kriging routines (#151) + +**Changes** + +* removed compat-layer for sklearn (#157) +* updated examples in documentation + + Version 1.5.0 ------------- *April 04, 2020* From 722f822335f34fc1303be7efbd37360ca8efc16e Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 19 Aug 2020 14:49:28 +0200 Subject: [PATCH 37/41] update examples and readme --- README.rst | 33 +++++++++++++-------------------- examples/02_kriging3D.py | 11 ++++++----- examples/04_krige_geometric.py | 4 ++-- 3 files changed, 21 insertions(+), 27 deletions(-) diff --git a/README.rst b/README.rst index f1c3e20..e2be898 100644 --- a/README.rst +++ b/README.rst @@ -26,6 +26,7 @@ PyKrige Kriging Toolkit for Python. + Purpose ^^^^^^^ @@ -39,7 +40,9 @@ point and all grid points. With the 'functional' drift capability, the user may of the spatial coordinates that define the drift(s). The package includes a module that contains functions that should be useful in working with ASCII grid files (`*.asc`). -See the documentation at `http://pykrige.readthedocs.io/ `_ for more details. +See the documentation at `http://pykrige.readthedocs.io/ `_ +for more details and examples. + Installation ^^^^^^^^^^^^ @@ -63,38 +66,27 @@ If you use conda, PyKrige can be installed from the `conda-forge` channel with, Features ^^^^^^^^ - Krigging algorithms ------------------- - -.. autosummary:: - :toctree: ./generated/ - - pykrige.ok.OrdinaryKriging - pykrige.uk.UniversalKriging - pykrige.ok3d.OrdinaryKriging3D - pykrige.uk3d.UniversalKriging3D - pykrige.rk.RegressionKriging +* ``OrdinaryKriging``: 2D ordinary kriging with estimated mean +* ``UniversalKriging``: 2D universal kriging providing drift terms +* ``OrdinaryKriging3D``: 3D ordinary kriging +* ``UniversalKriging3D``: 3D universal kriging +* ``RegressionKriging``: An implementation of Regression-Kriging Wrappers -------- -.. autosummary:: - :toctree: ./generated/ - - pykrige.rk.Krige +* ``rk.Krige``: A scikit-learn wrapper class for Ordinary and Universal Kriging Tools ----- -.. autosummary:: - :toctree: ./generated/ - - pykrige.kriging_tools.write_asc_grid - pykrige.kriging_tools.read_asc_grid +* ``kriging_tools.write_asc_grid``: Writes gridded data to ASCII grid file (*.asc) +* ``kriging_tools.read_asc_grid``: Reads ASCII grid file (*.asc) Kriging Parameters Tuning @@ -117,6 +109,7 @@ the ``OrdinaryKriging`` or the ``UniversalKriging`` class, and performs a correc A demonstration of the regression kriging is provided in the `corresponding example `_. + License ^^^^^^^ diff --git a/examples/02_kriging3D.py b/examples/02_kriging3D.py index 590ed2c..7cf14ea 100755 --- a/examples/02_kriging3D.py +++ b/examples/02_kriging3D.py @@ -91,13 +91,14 @@ # difference in the number of spatial coordinates for the passed drift functions). # See UniversalKriging.__doc__ for more information. -fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2) -ax1.imshow(k3d1[:, :, 0]) +fig, (ax1, ax2, ax3, ax4) = plt.subplots(4) +ax1.imshow(k3d1[:, :, 0], origin="lower") ax1.set_title("ordinary kriging") -ax2.imshow(k3d2[:, :, 0]) +ax2.imshow(k3d2[:, :, 0], origin="lower") ax2.set_title("regional lin. drift") -ax3.imshow(k3d3[:, :, 0]) +ax3.imshow(k3d3[:, :, 0], origin="lower") ax3.set_title("specified drift") -ax4.imshow(k3d3[:, :, 0]) +ax4.imshow(k3d3[:, :, 0], origin="lower") ax4.set_title("functional drift") +plt.tight_layout() plt.show() diff --git a/examples/04_krige_geometric.py b/examples/04_krige_geometric.py index 0756b8b..f641ff4 100644 --- a/examples/04_krige_geometric.py +++ b/examples/04_krige_geometric.py @@ -73,8 +73,8 @@ # weighted equally as latitude. fig, (ax1, ax2) = plt.subplots(1, 2) -ax1.imshow(z1, extent=[0, 360, -90, 90]) +ax1.imshow(z1, extent=[0, 360, -90, 90], origin="lower") ax1.set_title("geo-coordinates") -ax2.imshow(z2, extent=[0, 360, -90, 90]) +ax2.imshow(z2, extent=[0, 360, -90, 90], origin="lower") ax2.set_title("non geo-coordinates") plt.show() From 252b56217859d1210be146c35be58e2b0bd21cda Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 19 Aug 2020 14:53:38 +0200 Subject: [PATCH 38/41] update readme doc links --- README.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.rst b/README.rst index e2be898..c45e4e2 100644 --- a/README.rst +++ b/README.rst @@ -94,7 +94,7 @@ Kriging Parameters Tuning A scikit-learn compatible API for parameter tuning by cross-validation is exposed in `sklearn.model_selection.GridSearchCV `_. -See the `Krige CV `_ +See the `Krige CV `_ example for a more practical illustration. @@ -102,12 +102,12 @@ Regression Kriging ------------------ `Regression kriging `_ can be performed -with `pykrige.rk.RegressionKriging `_. +with `pykrige.rk.RegressionKriging `_. This class takes as parameters a scikit-learn regression model, and details of either either the ``OrdinaryKriging`` or the ``UniversalKriging`` class, and performs a correction steps on the ML regression prediction. A demonstration of the regression kriging is provided in the -`corresponding example `_. +`corresponding example `_. License From 72b9c9863d456ad11183f659a7df3e6fa16b2de7 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 19 Aug 2020 15:07:57 +0200 Subject: [PATCH 39/41] update gitignore --- .gitignore | 1 + examples/output.asc | 18 ------------------ 2 files changed, 1 insertion(+), 18 deletions(-) delete mode 100644 examples/output.asc diff --git a/.gitignore b/.gitignore index 353c781..9a5248a 100644 --- a/.gitignore +++ b/.gitignore @@ -121,3 +121,4 @@ pykrige/_version.py *.vtr examples/meuse_example_data/* *.ipynb_checkpoints/* +examples/output.asc diff --git a/examples/output.asc b/examples/output.asc deleted file mode 100644 index 0f3d43b..0000000 --- a/examples/output.asc +++ /dev/null @@ -1,18 +0,0 @@ -NCOLS 11 -NROWS 11 -XLLCENTER 0.00 -YLLCENTER 0.00 -DX 0.50 -DY 0.50 -NODATA_VALUE -999.00 -0.90 0.95 1.02 1.11 1.20 1.31 1.42 1.52 1.61 1.67 1.72 -0.83 0.89 0.96 1.05 1.15 1.27 1.40 1.51 1.61 1.69 1.73 -0.77 0.81 0.88 0.97 1.08 1.21 1.34 1.47 1.59 1.70 1.74 -0.70 0.74 0.79 0.88 1.00 1.13 1.27 1.40 1.52 1.64 1.69 -0.64 0.67 0.71 0.81 0.93 1.05 1.18 1.31 1.43 1.53 1.59 -0.58 0.61 0.66 0.75 0.85 0.97 1.09 1.21 1.32 1.42 1.48 -0.53 0.56 0.61 0.69 0.78 0.89 1.01 1.12 1.23 1.32 1.39 -0.49 0.51 0.56 0.62 0.71 0.81 0.92 1.04 1.14 1.23 1.30 -0.46 0.48 0.52 0.56 0.63 0.74 0.85 0.96 1.06 1.15 1.22 -0.45 0.46 0.49 0.53 0.58 0.69 0.80 0.90 1.00 1.08 1.15 -0.44 0.46 0.49 0.52 0.58 0.67 0.76 0.86 0.94 1.02 1.09 \ No newline at end of file From a13aa041982d723d95fcbd52c78301113778d072 Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 19 Aug 2020 16:03:13 +0200 Subject: [PATCH 40/41] bugfix in readme for PyPI --- README.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.rst b/README.rst index c45e4e2..584a1b3 100644 --- a/README.rst +++ b/README.rst @@ -85,8 +85,8 @@ Wrappers Tools ----- -* ``kriging_tools.write_asc_grid``: Writes gridded data to ASCII grid file (*.asc) -* ``kriging_tools.read_asc_grid``: Reads ASCII grid file (*.asc) +* ``kriging_tools.write_asc_grid``: Writes gridded data to ASCII grid file (\*.asc) +* ``kriging_tools.read_asc_grid``: Reads ASCII grid file (\*.asc) Kriging Parameters Tuning From 48bc43362218b283bcba51b2c0990b72f7e29fbb Mon Sep 17 00:00:00 2001 From: MuellerSeb Date: Wed, 19 Aug 2020 16:11:00 +0200 Subject: [PATCH 41/41] readme typo --- README.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.rst b/README.rst index 584a1b3..822de88 100644 --- a/README.rst +++ b/README.rst @@ -66,8 +66,8 @@ If you use conda, PyKrige can be installed from the `conda-forge` channel with, Features ^^^^^^^^ -Krigging algorithms -------------------- +Kriging algorithms +------------------ * ``OrdinaryKriging``: 2D ordinary kriging with estimated mean * ``UniversalKriging``: 2D universal kriging providing drift terms