如何根据压力数组插入不包含压力信息的数组
How to interpolate an array which doesn't contain pressure information based on a pressure array
假设我有一个变量x，形状为(Time: 127, bottom_top: 58, south_north: 76, west_east: 96)，其中bottom_top 没有明确包含原始压力值。也就是说，原来的压力水平来自一个新的变量p，形状与(Time: 127, bottom_top: 58, south_north: 76, west_east: 96).
相同
我想做的是在x的bottom_top方向上进行插值，从原始压力p到新的压力水平p_new。例如，
x[0,:,1,1].values =
array([404.13263, 404.13217, 404.13174, 404.11353, 404.11615, 404.1033 ,
404.0999 , 404.0954 , 404.0902 , 404.07092, 404.05087, 404.04208,
404.0345 , 404.01535, 403.98822, 403.95865, 403.92953, 403.90146,
403.88913, 403.88214, 403.87677, 403.8698 , 403.86185, 403.8513 ,
403.83362, 403.8124 , 403.791 , 403.76935, 403.7525 , 403.72525,
403.6679 , 403.61584, 403.52277, 403.44308, 403.36752, 403.2787 ,
403.17352, 403.0464 , 402.89468, 402.69656, 402.52823, 402.40897,
402.3129 , 402.2052 , 402.11267, 402.06073, 402.03674, 401.99417,
401.81857, 401.57474, 401.34845, 401.14847, 400.97156, 400.95645,
400.92035, 400.7942 , 400.16495, 398.33502], dtype=float32)
p[0,:,1,1].values =
array([86983.83 , 86819.19 , 86656.23 , 86487.69 , 86326.62 , 86120.52 ,
85831.95 , 85461.48 , 85050.03 , 84598. , 84103.71 , 83613.43 ,
83075.71 , 82500. , 81923.46 , 81346.94 , 80730.56 , 80030.49 ,
79289.44 , 78507.66 , 77642.914 , 76696.05 , 75707.52 , 74719.38 ,
73690.03 , 72619.51 , 71549.4 , 70438.445 , 69286.41 , 68093.85 ,
66819.19 , 65503.797 , 64148.29 , 62711.4 , 61234.523 , 59716.75 ,
58158.11 , 56557.582 , 54874.125 , 53109.348 , 51303.953 , 49458.043 ,
47530.21 , 45519.61 , 43423.426 , 41243.23 , 38992.438 , 36712.22 ,
34389.49 , 31879.646 , 29135.926 , 26400.299 , 23733.72 , 20932.85 ,
18026.195 , 15077.396 , 11595.89 , 7294.6274], dtype=float32)
p_new =
array([1.00758e+05, 9.00880e+04, 8.03100e+04, 7.13940e+04, 6.32770e+04,
5.59420e+04, 4.93250e+04, 4.33770e+04, 3.80270e+04, 3.32140e+04,
2.89180e+04, 2.50580e+04, 2.16130e+04, 1.85420e+04, 1.58370e+04,
1.34660e+04, 1.14090e+04, 9.64300e+03, 8.15100e+03, 6.90200e+03,
5.86000e+03, 4.98700e+03, 4.25500e+03, 3.63600e+03, 3.11200e+03,
2.66800e+03, 2.29100e+03, 1.97200e+03, 1.69900e+03, 1.46600e+03,
1.26700e+03, 1.09600e+03, 9.49000e+02, 8.23000e+02, 7.15000e+02,
6.21000e+02, 5.40000e+02, 4.70000e+02, 4.10000e+02, 3.58000e+02,
3.12000e+02, 2.73000e+02, 2.40000e+02, 2.10000e+02, 1.85000e+02,
1.63000e+02, 1.43000e+02, 1.26000e+02, 1.11000e+02, 9.80000e+01,
8.70000e+01, 7.70000e+01, 6.70000e+01, 5.90000e+01, 5.20000e+01,
4.60000e+01, 4.00000e+01, 3.50000e+01, 3.10000e+01, 2.70000e+01,
2.40000e+01, 2.10000e+01, 1.80000e+01, 1.60000e+01, 1.40000e+01,
1.20000e+01, 1.00000e+01, 9.00000e+00, 8.00000e+00, 7.00000e+00,
6.00000e+00])
有什么方法可以像
那样直接插值吗
x_new = interpolation(x,p,p_new).
提前致谢。
我需要非常频繁地沿各个轴进行插值。我使用以下作为一般解决方案。有更简单的特定解决方案。例如，使用嵌套的 for 循环沿每个轴迭代。我在 numpy 实用程序文件中有 iter_1d 函数。我不确定这是最好的方法，但对我来说效果很好。
# A general solution
import numpy as np
from scipy.interpolate import interp1d
# I've found this useful but never know whether there's a better way
def iter_1d( shape, axis=-1 ):
""" Iterator for all 1d slices along 'axis' in an array of 'shape'.
It returns a generic selection, not a selection of a specific array
Usage: for s in iterate_1d(arr.shape, 2):
print(arr[s], b[s], c[s])
"""
while axis<0:
axis+=len(shape)
limits = list(shape) # [ 127, 58, 76, 96 ]
current = [0]*len(shape) # [ 0, 0, 0, 0 ]
current[axis] = slice(None) # [ 0, :, 0, 0 ]
target = list(range(len(shape)-1, -1, -1)) # Reversed array of axis pointers
# e.g 4d case: target = [ 3, 2, 1, 0 ]
target.remove(axis) # Remove the 'axis' index from this list
# e,g, case: axis = 1 from above target = [ 3, 2, 0 ]
def incr_current():
""" Increments the 'right hand' index in current.
If this overflows, increnent the one before etc...
Returns True until axis 0 overflows, then returns False
"""
nonlocal current
for ix in target:
current[ix] += 1
if current[ix] < limits[ix]: return True
current[ix] = 0
return False
cont=True
while cont:
yield tuple(current)
cont = incr_current()
# Thos can iterate 1d slices along any axis. The question is for axis 1
def do_interpolate( x_new, x, y, axis=-1 ):
""" Will interpolate all 1D slices of the axis specified."""
result = np.empty_like(x_new)
for sel in iter_1d( x.shape, axis ):
do_interp = interp1d(x[sel], y[sel], fill_value='extrapolate')
result[sel] = do_interp(x_new[sel])
return result
x_new = do_interpolate( p_new, p, x, axis = 1 )
所有插值函数都从 1d 自变量开始工作，因此每组 p 都与不同的 x 集相关联，因此必须为每个切片生成插值对象。
对于问题中的具体要求，简单地用循环迭代可能更容易。
# Specific solution
def interpolate( x, p, p_new):
result = np.empty_like(p_new)
for t in range(127):
for ns in range(76):
for ew in range(96):
sel = np.s_[t, :, ns, ew]
do_interp = interp1d(p[sel], x[sel], fill_value='extrapolate')
result[sel] = do_interp(p_new[sel])
return result
x_new = interpolate( x, p, p_new )
假设我有一个变量x，形状为(Time: 127, bottom_top: 58, south_north: 76, west_east: 96)，其中bottom_top 没有明确包含原始压力值。也就是说，原来的压力水平来自一个新的变量p，形状与(Time: 127, bottom_top: 58, south_north: 76, west_east: 96).
我想做的是在x的bottom_top方向上进行插值，从原始压力p到新的压力水平p_new。例如，
x[0,:,1,1].values =
array([404.13263, 404.13217, 404.13174, 404.11353, 404.11615, 404.1033 ,
404.0999 , 404.0954 , 404.0902 , 404.07092, 404.05087, 404.04208,
404.0345 , 404.01535, 403.98822, 403.95865, 403.92953, 403.90146,
403.88913, 403.88214, 403.87677, 403.8698 , 403.86185, 403.8513 ,
403.83362, 403.8124 , 403.791 , 403.76935, 403.7525 , 403.72525,
403.6679 , 403.61584, 403.52277, 403.44308, 403.36752, 403.2787 ,
403.17352, 403.0464 , 402.89468, 402.69656, 402.52823, 402.40897,
402.3129 , 402.2052 , 402.11267, 402.06073, 402.03674, 401.99417,
401.81857, 401.57474, 401.34845, 401.14847, 400.97156, 400.95645,
400.92035, 400.7942 , 400.16495, 398.33502], dtype=float32)
p[0,:,1,1].values =
array([86983.83 , 86819.19 , 86656.23 , 86487.69 , 86326.62 , 86120.52 ,
85831.95 , 85461.48 , 85050.03 , 84598. , 84103.71 , 83613.43 ,
83075.71 , 82500. , 81923.46 , 81346.94 , 80730.56 , 80030.49 ,
79289.44 , 78507.66 , 77642.914 , 76696.05 , 75707.52 , 74719.38 ,
73690.03 , 72619.51 , 71549.4 , 70438.445 , 69286.41 , 68093.85 ,
66819.19 , 65503.797 , 64148.29 , 62711.4 , 61234.523 , 59716.75 ,
58158.11 , 56557.582 , 54874.125 , 53109.348 , 51303.953 , 49458.043 ,
47530.21 , 45519.61 , 43423.426 , 41243.23 , 38992.438 , 36712.22 ,
34389.49 , 31879.646 , 29135.926 , 26400.299 , 23733.72 , 20932.85 ,
18026.195 , 15077.396 , 11595.89 , 7294.6274], dtype=float32)
p_new =
array([1.00758e+05, 9.00880e+04, 8.03100e+04, 7.13940e+04, 6.32770e+04,
5.59420e+04, 4.93250e+04, 4.33770e+04, 3.80270e+04, 3.32140e+04,
2.89180e+04, 2.50580e+04, 2.16130e+04, 1.85420e+04, 1.58370e+04,
1.34660e+04, 1.14090e+04, 9.64300e+03, 8.15100e+03, 6.90200e+03,
5.86000e+03, 4.98700e+03, 4.25500e+03, 3.63600e+03, 3.11200e+03,
2.66800e+03, 2.29100e+03, 1.97200e+03, 1.69900e+03, 1.46600e+03,
1.26700e+03, 1.09600e+03, 9.49000e+02, 8.23000e+02, 7.15000e+02,
6.21000e+02, 5.40000e+02, 4.70000e+02, 4.10000e+02, 3.58000e+02,
3.12000e+02, 2.73000e+02, 2.40000e+02, 2.10000e+02, 1.85000e+02,
1.63000e+02, 1.43000e+02, 1.26000e+02, 1.11000e+02, 9.80000e+01,
8.70000e+01, 7.70000e+01, 6.70000e+01, 5.90000e+01, 5.20000e+01,
4.60000e+01, 4.00000e+01, 3.50000e+01, 3.10000e+01, 2.70000e+01,
2.40000e+01, 2.10000e+01, 1.80000e+01, 1.60000e+01, 1.40000e+01,
1.20000e+01, 1.00000e+01, 9.00000e+00, 8.00000e+00, 7.00000e+00,
6.00000e+00])
有什么方法可以像
那样直接插值吗x_new = interpolation(x,p,p_new).
提前致谢。
我需要非常频繁地沿各个轴进行插值。我使用以下作为一般解决方案。有更简单的特定解决方案。例如，使用嵌套的 for 循环沿每个轴迭代。我在 numpy 实用程序文件中有 iter_1d 函数。我不确定这是最好的方法，但对我来说效果很好。
# A general solution
import numpy as np
from scipy.interpolate import interp1d
# I've found this useful but never know whether there's a better way
def iter_1d( shape, axis=-1 ):
""" Iterator for all 1d slices along 'axis' in an array of 'shape'.
It returns a generic selection, not a selection of a specific array
Usage: for s in iterate_1d(arr.shape, 2):
print(arr[s], b[s], c[s])
"""
while axis<0:
axis+=len(shape)
limits = list(shape) # [ 127, 58, 76, 96 ]
current = [0]*len(shape) # [ 0, 0, 0, 0 ]
current[axis] = slice(None) # [ 0, :, 0, 0 ]
target = list(range(len(shape)-1, -1, -1)) # Reversed array of axis pointers
# e.g 4d case: target = [ 3, 2, 1, 0 ]
target.remove(axis) # Remove the 'axis' index from this list
# e,g, case: axis = 1 from above target = [ 3, 2, 0 ]
def incr_current():
""" Increments the 'right hand' index in current.
If this overflows, increnent the one before etc...
Returns True until axis 0 overflows, then returns False
"""
nonlocal current
for ix in target:
current[ix] += 1
if current[ix] < limits[ix]: return True
current[ix] = 0
return False
cont=True
while cont:
yield tuple(current)
cont = incr_current()
# Thos can iterate 1d slices along any axis. The question is for axis 1
def do_interpolate( x_new, x, y, axis=-1 ):
""" Will interpolate all 1D slices of the axis specified."""
result = np.empty_like(x_new)
for sel in iter_1d( x.shape, axis ):
do_interp = interp1d(x[sel], y[sel], fill_value='extrapolate')
result[sel] = do_interp(x_new[sel])
return result
x_new = do_interpolate( p_new, p, x, axis = 1 )
所有插值函数都从 1d 自变量开始工作，因此每组 p 都与不同的 x 集相关联，因此必须为每个切片生成插值对象。
对于问题中的具体要求，简单地用循环迭代可能更容易。
# Specific solution
def interpolate( x, p, p_new):
result = np.empty_like(p_new)
for t in range(127):
for ns in range(76):
for ew in range(96):
sel = np.s_[t, :, ns, ew]
do_interp = interp1d(p[sel], x[sel], fill_value='extrapolate')
result[sel] = do_interp(p_new[sel])
return result
x_new = interpolate( x, p, p_new )