FIX Geographic coordinates (#116)

Author: mjziebarth

pykrige/ok.py

@@ -459,9 +459,14 @@ def print_statistics(self):
     def _get_kriging_matrix(self, n):
         """Assembles the kriging matrix."""
 
-        xy = np.concatenate((self.X_ADJUSTED[:, np.newaxis],
-                             self.Y_ADJUSTED[:, np.newaxis]), axis=1)
-        d = cdist(xy, xy, 'euclidean')
+        if self.coordinates_type == 'euclidean':
+            xy = np.concatenate((self.X_ADJUSTED[:, np.newaxis],
+                                 self.Y_ADJUSTED[:, np.newaxis]), axis=1)
+            d = cdist(xy, xy, 'euclidean')
+        elif self.coordinates_type == 'geographic':
+            d = core.great_circle_distance(self.X_ADJUSTED[:,np.newaxis],
+                                           self.Y_ADJUSTED[:,np.newaxis],
+                                           self.X_ADJUSTED, self.Y_ADJUSTED)
         a = np.zeros((n+1, n+1))
         a[:n, :n] = - self.variogram_function(self.variogram_model_parameters,
                                               d)
@@ -660,6 +665,11 @@ def execute(self, style, xpoints, ypoints, mask=None, backend='vectorized',
             raise ValueError("style argument must be 'grid', "
                              "'points', or 'masked'")
 
+        if n_closest_points is not None and n_closest_points <= 1:
+            # If this is not checked, nondescriptive errors emerge
+            # later in the code.
+            raise ValueError("n_closest_points has to be at least two!")
+
         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]
@@ -699,28 +709,16 @@ def execute(self, style, xpoints, ypoints, mask=None, backend='vectorized',
                                                 [self.XCENTER, self.YCENTER],
                                                 [self.anisotropy_scaling],
                                                 [self.anisotropy_angle]).T
+            # Prepare for cdist:
             xy_data = np.concatenate((self.X_ADJUSTED[:, np.newaxis],
                                       self.Y_ADJUSTED[:, np.newaxis]), axis=1)
             xy_points = np.concatenate((xpts[:, np.newaxis],
                                         ypts[:, np.newaxis]), axis=1)
         elif self.coordinates_type == 'geographic':
-            # Quick version: Only difference between euclidean and spherical
-            # space regarding kriging is the distance metric.
-            # Since the relationship between three dimensional euclidean
-            # distance and great circle distance on the sphere is monotonous,
-            # use the existing (euclidean) infrastructure for nearest neighbour
-            # search and distance calculation in euclidean three space and
-            # convert distances to great circle distances afterwards.
-            lon_d = self.X_ADJUSTED[:, np.newaxis] * np.pi / 180.0
-            lat_d = self.Y_ADJUSTED[:, np.newaxis] * np.pi / 180.0
-            xy_data = np.concatenate((np.cos(lon_d) * np.cos(lat_d),
-                                      np.sin(lon_d) * np.cos(lat_d),
-                                      np.sin(lat_d)), axis=1)
-            lon_p = xpts[:, np.newaxis] * np.pi / 180.0
-            lat_p = ypts[:, np.newaxis] * np.pi / 180.0
-            xy_points = np.concatenate((np.cos(lon_p) * np.cos(lat_p),
-                                        np.sin(lon_p) * np.cos(lat_p),
-                                        np.sin(lat_p)), axis=1)
+            # In spherical coordinates, we do not correct for anisotropy.
+            # Also, we don't use scipy.spatial.cdist, so we do not have to
+            # format the input data accordingly.
+            pass
 
         if style != 'masked':
             mask = np.zeros(npt, dtype='bool')
@@ -741,13 +739,36 @@ def execute(self, style, xpoints, ypoints, mask=None, backend='vectorized',
             c_pars = {key: getattr(self, key) for key in ['Z', 'eps', 'variogram_model_parameters', 'variogram_function']}
 
         if n_closest_points is not None:
+            if self.coordinates_type == 'geographic':
+                # To make use of the KDTree, we have to convert the
+                # spherical coordinates into three dimensional Euclidean
+                # coordinates, since the standard KDTree cannot handle
+                # the periodicity.
+                # Do the conversion just for the step involving the KDTree:
+                lon_d = self.X_ADJUSTED[:, np.newaxis] * np.pi / 180.0
+                lat_d = self.Y_ADJUSTED[:, np.newaxis] * np.pi / 180.0
+                xy_data = np.concatenate((np.cos(lon_d) * np.cos(lat_d),
+                                          np.sin(lon_d) * np.cos(lat_d),
+                                          np.sin(lat_d)), axis=1)
+                lon_p = xpts[:, np.newaxis] * np.pi / 180.0
+                lat_p = ypts[:, np.newaxis] * np.pi / 180.0
+                xy_points = np.concatenate((np.cos(lon_p) * np.cos(lat_p),
+                                            np.sin(lon_p) * np.cos(lat_p),
+                                            np.sin(lat_p)), axis=1)
+
             from scipy.spatial import cKDTree
             tree = cKDTree(xy_data)
             bd, bd_idx = tree.query(xy_points, k=n_closest_points, eps=0.0)
+
             if self.coordinates_type == 'geographic':
-                # Convert euclidean distances to great circle distances:
-                bd = core.euclid3_to_great_circle(bd)
-            
+                # Between the nearest neighbours from Euclidean search,
+                # calculate the great circle distance using the standard method:
+                x_points = np.tile(xpts[:,np.newaxis],(1,n_closest_points))
+                y_points = np.tile(ypts[:,np.newaxis],(1,n_closest_points))
+                bd = core.great_circle_distance(x_points, y_points,
+                                                self.X_ADJUSTED[bd_idx],
+                                                self.Y_ADJUSTED[bd_idx])
+
             if backend == 'loop':
                 zvalues, sigmasq = \
                     self._exec_loop_moving_window(a, bd, mask, bd_idx)
@@ -760,11 +781,13 @@ def execute(self, style, xpoints, ypoints, mask=None, backend='vectorized',
                 raise ValueError('Specified backend {} for a moving window '
                                  'is not supported.'.format(backend))
         else:
-            bd = cdist(xy_points,  xy_data, 'euclidean')
-            if self.coordinates_type == 'geographic':
-                # Convert euclidean distances to great circle distances:
-                bd = core.euclid3_to_great_circle(bd)
-            
+            if self.coordinates_type == 'euclidean':
+                bd = cdist(xy_points,  xy_data, 'euclidean')
+            elif self.coordinates_type == 'geographic':
+                bd = core.great_circle_distance(xpts[:,np.newaxis],
+                                                ypts[:,np.newaxis],
+                                                self.X_ADJUSTED, self.Y_ADJUSTED)
+
             if backend == 'vectorized':
                 zvalues, sigmasq = self._exec_vector(a, bd, mask)
             elif backend == 'loop':

pykrige/tests/test_core.py

@@ -14,6 +14,7 @@
 import pytest
 from pytest import approx
 from numpy.testing import assert_allclose
+from scipy.spatial.distance import cdist
 
 from pykrige import kriging_tools as kt
 from pykrige import core
@@ -1996,7 +1997,7 @@ def test_geometric_code():
             rtol=1e-5)
 
 
-def test_ok_geometric():
+def test_ok_geographic():
     # Generate random data:
     np.random.seed(89239413)
     lon = 360.0*np.random.rand(50, 1)
@@ -2013,3 +2014,87 @@ def test_ok_geometric():
 
     # Execute on grid:
     z, ss = OK.execute('grid', grid_lon, grid_lat)
+
+
+def test_ok_geographic_vs_euclid():
+    # Generate some random data close to the north pole.
+    # Then we use a polar projected 2d euclidean coordinate
+    # system and compare the kriging results in that coordinate
+    # system with the geographic-option results.
+    # If data point distance to the north pole is small enough
+    # (choose maximum 0.01 degrees), the differences due to curvature
+    # should be negligible.
+    np.random.seed(89239413)
+    from_north = 1e-2*np.random.random(5)
+    lat = 90.0-from_north
+    lon = 360.0*np.random.random(5)
+    z = np.random.random(5)
+    z -= z.mean()
+    x = from_north * np.cos(np.deg2rad(lon))
+    y = from_north * np.sin(np.deg2rad(lon))
+
+    # Generate grids:
+    grid_lon = 360.0 * np.linspace(0, 1, 50)
+    grid_from_north = np.linspace(0,0.01,10)
+    grid_lat = 90.0 - grid_from_north
+    grid_x = grid_from_north[:,np.newaxis] * np.cos(np.deg2rad(grid_lon[np.newaxis,:]))
+    grid_y = grid_from_north[:,np.newaxis] * np.sin(np.deg2rad(grid_lon[np.newaxis,:]))
+    grid_lon, grid_lat = np.meshgrid(grid_lon, grid_lat, indexing='xy')
+
+    # Flatten the grids:
+    grid_x = grid_x.flatten()
+    grid_y = grid_y.flatten()
+    grid_lon = grid_lon.flatten()
+    grid_lat = grid_lat.flatten()
+
+    # Calculate and compare distance matrices ensuring that that part
+    # of the workflow works as intended (tested: 2e-9 is currently the
+    # match for this setup):
+    d_eucl = cdist(np.concatenate([x[:,np.newaxis],y[:,np.newaxis]],axis=1),
+                   np.concatenate([grid_x[:,np.newaxis],grid_y[:,np.newaxis]],axis=1))
+    d_geo = core.great_circle_distance(lon[:,np.newaxis], lat[:,np.newaxis],
+                 grid_lon[np.newaxis,:], grid_lat[np.newaxis,:])
+    assert_allclose(d_eucl,d_geo, rtol=2e-9)
+
+    # Create ordinary kriging objects:
+    OK_geo = OrdinaryKriging(lon, lat, z, variogram_model='linear', verbose=False,
+                             enable_plotting=False, coordinates_type='geographic')
+    OK_xy = OrdinaryKriging(x, y, z, variogram_model='linear', verbose=False,
+                            enable_plotting=False)
+    OK_wrong = OrdinaryKriging(lon, lat, z, variogram_model='linear', verbose=False,
+                               enable_plotting=False)
+
+    # Execute on grid:
+    zgeo, ss = OK_geo.execute('points', grid_lon, grid_lat)
+    zxy, ss = OK_xy.execute('points', grid_x, grid_y)
+    zwrong, ss = OK_wrong.execute('points', grid_lon, grid_lat)
+
+    # Assert equivalence / difference (tested: 2e-5 is currently the
+    # match for this setup):
+    assert_allclose(zgeo, zxy, rtol=2e-5)
+    assert not np.any(zgeo == 0)
+    assert np.abs((zgeo-zwrong)/zgeo).max() > 1.0
+
+
+def test_ok_geometric_closest_points():
+    # Generate random data:
+    np.random.seed(89239413)
+    lon = 360.0*np.random.rand(50, 1)
+    lat = 180.0*np.random.rand(50, 1) - 90.0
+    z = np.random.rand(50, 1)
+
+    # Generate grid:
+    grid_lon = 360.0*np.random.rand(120, 1)
+    grid_lat = 180.0*np.random.rand(120, 1) - 90.0
+
+    # Create ordinary kriging object:
+    OK = OrdinaryKriging(lon, lat, z, variogram_model='linear', verbose=False,
+                         enable_plotting=False, coordinates_type='geographic')
+
+    # Execute on grid:
+    with pytest.raises(ValueError):
+        # Test OK raising ValueError when closest_points == 1:
+        z, ss = OK.execute('grid', grid_lon, grid_lat, n_closest_points=1,
+                           backend='C')
+    z, ss = OK.execute('grid', grid_lon, grid_lat, n_closest_points=5,
+                       backend='C')

release_notes.md

@@ -1,5 +1,15 @@
 # PyKrige Changelog
 
+## Version 1.4.1
+*??? Date *
+
+### New features
+* Added method to obtain variogram model points. PR[#94](https://github.com/bsmurphy/PyKrige/pull/94) by [Daniel Mejía Raigosa](https://github.com/Daniel-M)
+
+### Bug fixes
+* Fixed OrdinaryKriging readme example. PR[#107](https://github.com/bsmurphy/PyKrige/pull/107) by [Harry Matchette-Downes](https://github.com/harrymd)
+* Fixed kriging matrix not being calculated correctly for geographic coordinates. PR[99](https://github.com/bsmurphy/PyKrige/pull/99) by [Mike Rilee](https://github.com/michaelleerilee)
+
 ## Version 1.4
 *April 24, 2018*
 

setup.py

@@ -14,7 +14,7 @@
 ext_errors = (CCompilerError, DistutilsExecError, DistutilsPlatformError)
 
 NAME = 'PyKrige'
-VERSION = '1.4.0'
+VERSION = '1.4.1'
 AUTHOR = 'Benjamin S. Murphy'
 EMAIL = '[email protected]'
 URL = 'https://github.com/bsmurphy/PyKrige'