在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!

基本上，您只是想在任意点插入2D数组.

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

scipy.ndimage.map_coordinates 是您想要的....

据我所知，您有一个二维的"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，您需要打开实际的coodinates:

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

希望这对您有帮助...

Hopefully that helps...

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