考虑到底图将被逐步淘汰,我正在从底图转向 cartopy。我之前使用过 basemap.interp 功能来插入数据,例如假设我有 1 度分辨率 (180x360) 的数据,我将运行以下命令来插值到 0.5 度。
import numpy as np
from mpl_toolkits import basemap
Old_Lon = np.linspace(-180,180,360)
Old_Lat = np.linspace(-90,90,180)
New_Lon = np.linspace(-180,180,720)
New_Lat = np.linspace(-90,90,360)
New_Lon,New_Lat = np.meshgrid(New_Lon,New_Lat)
New_Data = basemap.interp(Old_Data,Old_Lon,Old_Lat,New_Lon,New_Lat,order=0)
order
为我提供了从最近邻、双线性等中进行选择的选项。是否有其他方法可以以简单的方式执行此操作?我见过 scipy 有插值,但我不知道如何应用它。任何帮助将不胜感激!
我最终决定从 Basemap 中获取原始代码并将其变成一个独立的函数 - 我将向 cartopy 人员推荐它来实现它,因为它是一个有用的功能。在这里发帖可能对其他人有用:
def Interp(datain,xin,yin,xout,yout,interpolation='NearestNeighbour'):
"""
Interpolates a 2D array onto a new grid (only works for linear grids),
with the Lat/Lon inputs of the old and new grid. Can perfom nearest
neighbour interpolation or bilinear interpolation (of order 1)'
This is an extract from the basemap module (truncated)
"""
# Mesh Coordinates so that they are both 2D arrays
xout,yout = np.meshgrid(xout,yout)
# compute grid coordinates of output grid.
delx = xin[1:]-xin[0:-1]
dely = yin[1:]-yin[0:-1]
xcoords = (len(xin)-1)*(xout-xin[0])/(xin[-1]-xin[0])
ycoords = (len(yin)-1)*(yout-yin[0])/(yin[-1]-yin[0])
xcoords = np.clip(xcoords,0,len(xin)-1)
ycoords = np.clip(ycoords,0,len(yin)-1)
# Interpolate to output grid using nearest neighbour
if interpolation == 'NearestNeighbour':
xcoordsi = np.around(xcoords).astype(np.int32)
ycoordsi = np.around(ycoords).astype(np.int32)
dataout = datain[ycoordsi,xcoordsi]
# Interpolate to output grid using bilinear interpolation.
elif interpolation == 'Bilinear':
xi = xcoords.astype(np.int32)
yi = ycoords.astype(np.int32)
xip1 = xi+1
yip1 = yi+1
xip1 = np.clip(xip1,0,len(xin)-1)
yip1 = np.clip(yip1,0,len(yin)-1)
delx = xcoords-xi.astype(np.float32)
dely = ycoords-yi.astype(np.float32)
dataout = (1.-delx)*(1.-dely)*datain[yi,xi] + \
delx*dely*datain[yip1,xip1] + \
(1.-delx)*dely*datain[yip1,xi] + \
delx*(1.-dely)*datain[yi,xip1]
return dataout
--
SciPy 插值例程返回一个函数,您可以调用该函数来执行插值。对于规则网格上的最近邻插值,您可以使用
scipy.interpolate.RegularGridInterpolator
:
import numpy as np
from scipy.interpolate import RegularGridInterpolator
nearest_function = RegularGridInterpolator(
(old_lon, old_lat), old_data, method="nearest", bounds_error=False
)
new_data = np.array(
[[nearest_function([i, j]) for j in new_lat] for i in new_lon]
).squeeze()
但这并不完美,因为
lon=175
都是填充值。 (如果我没有设置 bounds_error=False
那么你会得到一个错误。)在这种情况下,你需要询问你想如何环绕日期线。一个简单的解决方案是将 lon=0
行复制到数组的末尾,并将其称为 lon=180
。
如果有一天您想要线性或高阶插值(如果您的数据是点而不是单元格,我建议您这样做),您可以使用
scipy.interpolate.RectBivariateSpline
:
import numpy as np
from scipy.interpolate import RectBivariateSpline
old_step = 10
old_lon = np.arange(-180, 180, old_step)
old_lat = np.arange(-90, 90, old_step)
old_data = np.random.random((len(old_lon), len(old_lat)))
interp_function = RectBivariateSpline(old_lon, old_lat, old_data, kx=1, ky=1)
new_lon = np.arange(-180, 180, new_step)
new_lat = np.arange(-90, 90, new_step)
new_data = interp_function(new_lon, new_lat)
cartopy
现在有这个名为 cartopy.img_transform.im_trans.regrid
的功能。
例如,
import numpy as np
from cartopy.img_transform import regrid
from cartopy import crs as ccrs
Old_Lon = np.linspace(-180,180,360)
Old_Lat = np.linspace(-90,90,180)
New_Lon = np.linspace(-180,180,720)
New_Lat = np.linspace(-90,90,360)
New_Lon,New_Lat = np.meshgrid(New_Lon,New_Lat)
source_cs = ccrs.PlateCarree()
target_proj = ccrs.PlateCarree()
New_Data = regrid(Old_Data,Old_Lon,Old_Lat, source_cs, target_proj, New_Lon,New_Lat)