Интерполяция Python / Scipy (map_coordinates) - PullRequest
Новая
       6

Интерполяция Python / Scipy (map_coordinates)

11 голосов
/ 26 февраля 2011

Я пытаюсь сделать некоторую интерполяцию со Сципи. Я рассмотрел много примеров, но не могу найти именно то, что хочу.

Допустим, у меня есть некоторые данные, где переменные строки и столбца могут варьироваться от 0 до 1. Дельта-изменения между каждой строкой и столбцом не всегда одинаковы (см. Ниже). 

      | 0.00   0.25  0.80  1.00
------|----------------------------
0.00  | 1.40   6.50  1.50  1.80
0.60  | 8.90   7.30  1.10  1.09
1.00  | 4.50   9.20  1.80  1.20

Теперь я хочу иметь возможность взять набор точек x, y и определить интерполированные значения. Я знаю, что могу сделать это с помощью map_coordinates. Мне интересно, есть ли какой-нибудь простой / умный способ сделать значение x, y соответствующим индексом в массиве данных.

Например, если я введу x, y = 0,60, 0,25, то я должен вернуть правильный индекс для интерполяции. В этом случае это будет 1,0, 1,0 с 0,60, 0,25 будет отображаться точно во второй строке и втором столбце. х = 0,3 соответствует 0,5, так как он находится на полпути между 0,00 и 0,60.

Я знаю, как получить желаемый результат, но я уверен, что есть очень быстрый / понятный одно или двухстрочный (или уже существующая функция), который может сделать это, чтобы сделать мой код более понятным. По сути, он должен кусочно интерполировать между некоторыми массивами.

Вот пример (в значительной степени основанный на коде из Интерполяция Сципи на массиве NumPy ) - я поместил TODO, куда пошла бы эта новая функция: 

from scipy.ndimage import map_coordinates
from numpy import arange
import numpy as np
#            0.000,  0.175,  0.817,  1.000
z = array([ [ 3.6,    6.5,    9.1,    11.5],    # 0.0000
            [ 3.9,   -7.3,    10.0,   13.1],    # 0.2620
            [ 1.9,    8.3,   -15.0,  -12.1],    # 0.6121
            [-4.5,    9.2,    12.2,   14.8] ])  # 1.0000
ny, nx = z.shape
xmin, xmax = 0., 1.
ymin, ymax = 0., 1.

xrange = array([0.000,  0.175,  0.817,  1.000 ])
yrange = array([0.0000, 0.2620, 0.6121, 1.0000])

# Points we want to interpolate at
x1, y1 = 0.20, 0.45
x2, y2 = 0.30, 0.85
x3, y3 = 0.95, 1.00

# To make our lives easier down the road, let's
# turn these into arrays of x & y coords
xi = np.array([x1, x2, x3], dtype=np.float)
yi = np.array([y1, y2, y3], dtype=np.float)

# Now, we'll set points outside the boundaries to lie along an edge
xi[xi > xmax] = xmax
xi[xi < xmin] = xmin

yi[yi > ymax] = ymax
yi[yi < ymin] = ymin

# We need to convert these to (float) indicies
#   (xi should range from 0 to (nx - 1), etc)
xi = (nx - 1) * (xi - xmin) / (xmax - xmin)
yi = (ny - 1) * (yi - ymin) / (ymax - ymin)
# TODO: Instead, xi and yi need to be mapped as described.  This can only work with
# even spacing...something like:
#xi = SomeInterpFunction(xi, xrange)
#yi = SomeInterpFunction(yi, yrange)

# Now we actually interpolate
# map_coordinates does cubic interpolation by default,
# use "order=1" to preform bilinear interpolation instead...
print xi
print yi
z1, z2, z3 = map_coordinates(z, [yi, xi], order=1)

# Display the results
for X, Y, Z in zip((x1, x2, x3), (y1, y2, y3), (z1, z2, z3)):
    print X, ',', Y, '-->', Z

1 Ответ

16 голосов
/ 26 февраля 2011

Я думаю, что вы хотите двумерный сплайн на прямоугольной структурированной сетке : 

import numpy
from scipy import interpolate
x = numpy.array([0.0, 0.60, 1.0])
y = numpy.array([0.0, 0.25, 0.80, 1.0])
z = numpy.array([ 
   [ 1.4 ,  6.5 ,  1.5 ,  1.8 ],
   [ 8.9 ,  7.3 ,  1.1 ,  1.09],
   [ 4.5 ,  9.2 ,  1.8 ,  1.2 ]])
# you have to set kx and ky small for this small example dataset
# 3 is more usual and is the default
# s=0 will ensure this interpolates.  s>0 will smooth the data
# you can also specify a bounding box outside the data limits
# if you want to extrapolate
sp = interpolate.RectBivariateSpline(x, y, z, kx=2, ky=2, s=0)

sp([0.60], [0.25])  # array([[ 7.3]])
sp([0.25], [0.60])  # array([[ 2.66427408]])
...