In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import xarray as xr
import xesmf as xe
In [2]:
infile = '/g/data/r87/dbi599/temp/thetao_Omon_IPSL-CM5A-LR_historical_r1i1p1_1850-01-01.nc'
#infile = '/g/data/ua6/DRSv3/CMIP5/IPSL-CM5A-LR/historical/mon/ocean/r1i1p1/thetao/latest/thetao_Omon_IPSL-CM5A-LR_historical_r1i1p1_185001-189912.nc'
In [3]:
ds = xr.open_dataset(infile, decode_times=False)
ds
Out[3]:
In [4]:
plt.figure(figsize=(12,2));
ax = plt.axes(projection=ccrs.PlateCarree());
ds['thetao'][0, 0, ::].plot.pcolormesh(ax=ax, x='lon', y='lat');
ax.coastlines();
In [10]:
ds_out = xe.util.grid_global(5, 4)
#ds_out = xr.Dataset({'lat': (['lat'], np.arange(-80, 81, 1)),
# 'lon': (['lon'], np.arange(2, 362, 2)),
# }
# )
In [6]:
regridder = xe.Regridder(ds, ds_out, 'bilinear')
The regular approach to remapping fails...
In [5]:
ds['lon']
Out[5]:
In [6]:
# multiple titles are stacked into a single 2D array
# just passinsg this will crash ESMPy
plt.scatter(ds['lon'], ds['lat'], s=0.2)
Out[6]:
In [7]:
# Get a more well-defined 2D mesh (subset of the full grid)
plt.scatter(ds['lon'][:80,:], ds['lat'][:80,], s=0.2)
Out[7]:
In [8]:
ds_subset = ds.isel(i=slice(0,80), j=slice(0,80))
In [11]:
regridder = xe.Regridder(ds_subset, ds_out, 'bilinear')
In [13]:
dr_out = regridder(ds['thetao'][:, :, 0:80, 0:80])
In [16]:
dr_out
Out[16]:
In [17]:
plt.figure(figsize=(12,2));
ax = plt.axes(projection=ccrs.PlateCarree());
dr_out[0, 0, ::].plot.pcolormesh(ax=ax, x='lon', y='lat');
ax.coastlines();
So the error could be fixed by breaking your full grid to several well-defined 2D tiles... or maybe it's easier to just use CDO...
In [19]:
cdo_infile = '/g/data/r87/dbi599/temp/thetao_Omon_IPSL-CM5A-LR_historical_r1i1p1_1850-01-01_susan-grid-cdo.nc'
In [20]:
ds_cdo = xr.open_dataset(cdo_infile, decode_times=False)
ds_cdo
Out[20]:
In [21]:
plt.figure(figsize=(12,2));
ax = plt.axes(projection=ccrs.PlateCarree());
ds_cdo['thetao'][0, 0, ::].plot.pcolormesh(ax=ax);
ax.coastlines();
In [29]:
ds_cdo['lev'].data[0]
Out[29]:
In [35]:
test_file = '/g/data/r87/dbi599/DRSv2/CMIP5/IPSL-CM5A-LR/historical/mon/ocean/r1i1p1/thetao/latest/thetao_Omon_IPSL-CM5A-LR_historical_r1i1p1_185001-189912_susan-horiz-grid.nc'
In [36]:
ds_test = xr.open_dataset(test_file, decode_times=False)
ds_test
Out[36]:
In [37]:
plt.figure(figsize=(12,2));
ax = plt.axes(projection=ccrs.PlateCarree());
ds_test['thetao'][0, 0, ::].plot.pcolormesh(ax=ax);
ax.coastlines();
In [38]:
type(ds_test['thetao'].data)
Out[38]:
In [39]:
ds_test['thetao'].
In [ ]: