import xarray as xr
import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from oggm import utils
import hvplot.pandas
import warnings
warnings.filterwarnings('ignore', r'(invalid value encountered in|Mean of empty slice)')
# Download the data
f_c20c = utils.file_downloader('https://cluster.klima.uni-bremen.de/~fmaussion/climate/cera-20c/cera-20c_t2m_1901-2010.nc')
# Read it with xarray and chunk it for better memory management
ds_c20c_o = xr.open_dataset(f_c20c)
ds_c20c_o = ds_c20c_o.chunk(chunks={'number':1, 'longitude':90})
ds_c20c_o
<xarray.Dataset> Dimensions: (latitude: 181, longitude: 360, number: 10, time: 1320) Coordinates: * longitude (longitude) float32 0.0 1.0 2.0 3.0 ... 356.0 357.0 358.0 359.0 * latitude (latitude) float32 90.0 89.0 88.0 87.0 ... -88.0 -89.0 -90.0 * number (number) int32 0 1 2 3 4 5 6 7 8 9 * time (time) datetime64[ns] 1901-01-01 1901-02-01 ... 2010-12-01 Data variables: t2m (time, number, latitude, longitude) float32 dask.array<shape=(1320, 10, 181, 360), chunksize=(1320, 1, 181, 90)> Attributes: Conventions: CF-1.6 history: 2019-05-26 18:33:08 GMT by grib_to_netcdf-2.12.0: grib_to_n...
# Roll the coords because [0-360] is not nice
ds_c20c = ds_c20c_o.roll(longitude=179, roll_coords=True)
ds_c20c['longitude'] = np.where(ds_c20c['longitude'] > 180, ds_c20c['longitude'] - 360, ds_c20c['longitude'])
ds_c20c = ds_c20c.chunk(chunks={'number':1, 'time':12})
# Convert to C°
ds_c20c['t2m'] = ds_c20c['t2m'] - 273.15
# Read CRU
ds_cru = xr.open_dataset(utils.get_cru_file('tmp', version='4.03')).chunk(chunks={'time':12})
ds_cru
<xarray.Dataset> Dimensions: (lat: 360, lon: 720, time: 1416) Coordinates: * lon (lon) float32 -179.75 -179.25 -178.75 ... 178.75 179.25 179.75 * lat (lat) float32 -89.75 -89.25 -88.75 -88.25 ... 88.75 89.25 89.75 * time (time) datetime64[ns] 1901-01-16 1901-02-15 ... 2018-12-16 Data variables: tmp (time, lat, lon) float32 dask.array<shape=(1416, 360, 720), chunksize=(12, 360, 720)> stn (time, lat, lon) float64 dask.array<shape=(1416, 360, 720), chunksize=(12, 360, 720)> Attributes: Conventions: CF-1.4 title: CRU TS4.03 Mean Temperature institution: Data held at British Atmospheric Data Centre, RAL, UK. source: Run ID = 1905011326. Data generated from:tmp.1905011321.dtb history: Wed 1 May 2019 15:42:51 BST : User ianharris : Program mak... references: Information on the data is available at http://badc.nerc.ac... comment: Access to these data is available to any registered CEDA user. contact: support@ceda.ac.uk
# Rename CRU's coordinates
ds_cru = ds_cru.rename_dims({'lon':'longitude', 'lat':'latitude'}).rename({'lon':'longitude', 'lat':'latitude'})
# CRU's time is horribly centred on month middle
ds_cru['time'] = pd.date_range('1901-01', '2018-12', freq='MS')
# Station network should have missing data as well
ds_cru['stn'] = ds_cru['stn'].where(ds_cru['stn']!=-999)
We don't care about annual cycle for now:
ds_cru = ds_cru.resample(time='AS').mean(dim='time')
ds_c20c = ds_c20c.resample(time='AS').mean(dim='time')
# RGI statistics as csv
frgi = utils.file_downloader('https://cluster.klima.uni-bremen.de/~fmaussion/misc/rgi62_allglaciers.csv')
odf = pd.read_csv(frgi, index_col=0,
converters={'Name':str, 'GLIMSId': str, 'BgnDate':str, 'EndDate':str,
'O1Region': str, 'O2Region':str,
'IsTidewater':bool, 'IsNominal':bool})
No known hash for cluster.klima.uni-bremen.de/~fmaussion/misc/rgi62_allglaciers.csv
lon_bins = np.linspace(-180, 180, 721)
lat_bins = np.linspace(-90., 90, 361)
odf['lon_id'] = np.digitize(odf['CenLon'], lon_bins)-1
odf['lat_id'] = np.digitize(odf['CenLat'], lat_bins)-1
odf['unique_id'] = ['{:03d}_{:03d}'.format(lon, lat) for (lon, lat) in zip(odf['lon_id'], odf['lat_id'])]
mdf = odf.drop_duplicates(subset='unique_id').set_index('unique_id')
mdf['Area'] = odf.groupby('unique_id').sum()['Area']
print(len(mdf))
4706
mask = np.zeros((360, 720), dtype=float) * np.NaN
mask[mdf['lat_id'], mdf['lon_id']] = mdf['Area']
# Now CRU doesn't have data everywhere - Only keep glaciers where we also have data
mask = np.where(np.isfinite(ds_cru.tmp.isel(time=0)), mask, np.NaN)
ds_cru['glacier_mask'] = (('latitude', 'longitude'), np.isfinite(mask))
ds_cru['glacier_area'] = (('latitude', 'longitude'), mask)
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
ds_cru['glacier_area'].plot(ax=ax, transform=ccrs.PlateCarree())
ax.coastlines(); ax.gridlines();
odf = pd.read_csv(frgi, index_col=0,
converters={'Name':str, 'GLIMSId': str, 'BgnDate':str, 'EndDate':str,
'O1Region': str, 'O2Region':str,
'IsTidewater':bool, 'IsNominal':bool})
lon_bins = np.linspace(-179.5, 180.5, 361)
lat_bins = np.linspace(90.5, -90.5, 182)
odf['lon_id'] = np.digitize(odf['CenLon'], lon_bins)-1
odf['lat_id'] = np.digitize(odf['CenLat'], lat_bins)-1
odf['unique_id'] = ['{:03d}_{:03d}'.format(lon, lat) for (lon, lat) in zip(odf['lon_id'], odf['lat_id'])]
mdf = odf.drop_duplicates(subset='unique_id').set_index('unique_id')
mdf['Area'] = odf.groupby('unique_id').sum()['Area']
print(len(mdf))
1976
mask = np.zeros((181, 360), dtype=float) * np.NaN
mask[mdf['lat_id'], mdf['lon_id']] = mdf['Area']
ds_c20c['glacier_mask'] = (('latitude', 'longitude'), np.isfinite(mask))
ds_c20c['glacier_area'] = (('latitude', 'longitude'), mask)
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
ds_c20c['glacier_area'].plot(ax=ax, transform=ccrs.PlateCarree())
ax.coastlines(); ax.gridlines();
We also compute the weight of glaciated cells (weighted per glaciated area).
# Weight
weight = np.cos(np.deg2rad(ds_cru.latitude))
weight = ds_cru.tmp.isel(time=0) * 0. + weight
ds_cru['weight'] = (('latitude', 'longitude'), weight / weight.sum())
ds_cru['weight'].plot();
# For glaciers we don't care about sphere, just about area per grid point
weight = ds_cru['weight'] * 0 + ds_cru['glacier_area']
ds_cru['weight_glacier'] = (('latitude', 'longitude'), weight / weight.sum())
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
(ds_cru['weight_glacier'] * 100).plot(ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={'label':'Glacier area (%)'})
ax.coastlines(); ax.gridlines(); plt.title('Glacier area per grid point (% of total)');
plt.savefig('cru_glacier_mask.png', dpi=150, bbox_inches='tight');
f, ax = plt.subplots(figsize=(4, 4))
(ds_cru['weight_glacier'] * 100).sum(dim='longitude').plot(y='latitude', ax=ax);
plt.ylim(-90, 90); plt.ylabel('Latitude'); plt.xlabel('Glacier area (% of total)');
plt.title('Zonal average of glacier area');
plt.savefig('zonal_avg.png', dpi=150, bbox_inches='tight');
weight = np.cos(np.deg2rad(ds_c20c.latitude)).clip(0)
weight = ds_c20c.t2m.isel(time=0, number=0) * 0. + weight
ds_c20c['weight'] = (('latitude', 'longitude'), weight / weight.sum())
ds_c20c['weight'].plot();
weight = ds_c20c['weight'] * 0 + ds_c20c['glacier_area']
ds_c20c['weight_glacier'] = (('latitude', 'longitude'), weight / weight.sum())
ds_c20c['weight_glacier'].plot();
In order to really compare CERA with CRU, we need to interpolate CERA on CRU - there is no way around it. This is a bit memory expensive.
# But we don't bother with ensemble here
c20c_t2m = ds_c20c.t2m.mean(dim='number').compute()
# Interpolate
ds_cru['cera_on_cru'] = c20c_t2m.interp(longitude=ds_cru.longitude, latitude=ds_cru.latitude)
ds_cru['cera_on_cru'] = ds_cru['cera_on_cru'].where(np.isfinite(ds_cru.tmp.isel(time=-1)))
OK! Now the fun things, finally!
The comparison is only partly meaningful because of the different reference topographies. For this reason we will stop working in absolute temp afterwards, but only talk about changes from the reference period '1961', '1990'.
ds_cru_ref = ds_cru.sel(time=slice('1961', '1990')).mean(dim='time')
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
ds_cru_ref['tmp'].plot(ax=ax, transform=ccrs.PlateCarree(), vmin=-31, vmax=31, cmap='RdBu_r', cbar_kwargs={'label':'°C'});
ax.coastlines(); ax.gridlines(); plt.title('Average temperature [1961-1990] CRU TS 4.03');
plt.savefig('refperiod_cru_map.png', dpi=150, bbox_inches='tight');
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
ds_cru_ref['cera_on_cru'].plot(ax=ax, transform=ccrs.PlateCarree(), vmin=-31, vmax=31, cmap='RdBu_r', cbar_kwargs={'label':'°C'})
ax.coastlines(); ax.gridlines(); plt.title('Average temperature [1961-1990] CERA-20C');
plt.savefig('refperiod_cera_map.png', dpi=150, bbox_inches='tight');
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
(ds_cru_ref['tmp'] - ds_cru_ref['cera_on_cru']).plot(ax=ax, transform=ccrs.PlateCarree(),
vmin=-12, vmax=12, cmap='RdBu_r', cbar_kwargs={'label':'°C'})
ax.coastlines(); ax.gridlines(); plt.title('Diff CRU - CERA-20C; average temperature [1961-1990]');
plt.savefig('refperiod_diff_cru_cera_map.png', dpi=150, bbox_inches='tight');
(ds_cru['tmp'] * ds_cru.weight).sum(dim=['longitude','latitude']).plot(label='CRU');
(ds_cru['cera_on_cru'] * ds_cru.weight).sum(dim=['longitude','latitude']).plot(label='CERA-20c on CRU');
plt.xlim('1901', '2009'); plt.ylim(11, 15); plt.legend(loc='best');
ds_cru['tmp'] = ds_cru['tmp'] - ds_cru_ref['tmp']
ds_cru['cera_on_cru'] = ds_cru['cera_on_cru'] - ds_cru_ref['cera_on_cru']
ts_cru_all = (ds_cru['tmp'] * ds_cru.weight).sum(dim=['longitude','latitude'])
df = ts_cru_all.to_pandas().to_frame(name='cru_world')
df['cera_on_cru_world'] = (ds_cru['cera_on_cru'] * ds_cru.weight).sum(dim=['longitude','latitude'])
df['cera_on_cru_world'] = np.where(df['cera_on_cru_world'] == 0, np.NaN, df['cera_on_cru_world'])
df.plot(color=['k', 'C0']);
plt.legend(['CRU TS 4.01', 'CERA-20C']);
plt.title('Average land temperatures, anomaly to 1961-1990');
plt.ylabel('°C'); plt.xlabel('');
plt.savefig('global_land_ts.png', dpi=150, bbox_inches='tight');
Much better!
df['cru_glaciers'] = (ds_cru['tmp'] * ds_cru.weight_glacier).sum(dim=['longitude','latitude'])
df['cera_on_cru_glaciers'] = (ds_cru['cera_on_cru'] * ds_cru.weight_glacier).sum(dim=['longitude','latitude'])
df['cera_on_cru_glaciers'] = np.where(df['cera_on_cru_glaciers'] == 0, np.NaN, df['cera_on_cru_glaciers'])
f, axs = plt.subplots(2, 2, figsize=(16, 9))
df[['cru_world', 'cera_on_cru_world']].plot(ax=axs.flatten()[0]);
df[['cru_glaciers', 'cera_on_cru_glaciers']].plot(ax=axs.flatten()[1]);
df[['cru_world', 'cru_glaciers']].plot(ax=axs.flatten()[2]);
df[['cera_on_cru_world', 'cera_on_cru_glaciers']].plot(ax=axs.flatten()[3]);
df[['cru_world', 'cru_glaciers']].plot(color=['k', 'grey']);
plt.legend(['CRU TS 4.01 (all land)', 'CRU TS 4.01 (glaciers)']);
plt.title('Average temperatures, anomaly to 1961-1990');
plt.ylabel('°C'); plt.xlabel(''); plt.ylim(-1.7, 2.8);
plt.savefig('global_cru_glacier_ts.png', dpi=150, bbox_inches='tight');
df[['cera_on_cru_world', 'cera_on_cru_glaciers']].plot(color=['C0', 'lightsteelblue']);
plt.legend(['CERA-20C (all land)', 'CERA-20C (glaciers)']);
plt.title('Average temperatures, anomaly to 1961-1990');
plt.ylabel('°C'); plt.xlabel(''); plt.ylim(-1.7, 2.8);
plt.savefig('global_cera_glacier_ts.png', dpi=150, bbox_inches='tight');
df[['cru_glaciers', 'cera_on_cru_glaciers']].plot(color=['grey', 'lightsteelblue']);
plt.legend(['CRU (glaciers)', 'CERA-20C (glaciers)']);
plt.title('Average temperatures, anomaly to 1961-1990');
plt.ylabel('°C'); plt.xlabel(''); plt.ylim(-1.7, 2.8);
plt.savefig('global_cru_cera_glacier_ts.png', dpi=150, bbox_inches='tight');
ds_cru_1920 = ds_cru.sel(time=slice('1921', '1950')).mean(dim='time')
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
ds_cru_1920['tmp'].plot(ax=ax, transform=ccrs.PlateCarree(),
vmin=-1.8, vmax=1.8, cmap='RdBu_r', cbar_kwargs={'label':'°C'})
ax.coastlines(); ax.gridlines(); plt.title('CRU: 1921-1950 temperature anomaly to 1961-1990');
plt.savefig('1920_cru_map.png', dpi=150, bbox_inches='tight');
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
ds_cru_1920['cera_on_cru'].plot(ax=ax, transform=ccrs.PlateCarree(),
vmin=-1.8, vmax=1.8, cmap='RdBu_r', cbar_kwargs={'label':'°C'})
ax.coastlines(); ax.gridlines(); plt.title('CERA: 1921-1950 temperature anomaly to 1961-1990');
plt.savefig('1920_cera_map.png', dpi=150, bbox_inches='tight');
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
(ds_cru_1920['tmp'] - ds_cru_1920['cera_on_cru']).plot(ax=ax, transform=ccrs.PlateCarree())
ax.coastlines(); ax.gridlines();
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
(ds_cru_1920['stn']).plot(ax=ax, transform=ccrs.PlateCarree(), cbar_kwargs={'label':'N stations (from 0 to 8)'})
ax.coastlines(); ax.gridlines(); plt.title('CRU: number of stations used in the interpolation 1921-2950');
plt.savefig('1920_cru_nstats.png', dpi=150, bbox_inches='tight');