Regridding model data with xESMF#

Import python packages#

# supress warnings
import warnings
warnings.filterwarnings('ignore') # don't output warnings

import os
# import packages
import xarray as xr
xr.set_options(display_style='html')
import intake
import cftime
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib.path as mpath
import numpy as np
import xesmf as xe
from cmcrameri import cm

%matplotlib inline

Open CMIP6 online catalog#

cat_url = "https://storage.googleapis.com/cmip6/pangeo-cmip6.json"
col = intake.open_esm_datastore(cat_url)
col

Get data in xarray#

Search od550aer variable dataset#

cat = col.search(experiment_id=['historical'], variable_id='od550aer', member_id=['r1i1p1f1'], grid_label='gn')
cat.df

cat.df['source_id'].unique()

Create dictionary from the list of datasets we found#

  • This step may take several minutes so be patient!

dset_dict = cat.to_dataset_dict(zarr_kwargs={'use_cftime':True})

list(dset_dict.keys())

Select model and visualize a single date#

  • Use data as xarray to make a simple plot

ds = dset_dict['CMIP.NCC.NorESM2-LM.historical.AERmon.gn']
ds

Plot on NorthPolarStereo and set the latitude limit#

def polarCentral_set_latlim(lat_lims, ax):
    ax.set_extent([-180, 180, lat_lims[0], lat_lims[1]], ccrs.PlateCarree())
    # Compute a circle in axes coordinates, which we can use as a boundary
    # for the map. We can pan/zoom as much as we like - the boundary will be
    # permanently circular.
    theta = np.linspace(0, 2*np.pi, 100)
    center, radius = [0.5, 0.5], 0.5
    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
    circle = mpath.Path(verts * radius + center)

    ax.set_boundary(circle, transform=ax.transAxes)

fig = plt.figure(1, figsize=[10,10])

# Fix extent
minval = 0
maxval = 0.3

ax = plt.subplot(1, 1, 1, projection=ccrs.NorthPolarStereo())
ax.coastlines()
ax.gridlines()
polarCentral_set_latlim([50,90], ax)
ds['od550aer'].sel(time=cftime.DatetimeNoLeap(1985, 1, 16, 12, 0, 0, 0)).plot(ax=ax, vmin=minval, vmax=maxval, transform=ccrs.PlateCarree(), cmap=cm.oslo_r)

Get attributes (unique identifier)#

ds.attrs['tracking_id']

Regrid CMIP6 data to common NorESM2-LM grid#

  • Select a time range

  • we use squeeze to remove dimension with one element only e.g. member_id=’r1i1p1f1’

starty = 1985; endy = 1986
year_range = range(starty, endy+1)

# Read in the output grid from NorESM
ds_out = ds.sel(time = ds.time.dt.year.isin(year_range)).squeeze()
ds_out

# create dictionary for reggridded data
ds_regrid_dict = dict()
for key in dset_dict.keys():
    print(key)
    ds_in = dset_dict[keys]
    ds_in = ds_in.sel(time = ds_in.time.dt.year.isin(year_range)).squeeze()
    regridder = xe.Regridder(ds_in, ds_out, 'bilinear')
    # Apply regridder to data
    # the entire dataset can be processed at once
    ds_in_regrid = regridder(ds_in, keep_attrs=True)
    # Save to netcdf file
    model = key.split('.')[2]
    filename = 'od550aer_AERmon.nc'
    savepath = 'CMIP6_hist/{}'.format(model)
    nc_out = os.path.join(savepath, filename)
    os.makedirs(savepath, exist_ok=True) 
    ds_in_regrid.to_netcdf(nc_out)
    # create dataset with all models
    ds_regrid_dict[model] = ds_in_regrid
    print('file written: {}'.format(nc_out))

Concatenate all models#

_ds = list(ds_regrid_dict.values())
_coord = list(ds_regrid_dict.keys())
ds_out_regrid = xr.concat(objs=_ds, dim=_coord, coords="all").rename({'concat_dim':'model'})
ds_out_regrid

Compute seasonal mean of all regridded models#

ds_seas = ds_out_regrid.mean('model', keep_attrs=True, skipna = True).groupby('time.season').mean('time', keep_attrs=True, skipna = True)
ds_seas

ds_seas['od550aer'].min().compute(), ds_seas['od550aer'].max().compute()

Save seasonal mean in a new netCDF file and in the current Galaxy history#

ds_seas.to_netcdf('CMIP6_hist/od550aer_seasonal.nc')

!put -p CMIP6_hist/od550aer_seasonal.nc -t netcdf

Visualize final results (seasonal mean for all models)#

import matplotlib
proj_plot = ccrs.Mercator()

p = ds_seas['od550aer'].plot(x='lon', y='lat', transform=ccrs.PlateCarree(),
                             aspect=ds_seas.dims["lon"] / ds_seas.dims["lat"],  # for a sensible figsize
                             subplot_kws={"projection": proj_plot},
                             col='season', col_wrap=2, robust=True, cmap='PiYG')
# We have to set the map's options on all four axes
for ax,i in zip(p.axes.flat,  ds_seas.season.values):
    ax.coastlines()
    ax.set_title('Season '+i, fontsize=18)
fig = matplotlib.pyplot.gcf()
fig.set_size_inches(18.5, 10.5)
fig.savefig('od550aer_seasonal_mean.png', dpi=100)