使用底图重新网格化 Python 中的 3-D 卫星数据时出错,2-D 工作

Error while regridding 3-D satellite data in Python with Basemap, 2-D works

我有以下问题:

我有来自不同卫星的云层数据集,我想将其重新网格化到气候模型的网格上,以便在模型输出和观测到的卫星数据之间进行比较。

现在我使用的是底图中的 interp 函数,它非常适合以下形状的数组:1 x longitude x latitude,但它不适用于 n x longitude x latitude 形状的数组。重新网格化这些 3-D 阵列的最佳方法是什么?

from mpl_toolkits.basemap import interp
import xarray as xr 
def regrid_MAR10km(x_in,y_in,data_in): 
    mar_10km = xr.open_dataset('/media/..../MAR_10km /MARv3.5.2-10km-ERA-2008.nc')
    lat = mar_10km['LAT']
    lon = mar_10km['LON']
    result = interp(data_in, x_in, y_in,lon,lat)
    return result  

我的问题是,当我尝试使用 3-D 数据时,我收到以下形式的错误消息,在我的例子中,数组的形式是 161(这是每个月的云量!) x经度 x 纬度

<xarray.DataArray 'Cloud_Fraction_Mask_Total_Mean' (Begin_Date: 161, lat: 28, lon: 110)>
[495880 values with dtype=float64]
Coordinates:
* Begin_Date  (Begin_Date) datetime64[ns] 2002-07-01 2002-08-01 2002-09-01 ...
* lat         (lat) float32 85.5 84.5 83.5 82.5 81.5 80.5 79.5 78.5 77.5 ...
* lon         (lon) float32 -94.5 -93.5 -92.5 -91.5 -90.5 -89.5 -88.5 ...

这就是它产生的错误:

 IndexError                                Traceback (most recent call   last)
<ipython-input-19-5117a62e570e> in <module>()
----> 1 cloud_regrid =      fc.regrid_MAR10km(longitude,latitude,cloud_data_UD)

/media/sf_Shared/Black_and_bloom/CODE/functions.py in regrid_MAR10km(x_in, y_in, data_in)
 20         lat = mar_10km['LAT']
 21         lon = mar_10km['LON']
---> 22         result = interp(data_in, x_in, y_in,lon,lat)
 23         return result

/home/sh16450/miniconda3/envs/snowflakes/lib/python3.5/site-packages/mpl_toolkits/basemap/__init__.py in interp(datain, xin, yin, xout, yout,  checkbounds, masked, order)
4958         dataout = (1.-delx)*(1.-dely)*datain[yi,xi] + \
4959                   delx*dely*datain[yip1,xip1] + \
-> 4960                   (1.-delx)*dely*datain[yip1,xi] + \
4961                   delx*(1.-dely)*datain[yi,xip1]
4962     elif order == 0:

IndexError: index 41 is out of bounds for axis 1 with size 28

Basemap docs 开始,interp 方法专门处理二维数组,其中第一个维度对应于 y 维度(纬度),第二个维度对应于 x(经度)。如果你给它任何其他东西,恐怕你最终会得到糟糕的结果。

实际发生的情况是,由于 Basemap 假定您的数据是二维的,因此它会尝试在您的时间和纬度 (161,28) 轴上执行插值,就好像它们分别是纬度和经度一样。由于您传递给该方法的经度数组 (110,) 现在大于轴 Basemap 假定为经度 (28,),因此您最终会导致形状不匹配,从而导致此 IndexError。

我不太确定为什么 (1,28,110) 有效,但我打赌 NumPy 在幕后执行的一些数组 flattening or squeezing 最终基本上忽略了那个空维度。

我建议要么遍历每个时间步长并对二维数据切片执行插值,要么搜索 scipy.interpolate 文档以查看它们是否提供了一个函数来执行您要执行的操作正在寻找。

祝你好运!