numpy 数组上的 Scipy 插值 [英] Scipy interpolation on a numpy array

问题描述

我有一个按以下方式定义的查找表:

I have a lookup table that is defined the following way:

       | <1    2    3    4    5+
-------|----------------------------
<10000 | 3.6   6.5  9.1  11.5 13.8
20000  | 3.9   7.3  10.0 13.1 15.9
20000+ | 4.5   9.2  12.2 14.8 18.2


TR_ua1 = np.array([ [3.6, 6.5, 9.1, 11.5, 13.8],
                    [3.9, 7.3, 10.0, 13.1, 15.9],
                    [4.5, 9.2, 12.2, 14.8, 18.2] ])

• 标题行元素是 (hh) <1,2,3,4,5+
• 标题列 (inc) 元素为 <10000、20000、20001+

用户将输入一个值，例如 (1.3, 25,000)、(0.2, 50,000) 等等.scipy.interpolate() 应该进行插值以确定正确的值.

The user will input a value example (1.3, 25,000), (0.2, 50,000), so on. scipy.interpolate() should interpolate to determine the correct value.

目前，我能做到这一点的唯一方法是使用一堆 if/elifs，如下所示.我很确定有更好、更有效的方法来做到这一点

Currently, the only way I can do this is with a bunch of if/elifs as exemplified below. I'm pretty sure there is a better, more efficient way of doing this

这是我到目前为止所得到的:

Here's what I've got so far:

import numpy as np
from scipy import interpolate

if (ua == 1):
    if (inc <= low_inc):  # low_inc = 10,000
      if (hh <= 1):
        return TR_ua1[0][0]
      elif (hh >= 1 & hh < 2):
        return interpolate( (1, 2), (TR_ua1[0][1], TR_ua1[0][2]) )

推荐答案

更新内容以反映您在上面的说明.你的问题现在清楚多了，谢谢！

Updated things to reflect your clarifications above. Your question is much clearer now, thanks!

基本上，您只想在任意点插入二维数组.

Basically, you're just wanting to interpolate a 2D array at an arbitrary point.

scipy.ndimage.map_coordinates 就是你想要的....

scipy.ndimage.map_coordinates is what you want....

据我所知，您有一个z"值的二维数组，范围从某个 xmin 到 xmax，以及每个方向上的 ymin 到 ymax.

As I understand your question, you have a 2D array of "z" values that ranges from some xmin to xmax, and ymin to ymax in each direction.

要从数组的边缘返回值的那些边界坐标之外的任何内容.

Anything outside of those bounding coordinates you want to return values from the edges of the array.

map_coordinates 有几个选项来处理网格边界外的点，但它们都不是你想要的.相反，我们可以将边界外的任何东西设置为位于边缘，并像往常一样使用 map_coordinates.

map_coordinates has several options to handle points outside the boundaries of the grid, but none of them do exactly what you want. Instead, we can just set anything outside the boundaries to lie on the edge, and use map_coordinates as usual.

因此，要使用 map_coordinates，您需要转换实际坐标:

So, to use map_coordinates, you need to turn your actual coodinates:

       | <1    2    3    4    5+
-------|----------------------------
<10000 | 3.6   6.5  9.1  11.5 13.8
20000  | 3.9   7.3  10.0 13.1 15.9
20000+ | 4.5   9.2  12.2 14.8 18.2

进入索引坐标:

       |  0     1    2    3    4
-------|----------------------------
   0   | 3.6   6.5  9.1  11.5 13.8
   1   | 3.9   7.3  10.0 13.1 15.9
   2   | 4.5   9.2  12.2 14.8 18.2

请注意，您的边界在每个方向上的表现都不同……在 x 方向上，事情表现得很顺利，但在 y 方向上，您要求硬"中断，其中 y=20000 -->3.9，但 y=20000.000001 --> 4.5.

Note that your boundaries behave differently in each direction... In the x-direction, things behave smoothly, but in the y-direction, you're asking for a "hard" break, where y=20000 --> 3.9, but y=20000.000001 --> 4.5.

举个例子:

import numpy as np
from scipy.ndimage import map_coordinates

#-- Setup ---------------------------
z = np.array([ [3.6, 6.5, 9.1, 11.5, 13.8],
               [3.9, 7.3, 10.0, 13.1, 15.9],
               [4.5, 9.2, 12.2, 14.8, 18.2] ])
ny, nx = z.shape
xmin, xmax = 1, 5
ymin, ymax = 10000, 20000

# Points we want to interpolate at
x1, y1 = 1.3, 25000
x2, y2 = 0.2, 50000
x3, y3 = 2.5, 15000

# 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

# To deal with the "hard" break, we'll have to treat y differently, 
# so we're ust setting the min here...
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)

# Deal with the "hard" break in the y-direction
yi = (ny - 2) * (yi - ymin) / (ymax - ymin)
yi[yi > 1] = 2.0

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

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

这产生:

1.3 , 25000 --> 5.1807375
0.2 , 50000 --> 4.5
2.5 , 15000 --> 8.12252371652

希望有帮助...

这篇关于numpy 数组上的 Scipy 插值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！