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)
我的目标是将邻近气象站的气象数据插入到具有精确坐标的点中。在 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)