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
import intake
import cftime
import 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 = ""
col = intake.open_esm_datastore(cat_url)

Get data in xarray#

Search od550aer variable dataset#

cat =['historical'], variable_id='od550aer', member_id=['r1i1p1f1'], grid_label='gn')

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})

Select model and visualize a single date#

  • Use data as xarray to make a simple plot

ds = dset_dict['']

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())
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)#


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()
# create dictionary for reggridded data
ds_regrid_dict = dict()
for key in dset_dict.keys():
    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 = ''
    savepath = 'CMIP6_hist/{}'.format(model)
    nc_out = os.path.join(savepath, filename)
    os.makedirs(savepath, exist_ok=True) 
    # 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'})

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['od550aer'].min().compute(), ds_seas['od550aer'].max().compute()

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

!put -p CMIP6_hist/ -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.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)