Ошибка при regridding спутниковых данных 3-D в Python с BaseMap, 2-D работает

Question

У меня есть следующий вопрос:

У меня есть облачный покров наборов данных из разных спутников, что я хочу regrid на сетку климата модель для сравнения между выходом модели и наблюдаемыми спутниковыми данными.

На данный момент я использовал функцию interp из baseemap, которая отлично работает для массивов в форме: 1 x долгота x широты, но она не работает для массивов формы n x долготы 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 (который облачный покров за каждый месяц!) х долгота х ш

<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 

Answer 1

Из Basemap docs на метод интерполяцией сделки исключительно с 2-D массивов, где первое измерение соответствует ярд (широты) и от 2 до x (долготы). Если вы дадите что-нибудь еще, я боюсь, что вы получите плохие результаты.

На самом деле происходит то, что, поскольку Basemap предполагает, что ваши данные являются 2-D, он пытается выполнить интерполяцию на оси времени и широты (161,28), как если бы они были широтой и долготой соответственно. Поскольку массив долготы (110), который вы передали методу, теперь больше, чем ось. Базовая карта предполагает, что это долготы (28), вы оказываете несоответствие формы, которое приводит к этому индексу.

Я не совсем уверен, почему работает (1,28,110), но я утверждаю, что есть какой-то массив flattening или squeezing, который NumPy делает под обложками, что в итоге игнорирует это пустое измерение.

Я рекомендую либо выполнять итерацию через каждый из этапов времени, и выполнять интерполяцию на этом двумерном фрагменте данных, либо искать в документах scipy.interpolate, чтобы узнать, предлагают ли они функцию для выполнения того, что вы ищете.

Удачи вам!