Python地理空间插值(气象数据)

Python geospatial interpolation (meteorological data)

我的目标是将邻近气象站的气象数据插入到具有精确坐标的点中。在 SciPy 文档中,我找到了有关多维插值的信息 (from scipy.interpolate import griddata)。但老实说,我不明白如何让这段代码适合我的任务。

我有 df 以站点坐标和大气压力值作为输入(以及没有站点的地方的坐标)

有人可以帮我解决这个问题吗?

(https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html#id4) - 多元数据插值 (griddata)¶

据我所知,SciPy 的 griddata 非常适合在网格上插入多个值(顾名思义),但如果您只有一点你想获得一个值,我建议使用 SciPy 的 interpolate.interp2d 函数 (see docs)。在任何一种情况下,您仍然会设置一个 2D 网格来导出插值函数,但我认为使用后一种方法更容易获得您想要的结果。

这是一个简单的示例,它使用了我想象中您拥有的那种数据。

from scipy.interpolate import interp2d
import numpy as np
import pandas as pd

# Set up your dataframe with station data (you might have lon and lat instead of x and y)
df = pd.DataFrame({'x':[1,2,1,2], 'y':[1,1,2,2], 'Pa':[10,10,20,20]})

# Create your linearly-spaced 2D grid, the higher num_pts the higher the resolution 
num_pts = 10
x_grid = np.linspace(min(df['x']), max(df['x']), num_pts)
y_grid = np.linspace(min(df['y']), max(df['y']), num_pts)

# Set up your interpolation function (I would recommend linear interpolation 
# for your problem, but you can change this if you want with the 'kind' parameter)
f = interp2d(df['x'], df['y'], df['Pa'], kind='linear')

# Let's say we want to interpolate "Pa" values at:
x_prime, y_prime = (1.5, 1.5)

# Just plug the new coords into your function
pa_new = f(x_prime, y_prime)
pa_new 

>>> array([15.])

让您了解 griddata 可以为您做什么:

from scipy.interpolate import griddata

X,Y = np.meshgrid(x_grid, y_grid)

Z = griddata([(x,y) for x,y in zip(df['x'],df['y'])], df['pa'], (X, Y), method='linear')

plt.subplot(111)
plt.scatter(df['x'], df['y'], 100, 'r', edgecolor='w')
plt.imshow(Z, extent=(min(df['x']-0.1),max(df['x']+0.1),min(df['y']-0.1),max(df['y']+0.1)), origin='lower')
plt.colorbar()

Z
>>> array([[10.        , 10.        , 10.        , 10.        , 10.        ,
        10.        , 10.        , 10.        , 10.        , 10.        ],
       [11.11111111, 11.11111111, 11.11111111, 11.11111111, 11.11111111,
        11.11111111, 11.11111111, 11.11111111, 11.11111111, 11.11111111],
       [12.22222222, 12.22222222, 12.22222222, 12.22222222, 12.22222222,
        12.22222222, 12.22222222, 12.22222222, 12.22222222, 12.22222222],
       [13.33333333, 13.33333333, 13.33333333, 13.33333333, 13.33333333,
        13.33333333, 13.33333333, 13.33333333, 13.33333333, 13.33333333],
       [14.44444444, 14.44444444, 14.44444444, 14.44444444, 14.44444444,
        14.44444444, 14.44444444, 14.44444444, 14.44444444, 14.44444444],
       [15.55555556, 15.55555556, 15.55555556, 15.55555556, 15.55555556,
        15.55555556, 15.55555556, 15.55555556, 15.55555556, 15.55555556],
       [16.66666667, 16.66666667, 16.66666667, 16.66666667, 16.66666667,
        16.66666667, 16.66666667, 16.66666667, 16.66666667, 16.66666667],
       [17.77777778, 17.77777778, 17.77777778, 17.77777778, 17.77777778,
        17.77777778, 17.77777778, 17.77777778, 17.77777778, 17.77777778],
       [18.88888889, 18.88888889, 18.88888889, 18.88888889, 18.88888889,
        18.88888889, 18.88888889, 18.88888889, 18.88888889, 18.88888889],
       [20.        , 20.        , 20.        , 20.        , 20.        ,
        20.        , 20.        , 20.        , 20.        , 20.        ]])

在此处找到解决方案:

Inverse Distance Weighted (IDW) Interpolation with Python

IDW 插值对我来说已经足够了,但是@user6386471,感谢您的贡献!

def linear_rbf(x, y, z, xi, yi):
dist = distance_matrix(x,y, xi,yi)

# Mutual pariwise distances between observations
internal_dist = distance_matrix(x,y, x,y)

# Now solve for the weights such that mistfit at the observations is minimized
weights = np.linalg.solve(internal_dist, z)

# Multiply the weights for each interpolated point by the distances
zi =  np.dot(dist.T, weights)
return zi

(这里使用distance_matrix函数:)

def distance_matrix(x0, y0, x1, y1):
obs = np.vstack((x0, y0)).T
interp = np.vstack((x1, y1)).T

# Make a distance matrix between pairwise observations
# Note: from <
# (Yay for ufuncs!)
d0 = np.subtract.outer(obs[:,0], interp[:,0])
d1 = np.subtract.outer(obs[:,1], interp[:,1])

return np.hypot(d0, d1)