CODEWITHSUNDEEP
pythonpython-2.7numpy
user2201041
user2201041

Reputation:

How to interpolate into a rotated grid?

I have

[[0, 1, 2, 3, 4],
 [1, 2, 3, 4, 5],
 [2, 3, 4, 5, 6],
 [3, 4, 5, 6, 7],
 [4, 5, 6, 7, 8]]

and I want to interpolate it into a rotated grid that has a corner on the left edge. Similar to:

[[2, ~2, 2],
 [~4, ~4, ~4],
 [6, ~6, 6]]

(I use ~ to denote approximate values. )

(Of course, my actual data is more complex. The scenario is that I want to map DEM data by pixel onto a rotated image.)

Here is the setup:

import numpy
from scipy import interpolate as interp

grid = numpy.ndarray((5, 5))
for I in range(grid.shape[0]):
    for j in range(grid.shape[1]):
        grid[I, j] = I + j

grid = ndimage.interpolation.shift(
    ndimage.interpolation.rotate(grid, -45, reshape=False),
    -1)

source_x, source_y = numpy.meshgrid(
    numpy.arange(0, 5), numpy.arange(0, 5))
target_x, target_y = numpy.meshgrid(
    numpy.arange(0, 2), numpy.arange(0, 2))

print(interp.griddata(
    numpy.array([source_x.ravel(), source_y.ravel()]).T,
    grid.ravel(),
    target_x, target_y))

This is giving me:

[[2.4467   2.6868 2.4467]
 [4.       4.     4.    ]
 [5.5553   5.3132 5.5553]]

This is promising. However, the rotation and shift values are hard-coded, and I should at least be able to get the upper-left corner exact.

I do know the indices of the corners of the grid I wish to interpolate to. That is, I have

upper_left = 2, 0
upper_right = 0, 2
lower_right = 4, 2
lower_left = 2, 4

Upvotes: 2

Views: 1452

Answers (2)

Paul Panzer
Paul Panzer

Reputation: 53029

This may not be builtin enough for your taste, but here is a method that uses your starting point (grid corners) more directly, and applies spline interpolation (cubic per default).

import numpy as np
from scipy.interpolate import RectBivariateSpline

# input data
data = np.array([[0, 1, 2, 3, 4],
                 [1, 2, 3, 4, 5],
                 [2, 3, 4, 5, 6],
                 [3, 4, 5, 6, 7],
                 [4, 5, 6, 7, 8]])

upper_left = 2, 0
upper_right = 0, 2
lower_right = 2, 4   # note that I swapped this
lower_left = 4, 2    # and this
n_steps = 3, 3


# build interpolator
m, n = data.shape
x, y = np.arange(m), np.arange(n)

interpolator = RectBivariateSpline(x, y, data)

# build grid
ul,ur,ll,lr = map(np.array, (upper_left,upper_right,lower_left,lower_right))
assert np.allclose(ul + lr, ur + ll)    # make sure edges are parallel

x, y = ul[:, None, None] \
       + np.outer(ll-ul, np.linspace(0.0, 1.0, n_steps[0]))[:, :, None] \
       + np.outer(ur-ul, np.linspace(0.0, 1.0, n_steps[1]))[:, None, :]

# intepolate on grid
print(interpolator.ev(x, y))

Prints:

[[2. 2. 2.]
 [4. 4. 4.]
 [6. 6. 6.]]

Upvotes: 1

user2201041
user2201041

Reputation:

While this does answer my question, I hope there is a more built-in way of doing things. I am still seeking a better approach.

This is really two problems: rotating the original grid, and then interpolating. Rotating the grid and then translating it to the correct upper-left corner can be done with an Affine Transformation.

skimage provides an easy function for this purpose.

from skimage.transform import AffineTransform, warp
# Set up grid as in question
transform = AffineTransform(rotation=-math.pi / 4,
                            scale=(math.sqrt(2)/2, math.sqrt(2)/2),
                            translations=(0,2))
grid = warp(grid, transform)

The result of this is

[[2. 2. 2. 2. 2.]
 [3. 3. 3. 3. 3.]
 [4. 4. 4. 4. 4.]
 [5. 5. 5. 5. 5.]
 [6. 6. 6. 6. 6.]]

which may then be simply resampled as desired.

In general, if we have a grid with dimensions x and y, and coordinates p1, p2, p3, p4 (starting from upper-left, going clockwise) we want to rotate to, we have

rotation = math.atan2(p4.x - p1.x, p4.y - p1.y)
scale = (math.sqrt((p2.y - p1.y) ** 2 + (p2.x - p1.x) ** 2) / x,
         math.sqrt((p4.y - p1.y) ** 2 + (p4.x - p1.x) ** 2) / y)
translation = (0, p1.y)

Upvotes: 0

Related Questions