matplotlib 底图插值备选方案-重新网格化数据

ffx8fchx  于 2023-10-24  发布在  其他
关注(0)|答案(3)|浏览(92)

我将从底图转向制图,因为底图将被淘汰。我以前使用过basemap.interp功能来插值数据,例如,假设我有1度分辨率(180 x360)的数据,我将运行以下命令来插值到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有插值,但我不知道如何应用它。任何帮助都将不胜感激!

4si2a6ki

4si2a6ki1#

我最终决定从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

--

rhfm7lfc

rhfm7lfc2#

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)
a5g8bdjr

a5g8bdjr3#

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)

相关问题