Skip to content

Commit

Permalink
Adapted the _adjust_for_anisotropy code
Browse files Browse the repository at this point in the history
  • Loading branch information
@rth
rth committed Jan 8, 2017
1 parent 4ba0af5 commit be51a38
Show file tree
Hide file tree
Showing 5 changed files with 43 additions and 93 deletions.
62 changes: 0 additions & 62 deletions pykrige/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,68 +182,6 @@ def _adjust_for_anisotropy(X, center, scaling, angle):
return X_adj


def adjust_for_anisotropy(x, y, xcenter, ycenter, 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


def initialize_variogram_model(x, y, z, variogram_model, variogram_model_parameters,
variogram_function, nlags, weight, coordinates_type):
"""Initializes the variogram model for kriging according
Expand Down
8 changes: 5 additions & 3 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
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
8 changes: 5 additions & 3 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
30 changes: 17 additions & 13 deletions pykrige/uk3d.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
import warnings


Expand Down Expand Up @@ -265,10 +266,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 @@ -377,10 +379,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 @@ -778,10 +781,11 @@ def execute(self, style, xpoints, ypoints, zpoints, mask=None, backend='vectoriz
warnings.warn("Provided specified drift values, but 'specified' drift was not initialized during "
"instantiation of UniversalKriging3D class.", RuntimeWarning)

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 All @@ -806,4 +810,4 @@ def execute(self, style, xpoints, ypoints, zpoints, mask=None, backend='vectoriz
kvalues = kvalues.reshape((nz, ny, nx))
sigmasq = sigmasq.reshape((nz, ny, nx))

return kvalues, sigmasq
return kvalues, sigmasq

0 comments on commit be51a38

Please sign in to comment.