Skip to content

Commit

Permalink
ND refactoring of adjust_for_anisotropy (#33)
Browse files Browse the repository at this point in the history
  • Loading branch information
@rth
rth committed Jan 27, 2017
1 parent 82d7245 commit 9b1af38
Show file tree
Hide file tree
Showing 6 changed files with 139 additions and 133 deletions.
118 changes: 54 additions & 64 deletions pykrige/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,11 +11,8 @@
scipy (scipy.optimize.minimize())
Functions:
adjust_for_anisotropy(x, y, xcenter, ycenter, scaling, angle):
Returns X and Y arrays of adjusted data coordinates. Angle is CCW.
adjust_for_anisotropy_3d(x, y, z, xcenter, ycenter, zcenter, scaling_y,
scaling_z, angle_x, angle_y, angle_z):
Returns X, Y, Z arrays of adjusted data coordinates. Angles are CCW about
_adjust_for_anisotropy(X, y, center, scaling, angle):
Returns X_adj array of adjusted data coordinates. Angles are CCW about
specified axes. Scaling is applied in rotated coordinate system.
initialize_variogram_model(x, y, z, variogram_model, variogram_model_parameters,
variogram_function, nlags):
Expand Down Expand Up @@ -127,66 +124,59 @@ def euclid3_to_great_circle(euclid3_distance):
return 180.0 - 360.0/np.pi*np.arccos(0.5*euclid3_distance)


def adjust_for_anisotropy(x, y, xcenter, ycenter, scaling, angle):
def _adjust_for_anisotropy(X, center, scaling, angle):
"""Adjusts data coordinates to take into account anisotropy.
Can also be used to take into account data scaling."""

x -= xcenter
y -= ycenter
xshape = x.shape
yshape = y.shape
x = x.flatten()
y = y.flatten()

coords = np.vstack((x, y))
stretch = np.array([[1, 0], [0, scaling]])
rotate = np.array([[np.cos(-angle * np.pi/180.0), -np.sin(-angle * np.pi/180.0)],
[np.sin(-angle * np.pi/180.0), np.cos(-angle * np.pi/180.0)]])
rotated_coords = np.dot(stretch, np.dot(rotate, coords))
x = rotated_coords[0, :].reshape(xshape)
y = rotated_coords[1, :].reshape(yshape)
x += xcenter
y += ycenter

return x, y


def adjust_for_anisotropy_3d(x, y, z, xcenter, ycenter, zcenter, scaling_y,
scaling_z, angle_x, angle_y, angle_z):
"""Adjusts data coordinates to take into account anisotropy.
Can also be used to take into account data scaling."""

x -= xcenter
y -= ycenter
z -= zcenter
xshape = x.shape
yshape = y.shape
zshape = z.shape
x = x.flatten()
y = y.flatten()
z = z.flatten()

coords = np.vstack((x, y, z))
stretch = np.array([[1., 0., 0.], [0., scaling_y, 0.], [0., 0., scaling_z]])
rotate_x = np.array([[1., 0., 0.],
[0., np.cos(-angle_x * np.pi/180.), -np.sin(-angle_x * np.pi/180.)],
[0., np.sin(-angle_x * np.pi/180.), np.cos(-angle_x * np.pi/180.)]])
rotate_y = np.array([[np.cos(-angle_y * np.pi/180.), 0., np.sin(-angle_y * np.pi/180.)],
[0., 1., 0.],
[-np.sin(-angle_y * np.pi/180.), 0., np.cos(-angle_y * np.pi/180.)]])
rotate_z = np.array([[np.cos(-angle_z * np.pi/180.), -np.sin(-angle_z * np.pi/180.), 0.],
[np.sin(-angle_z * np.pi/180.), np.cos(-angle_z * np.pi/180.), 0.],
[0., 0., 1.]])
rot_tot = np.dot(rotate_z, np.dot(rotate_y, rotate_x))
rotated_coords = np.dot(stretch, np.dot(rot_tot, coords))
x = rotated_coords[0, :].reshape(xshape)
y = rotated_coords[1, :].reshape(yshape)
z = rotated_coords[2, :].reshape(zshape)
x += xcenter
y += ycenter
z += zcenter

return x, y, z
Can also be used to take into account data scaling.
Parameters
----------
X: ndarray
float array [n_samples, n_dim], the input array of coordinates
center: ndarray
float array [n_dim], the coordinate of centers
scaling: ndarray
float array [n_dim - 1], the scaling of last two dimensions
angle : ndarray
float array [2*n_dim - 3], the anysotropy angle (degrees)
Returns
-------
X_adj : ndarray
float array [n_samples, n_dim], the X array adjusted for anisotropy.
"""

center = np.asarray(center)[None, :]
angle = np.asarray(angle)*np.pi/180

X -= center

Ndim = X.shape[1]

if Ndim == 1:
raise NotImplementedError('Not implemnented yet?')
elif Ndim == 2:
stretch = np.array([[1, 0], [0, scaling[0]]])
rot_tot = np.array([[np.cos(-angle[0]), -np.sin(-angle[0])],
[np.sin(-angle[0]), np.cos(-angle[0])]])
elif Ndim == 3:
stretch = np.array([[1., 0., 0.], [0., scaling[0], 0.], [0., 0., scaling[1]]])
rotate_x = np.array([[1., 0., 0.],
[0., np.cos(-angle[0]), -np.sin(-angle[0])],
[0., np.sin(-angle[0]), np.cos(-angle[0])]])
rotate_y = np.array([[np.cos(-angle[1]), 0., np.sin(-angle[1])],
[0., 1., 0.],
[-np.sin(-angle[1]), 0., np.cos(-angle[1])]])
rotate_z = np.array([[np.cos(-angle[2]), -np.sin(-angle[2]), 0.],
[np.sin(-angle[2]), np.cos(-angle[2]), 0.],
[0., 0., 1.]])
rot_tot = np.dot(rotate_z, np.dot(rotate_y, rotate_x))
else:
raise ValueError("Adjust for anysotropy function doesn't support ND spaces where N>3")
X_adj = np.dot(stretch, np.dot(rot_tot, X.T)).T

X_adj += center

return X_adj


def initialize_variogram_model(x, y, z, variogram_model, variogram_model_parameters,
Expand Down
23 changes: 13 additions & 10 deletions pykrige/ok.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@
import matplotlib.pyplot as plt
from . import variogram_models
from . import core
from .core import _adjust_for_anisotropy


class OrdinaryKriging:
Expand Down Expand Up @@ -220,9 +221,10 @@ def __init__(self, x, y, z, variogram_model='linear', variogram_parameters=None,
if self.verbose:
print("Adjusting data for anisotropy...")
self.X_ADJUSTED, self.Y_ADJUSTED = \
core.adjust_for_anisotropy(np.copy(self.X_ORIG), np.copy(self.Y_ORIG),
self.XCENTER, self.YCENTER,
self.anisotropy_scaling, self.anisotropy_angle)
_adjust_for_anisotropy(np.vstack((self.X_ORIG, self.Y_ORIG)).T,
[self.XCENTER, self.YCENTER],
[self.anisotropy_scaling],
[self.anisotropy_angle]).T
elif coordinates_type == 'geographic':
# Leave everything as is in geographic case.
# May be open to discussion?
Expand Down Expand Up @@ -307,11 +309,10 @@ def update_variogram_model(self, variogram_model, variogram_parameters=None,
self.anisotropy_scaling = anisotropy_scaling
self.anisotropy_angle = anisotropy_angle
self.X_ADJUSTED, self.Y_ADJUSTED = \
core.adjust_for_anisotropy(np.copy(self.X_ORIG),
np.copy(self.Y_ORIG),
self.XCENTER, self.YCENTER,
self.anisotropy_scaling,
self.anisotropy_angle)
_adjust_for_anisotropy(np.vstack((self.X_ORIG, self.Y_ORIG)).T,
[self.XCENTER, self.YCENTER],
[self.anisotropy_scaling],
[self.anisotropy_angle]).T
elif self.coordinates_type == 'geographic':
self.anisotropy_scaling = 1.0
self.anisotropy_angle = 0.0
Expand Down Expand Up @@ -619,8 +620,10 @@ def execute(self, style, xpoints, ypoints, mask=None, backend='vectorized', n_cl
raise ValueError("style argument must be 'grid', 'points', or 'masked'")

if self.coordinates_type == 'euclidean':
xpts, ypts = core.adjust_for_anisotropy(xpts, ypts, self.XCENTER, self.YCENTER,
self.anisotropy_scaling, self.anisotropy_angle)
xpts, ypts = _adjust_for_anisotropy(np.vstack((xpts, ypts)).T,
[self.XCENTER, self.YCENTER],
[self.anisotropy_scaling],
[self.anisotropy_angle]).T
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':
Expand Down
28 changes: 16 additions & 12 deletions pykrige/ok3d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
import matplotlib.pyplot as plt
from . import variogram_models
from . import core
from .core import _adjust_for_anisotropy


class OrdinaryKriging3D:
Expand Down Expand Up @@ -230,10 +231,11 @@ def __init__(self, x, y, z, val, variogram_model='linear', variogram_parameters=
if self.verbose:
print("Adjusting data for anisotropy...")
self.X_ADJUSTED, self.Y_ADJUSTED, self.Z_ADJUSTED = \
core.adjust_for_anisotropy_3d(np.copy(self.X_ORIG), np.copy(self.Y_ORIG), np.copy(self.Z_ORIG),
self.XCENTER, self.YCENTER, self.ZCENTER, self.anisotropy_scaling_y,
self.anisotropy_scaling_z, self.anisotropy_angle_x, self.anisotropy_angle_y,
self.anisotropy_angle_z)
_adjust_for_anisotropy(np.vstack((self.X_ORIG, self.Y_ORIG, self.Z_ORIG)).T,
[self.XCENTER, self.YCENTER, self.ZCENTER],
[self.anisotropy_scaling_y, self.anisotropy_scaling_z],
[self.anisotropy_angle_x,
self.anisotropy_angle_y, self.anisotropy_angle_z]).T

self.variogram_model = variogram_model
if self.variogram_model not in self.variogram_dict.keys() and self.variogram_model != 'custom':
Expand Down Expand Up @@ -301,10 +303,11 @@ def update_variogram_model(self, variogram_model, variogram_parameters=None, var
self.anisotropy_angle_y = anisotropy_angle_y
self.anisotropy_angle_z = anisotropy_angle_z
self.X_ADJUSTED, self.Y_ADJUSTED, self.Z_ADJUSTED = \
core.adjust_for_anisotropy_3d(np.copy(self.X_ORIG), np.copy(self.Y_ORIG), np.copy(self.Z_ORIG),
self.XCENTER, self.YCENTER, self.ZCENTER, self.anisotropy_scaling_y,
self.anisotropy_scaling_z, self.anisotropy_angle_x,
self.anisotropy_angle_y, self.anisotropy_angle_z)
_adjust_for_anisotropy(np.vstack((self.X_ORIG, self.Y_ORIG, self.Z_ORIG)).T,
[self.XCENTER, self.YCENTER, self.ZCENTER],
[self.anisotropy_scaling_y, self.anisotropy_scaling_z],
[self.anisotropy_angle_x,
self.anisotropy_angle_y, self.anisotropy_angle_z]).T

self.variogram_model = variogram_model
if self.variogram_model not in self.variogram_dict.keys() and self.variogram_model != 'custom':
Expand Down Expand Up @@ -615,10 +618,11 @@ def execute(self, style, xpoints, ypoints, zpoints, mask=None, backend='vectoriz
else:
raise ValueError("style argument must be 'grid', 'points', or 'masked'")

xpts, ypts, zpts = core.adjust_for_anisotropy_3d(xpts, ypts, zpts, self.XCENTER, self.YCENTER, self.ZCENTER,
self.anisotropy_scaling_y, self.anisotropy_scaling_z,
self.anisotropy_angle_x, self.anisotropy_angle_y,
self.anisotropy_angle_z)
xpts, ypts, zpts = _adjust_for_anisotropy(np.vstack((xpts, ypts, zpts)).T,
[self.XCENTER, self.YCENTER, self.ZCENTER],
[self.anisotropy_scaling_y, self.anisotropy_scaling_z],
[self.anisotropy_angle_x, self.anisotropy_angle_y,
self.anisotropy_angle_z]).T

if style != 'masked':
mask = np.zeros(npt, dtype='bool')
Expand Down
42 changes: 22 additions & 20 deletions pykrige/test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,29 +57,31 @@ def setUp(self):

def test_core_adjust_for_anisotropy(self):

x = np.array([1.0, 0.0, -1.0, 0.0])
y = np.array([0.0, 1.0, 0.0, -1.0])
rotated_x, rotated_y = core.adjust_for_anisotropy(x, y, 0.0, 0.0, 2.0, 90.0)
self.assertTrue(np.allclose(rotated_x, np.array([0.0, 1.0, 0.0, -1.0])))
self.assertTrue(np.allclose(rotated_y, np.array([-2.0, 0.0, 2.0, 0.0])))
X = np.array([[1.0, 0.0, -1.0, 0.0],
[0.0, 1.0, 0.0, -1.0]]).T
X_adj = core._adjust_for_anisotropy(X, [0.0, 0.0], [2.0], [90.0])
self.assertTrue(np.allclose(X_adj[:, 0], np.array([0.0, 1.0, 0.0, -1.0])))
self.assertTrue(np.allclose(X_adj[:, 1], np.array([-2.0, 0.0, 2.0, 0.0])))

def test_core_adjust_for_anisotropy_3d(self):

x = np.array([1.0, 0.0, 0.0])
y = np.array([0.0, 1.0, 0.0])
z = np.array([0.0, 0.0, 1.0])
rotated_x, rotated_y, rotated_z = core.adjust_for_anisotropy_3d(x, y, z, 0., 0., 0., 2., 2., 90., 0., 0.)
self.assertTrue(np.allclose(rotated_x, np.array([1., 0., 0.])))
self.assertTrue(np.allclose(rotated_y, np.array([0., 0., 2.])))
self.assertTrue(np.allclose(rotated_z, np.array([0., -2., 0.])))
rotated_x, rotated_y, rotated_z = core.adjust_for_anisotropy_3d(x, y, z, 0., 0., 0., 2., 2., 0., 90., 0.)
self.assertTrue(np.allclose(rotated_x, np.array([0., 0., -1.])))
self.assertTrue(np.allclose(rotated_y, np.array([0., 2., 0.])))
self.assertTrue(np.allclose(rotated_z, np.array([2., 0., 0.])))
rotated_x, rotated_y, rotated_z = core.adjust_for_anisotropy_3d(x, y, z, 0., 0., 0., 2., 2., 0., 0., 90.)
self.assertTrue(np.allclose(rotated_x, np.array([0., 1., 0.])))
self.assertTrue(np.allclose(rotated_y, np.array([-2., 0., 0.])))
self.assertTrue(np.allclose(rotated_z, np.array([0., 0., 2.])))
# this is a bad examples, as the X matrix is symmetric
# and insensitive to transpositions
X = np.array([[1.0, 0.0, 0.0],
[0.0, 1.0, 0.0],
[0.0, 0.0, 1.0]]).T
X_adj = core._adjust_for_anisotropy(X, [0., 0., 0.], [2., 2.], [90., 0., 0.])
self.assertTrue(np.allclose(X_adj[:, 0], np.array([1., 0., 0.])))
self.assertTrue(np.allclose(X_adj[:, 1], np.array([0., 0., 2.])))
self.assertTrue(np.allclose(X_adj[:, 2], np.array([0., -2., 0.])))
X_adj = core._adjust_for_anisotropy(X, [0., 0., 0.], [2., 2.], [0., 90., 0.])
self.assertTrue(np.allclose(X_adj[:, 0], np.array([0., 0., -1.])))
self.assertTrue(np.allclose(X_adj[:, 1], np.array([0., 2., 0.])))
self.assertTrue(np.allclose(X_adj[:, 2], np.array([2., 0., 0.])))
X_adj = core._adjust_for_anisotropy(X, [0., 0., 0.], [2., 2.], [0., 0., 90.])
self.assertTrue(np.allclose(X_adj[:, 0], np.array([0., 1., 0.])))
self.assertTrue(np.allclose(X_adj[:, 1], np.array([-2., 0., 0.])))
self.assertTrue(np.allclose(X_adj[:, 2], np.array([0., 0., 2.])))

def test_core_initialize_variogram_model(self):

Expand Down
31 changes: 17 additions & 14 deletions pykrige/uk.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@
import matplotlib.pyplot as plt
from . import variogram_models
from . import core
from .core import _adjust_for_anisotropy
import warnings


Expand Down Expand Up @@ -265,9 +266,10 @@ def __init__(self, x, y, z, variogram_model='linear', variogram_parameters=None,
if self.verbose:
print("Adjusting data for anisotropy...")
self.X_ADJUSTED, self.Y_ADJUSTED = \
core.adjust_for_anisotropy(np.copy(self.X_ORIG), np.copy(self.Y_ORIG),
self.XCENTER, self.YCENTER,
self.anisotropy_scaling, self.anisotropy_angle)
_adjust_for_anisotropy(np.vstack((self.X_ORIG, self.Y_ORIG)).T,
[self.XCENTER, self.YCENTER],
[self.anisotropy_scaling],
[self.anisotropy_angle]).T

self.variogram_model = variogram_model
if self.variogram_model not in self.variogram_dict.keys() and self.variogram_model != 'custom':
Expand Down Expand Up @@ -366,10 +368,10 @@ def __init__(self, x, y, z, variogram_model='linear', variogram_parameters=None,
point_log = np.atleast_2d(np.squeeze(np.array(point_drift, copy=True)))
self.point_log_array = np.zeros(point_log.shape)
self.point_log_array[:, 2] = point_log[:, 2]
x_adj, y_adj = core.adjust_for_anisotropy(point_log[:, 0], point_log[:, 1], self.XCENTER, self.YCENTER,
self.anisotropy_scaling, self.anisotropy_angle)
self.point_log_array[:, 0] = x_adj
self.point_log_array[:, 1] = y_adj
self.point_log_array[:, :2] = _adjust_for_anisotropy(np.vstack((point_log[:, 0], point_log[:, 1])).T,
[self.XCENTER, self.YCENTER],
[self.anisotropy_scaling],
[self.anisotropy_angle])
if self.verbose:
print("Implementing external point-logarithmic drift; number of points =",
self.point_log_array.shape[0], '\n')
Expand Down Expand Up @@ -512,11 +514,10 @@ def update_variogram_model(self, variogram_model, variogram_parameters=None,
self.anisotropy_scaling = anisotropy_scaling
self.anisotropy_angle = anisotropy_angle
self.X_ADJUSTED, self.Y_ADJUSTED = \
core.adjust_for_anisotropy(np.copy(self.X_ORIG),
np.copy(self.Y_ORIG),
self.XCENTER, self.YCENTER,
self.anisotropy_scaling,
self.anisotropy_angle)
_adjust_for_anisotropy(np.vstack((self.X_ORIG, self.Y_ORIG)).T,
[self.XCENTER, self.YCENTER],
[self.anisotropy_scaling],
[self.anisotropy_angle]).T

self.variogram_model = variogram_model
if self.variogram_model not in self.variogram_dict.keys() and self.variogram_model != 'custom':
Expand Down Expand Up @@ -938,8 +939,10 @@ def execute(self, style, xpoints, ypoints, mask=None, backend='vectorized', spec
"instantiation of UniversalKriging class.", RuntimeWarning)

xy_points_original = np.concatenate((xpts[:, np.newaxis], ypts[:, np.newaxis]), axis=1)
xpts, ypts = core.adjust_for_anisotropy(xpts, ypts, self.XCENTER, self.YCENTER,
self.anisotropy_scaling, self.anisotropy_angle)
xpts, ypts = _adjust_for_anisotropy(np.vstack((xpts, ypts)).T,
[self.XCENTER, self.YCENTER],
[self.anisotropy_scaling],
[self.anisotropy_angle]).T
xy_points = np.concatenate((xpts[:, np.newaxis], ypts[:, np.newaxis]), axis=1)
xy_data = np.concatenate((self.X_ADJUSTED[:, np.newaxis], self.Y_ADJUSTED[:, np.newaxis]), axis=1)

Expand Down

0 comments on commit 9b1af38

Please sign in to comment.