Tuesday, 9 October 2018

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



from How to interpolate into a rotated grid?

No comments:

Post a Comment