Блог

Я интерполирую данные из океанического компонента моделей CMIP6 в сетку 2×2. Поле имеет значения dim (время, nav_lat, nav_lon) и nan в континенте. Здесь nav_lon и nav_lat представляют собой двумерную криволинейную сетку. Я могу выполнить интерполяцию с помощью griddata из scipy, но мне придется использовать цикл с течением времени. Цикл делает его довольно медленным, если данные содержат тысячи записей о времени. Мой вопрос в том, как векторизовать интерполяцию с течением времени.

Ниже приведен мой код:

 import xarray as xr
import numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt

source = xr.open_dataset('data/zos_2850.nc',decode_times=False)

# obtain old lon and lat (and put lon in 0-360 range)
# nav_lon is from -180 to 180, not in 0-360 range
loni, lati = source.nav_lon.values60, source.nav_lat.values

# flatten the source coordinates
loni_flat, lati_flat = loni.flatten(), lati.flatten()

# define a 2x2 lon-lat grid
lon, lat = np.linspace(0,360,181), np.linspace(-90,90,91) 

# create mesh
X, Y = np.meshgrid(lon,lat)

# loop over time
ntime = len(source.time)
tmp = []
for t in range(ntime):  
    print(t)
    var_s = source.zos[t].values
    var_s_flat = var_s.flatten()

    # index indicates where they are valid values
    index = np.where(~np.isnan(var_s_flat))
    # remap the valid values to the new grid 
    var_t = griddata((loni_flat[index],lati_flat[index]),var_s_flat[index], (X,Y), 
    method='cubic')
    
    # interpolate mask using nearest
    maskinterp = griddata((loni_flat,lati_flat),var_s_flat, (X,Y), method='nearest')
    # re-mask interpolated data
    var_t[np.isnan(maskinterp)] = np.nan

tmp.append(var_t)

# convert to DataArray
da = xr.DataArray(data=tmp, 
              dims=["time","lat","lon"], 
              coords=dict(lon=(["lon"], lon),lat=(["lat"], lat),time=source['time']))