1

我有以下问题：

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

现在我正在使用底图中的 interp 函数，它对以下形状的数组非常有效：1 x 经度 x 纬度，但它不适用于形状为 nx 经度 x 纬度的数组。重新排列这些 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 lon x lat

<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
4

1 回答 1

0

从Basemap docs interp 方法专门处理二维数组，其中第一个维度对应于 y 维度（纬度），第二个维度对应于 x（经度）。如果你给它其他任何东西，我担心你最终会得到糟糕的结果。

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

我不确定 (1,28,110) 为何有效，但我敢打赌，NumPy 在幕后会进行一些数组展平或挤压，最终基本上忽略了那个空维度。

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

祝你好运！

于 2016-05-12T17:52:45.120 回答