Code and data used to generate the plots shown in Fabien's post: https://fabienmaussion.info/2019/08/29/era5
# Data wrangling libs
import xarray as xr
import numpy as np
import pandas as pd
from scipy import stats
# Just for download and fancy bar
from oggm import utils
import progressbar
# Plots
import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import seaborn as sns
Data downloaded from ECMWF and stored here for convenience:
f_era = utils.file_downloader('https://cluster.klima.uni-bremen.de/~fmaussion/climate/era5/monthly/v1.0/era5_monthly_t2m_1979-2018.nc')
ds = xr.open_dataset(f_era).chunk(chunks={'time':12})
ds
No known hash for cluster.klima.uni-bremen.de/~fmaussion/climate/era5/monthly/v1.0/era5_monthly_t2m_1979-2018.nc
<xarray.Dataset> Dimensions: (latitude: 721, longitude: 1440, time: 480) Coordinates: * longitude (longitude) float32 0.0 0.25 0.5 0.75 ... 359.25 359.5 359.75 * latitude (latitude) float32 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0 * time (time) datetime64[ns] 1979-01-01 1979-02-01 ... 2018-12-01 Data variables: t2m (time, latitude, longitude) float32 dask.array<shape=(480, 721, 1440), chunksize=(12, 721, 1440)> Attributes: Conventions: CF-1.6 history: 2019-05-26 11:47:21 GMT by grib_to_netcdf-2.10.0: /opt/ecmw...
We work we annual averages only:
ds = ds.resample(time='AS').mean(dim='time')
For analysis of changes, we compute the anomaly to the first decade in the dataset as reference.
ds['t2m_ano'] = ds['t2m'] - ds['t2m'].sel(time=slice('1979', '1988')).mean(dim='time')
f_era = utils.file_downloader('https://cluster.klima.uni-bremen.de/~fmaussion/climate/era5/monthly/v1.0/era5_invariant.nc')
ds_inv = xr.open_dataset(f_era).isel(time=0)
ds_inv
No known hash for cluster.klima.uni-bremen.de/~fmaussion/climate/era5/monthly/v1.0/era5_invariant.nc
<xarray.Dataset> Dimensions: (latitude: 721, longitude: 1440) Coordinates: * longitude (longitude) float32 0.0 0.25 0.5 0.75 ... 359.25 359.5 359.75 * latitude (latitude) float32 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0 time datetime64[ns] 1979-01-01 Data variables: lsm (latitude, longitude) float32 ... z (latitude, longitude) float32 ... Attributes: Conventions: CF-1.6 history: 2019-08-19 17:31:31 GMT by grib_to_netcdf-2.10.0: /opt/ecmw...
# Convert geop to elevation
ds_inv['z'] = (ds_inv['z'] / 9.80665)
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
ds_inv['z'].plot(ax=ax, transform=ccrs.PlateCarree(), cmap='terrain', vmin=0, cbar_kwargs={'label':'Elevation (m a.s.l.)'})
ax.coastlines(); ax.gridlines(); ax.set_title('ERA5 elevation');
plt.savefig('ERA5/era5_map_elevation.png', dpi=150, bbox_inches='tight');
plt.savefig('ERA5/era5_map_elevation_l.png', dpi=300, bbox_inches='tight');
We download data from the Randolph Glacier Inventory (RGI 6.0), pre-processed in tabular form:
# 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
Some data wrangling to locate the ERA5 grid points where we have glaciers:
nx, ny = ds.dims['longitude'], ds.dims['latitude']
# Nearest neighbor lookup
cenlon_for_bins = np.where(odf['CenLon'] < -0.125, odf['CenLon']+360, odf['CenLon'])
lon_bins = np.linspace(-0.125, 359.75+0.125, nx)
lat_bins = np.linspace(90+0.125, -90-0.125, ny)
odf['lon_id'] = np.digitize(cenlon_for_bins, lon_bins)-1
odf['lat_id'] = np.digitize(odf['CenLat'], lat_bins)-1
# Use unique grid points as index and compute the area per location
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('Total number of glaciers: {} and number of ERA5 gridpoints with glaciers in them: {}'.format(len(odf), len(mdf)))
Total number of glaciers: 216502 and number of ERA5 gridpoints with glaciers in them: 11323
Glacier locations and corresponding weights:
mask = np.full((ny, nx), np.NaN)
mask[mdf['lat_id'], mdf['lon_id']] = mdf['Area']
ds['glacier_mask'] = (('latitude', 'longitude'), np.isfinite(mask))
ds['glacier_area'] = (('latitude', 'longitude'), mask)
ds['weight_glacier'] = (('latitude', 'longitude'), mask / np.nansum(mask))
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
ds['glacier_area'].plot(ax=ax, transform=ccrs.PlateCarree(), vmax=500,
cbar_kwargs={'label':'Glacier area (km$^{2}$)'})
ax.coastlines(); ax.gridlines(); ax.set_title('Glacier area per ERA5 grid point (km$^{2}$)');
plt.savefig('ERA5/era5_map_area.png', dpi=150, bbox_inches='tight');
plt.savefig('ERA5/era5_map_area_l.png', dpi=300, bbox_inches='tight');
The Earth is not flat: take this into account as well:
# Weight
weight = np.cos(np.deg2rad(ds.latitude)).clip(0)
weight = ds.t2m.isel(time=0) * 0. + weight
ds['weight'] = (('latitude', 'longitude'), weight / weight.sum())
weight = ds['weight'].where(ds_inv['lsm'] > 0.5)
ds['weight_land'] = weight / weight.sum().data
weight = ds['weight'].where(ds_inv['lsm'] < 0.5)
ds['weight_ocean'] = weight / weight.sum().data
# Data containers
trends = np.full((ny, nx), np.NaN)
p_values = np.full((ny, nx), np.NaN)
trends_std_err = np.full((ny, nx), np.NaN)
# Things we are going to regress onto
x = np.asarray(ds['time.year'])
y = ds['t2m'].values
# Good old for loop
for j in progressbar.progressbar(range(ny)):
for i in range(nx):
slope, intercept, r_value, p_value, std_err = stats.linregress(x, y[:, j, i])
trends[j, i] = slope
p_values[j, i] = p_value
trends_std_err[j, i] = std_err
# Back to xarray
ds['trends'] = (('latitude', 'longitude'), trends * 10)
ds['trends_std_err'] = (('latitude', 'longitude'), trends_std_err * 10)
ds['p_values'] = (('latitude', 'longitude'), p_values)
100% (721 of 721) |######################| Elapsed Time: 0:03:13 Time: 0:03:13
fig = plt.figure(figsize=(12, 5))
ax = plt.axes(projection=ccrs.Robinson())
ds['trends'].plot(ax=ax, transform=ccrs.PlateCarree(), vmax=1.5,
cbar_kwargs={'label':'Temperature trend (°C decade$^{-1}$)'})
ax.coastlines(); ax.gridlines(); ax.set_title('ERA5 40yr temperature trend (°C per decade)');
plt.savefig('ERA5/era5_map_trends.png', dpi=150, bbox_inches='tight');
plt.savefig('ERA5/era5_map_trends_l.png', dpi=300, bbox_inches='tight');
ts_global = (ds['t2m_ano'] * ds.weight).sum(dim=['longitude','latitude']).to_dataframe(name='Global')
ts_global['Glaciers'] = (ds['t2m_ano'] * ds.weight_glacier).sum(dim=['longitude','latitude']).to_series()
ts_global.plot();
f, ax = plt.subplots()
ax.plot(ts_global.index.year, ts_global.Global, color='C0', label='Global')
sns.regplot(x=np.array(ts_global.index.year), y=ts_global.Global, color='C0', ax=ax, truncate=True)
ax.plot(ts_global.index.year, ts_global.Glaciers, color='C3', label='Glaciers')
sns.regplot(x=np.array(ts_global.index.year), y=ts_global.Glaciers, color='C3', ax=ax, truncate=True)
plt.legend(loc='best'); plt.ylabel('Temperature anomaly to 1979-88 (K)');
plt.savefig('ERA5/era5_ts_global_glaciers.png', dpi=150, bbox_inches='tight');
lin1 = stats.linregress(ts_global.index.year, ts_global['Global'])
lin2 = stats.linregress(ts_global.index.year, ts_global['Glaciers'])
print('Global warming trend is {:.2f} ± {:.2f} K per decade, '
'Glacier trend is {:.2f} ± {:.2f} K per decade'.format(lin1[0]*10, lin1[-1]*10, lin2[0]*10, lin2[-1]*10))
print('This is a factor of {:.2f}'.format(lin2[0]/lin1[0]))
Global warming trend is 0.18 ± 0.02 K per decade, Glacier trend is 0.44 ± 0.05 K per decade This is a factor of 2.41
ts_global['Land'] = (ds['t2m_ano'] * ds.weight_land).sum(dim=['longitude','latitude']).to_series()
f, ax = plt.subplots()
ax.plot(ts_global.index.year, ts_global.Land, color='C2', label='Land areas')
sns.regplot(x=np.array(ts_global.index.year), y=ts_global.Land, color='C2', ax=ax, truncate=True)
ax.plot(ts_global.index.year, ts_global.Glaciers, color='C3', label='Glaciers')
sns.regplot(x=np.array(ts_global.index.year), y=ts_global.Glaciers, color='C3', ax=ax, truncate=True)
plt.legend(loc='best'); plt.ylabel('Temperature anomaly to 1979-88 (K)');
plt.savefig('ERA5/era5_ts_land_glaciers.png', dpi=150, bbox_inches='tight');
lin1 = stats.linregress(ts_global.index.year, ts_global['Land'])
lin2 = stats.linregress(ts_global.index.year, ts_global['Glaciers'])
print('Land areas warming trend is {:.2f} ± {:.2f} K per decade, '
'Glacier trend is {:.2f} ± {:.2f} K per decade'.format(lin1[0]*10, lin1[-1]*10, lin2[0]*10, lin2[-1]*10))
print('This is a factor of {:.2f}'.format(lin2[0]/lin1[0]))
Land areas warming trend is 0.30 ± 0.03 K per decade, Glacier trend is 0.44 ± 0.05 K per decade This is a factor of 1.46
df_trends = ds.trends.mean(dim='longitude').to_series().to_frame(name='All')
df_trends['Land'] = ds.trends.where(ds_inv.lsm > 0.5).mean(dim='longitude').to_series()
df_trends['Ocean'] = ds.trends.where(ds_inv.lsm < 0.5).mean(dim='longitude').to_series()
df_trends['Glaciers'] = ds.trends.where(ds['weight_glacier'] > 0).mean(dim='longitude').to_series()
df_trends['Above 1000'] = ds.trends.where(ds_inv['z'] > 1000).mean(dim='longitude').to_series()
df_trends['Above 2000'] = ds.trends.where(ds_inv['z'] > 2000).mean(dim='longitude').to_series()
df_trends['Global average'] = (ds['trends'] * ds.weight).sum(dim=['longitude','latitude']).values
df_trends['Land average'] = (ds['trends'] * ds.weight_land).sum(dim=['longitude','latitude']).values
df_trends['Ocean average'] = (ds['trends'] * ds.weight_ocean).sum(dim=['longitude','latitude']).values
df_trends['Glacier average'] = (ds['trends'] * ds.weight_glacier).sum(dim=['longitude','latitude']).values
df_trends['Weight All'] = ds.weight.sum(dim='longitude').to_series()
df_trends['Weight Land'] = ds.weight_land.sum(dim='longitude').to_series()
df_trends['Weight Ocean'] = ds.weight_ocean.sum(dim='longitude').to_series()
df_trends['Weight Glaciers'] = ds.weight_glacier.sum(dim='longitude').to_series()
/home/mowglie/.pyvirtualenvs/py3/lib/python3.5/site-packages/xarray/core/nanops.py:160: RuntimeWarning: Mean of empty slice return np.nanmean(a, axis=axis, dtype=dtype)
f, ax = plt.subplots(figsize=(16, 7))
df_trends[['All', 'Land', 'Ocean']].plot(ax=ax);
f, ax = plt.subplots(figsize=(16, 7))
df_trends[['Land', 'Glaciers', 'Above 2000']].plot(ax=ax);
f, ax = plt.subplots(figsize=(16, 7))
df_trends[['Weight All', 'Weight Land', 'Weight Ocean', 'Weight Glaciers']].plot(ax=ax);
We smooth just a little bit:
dfs = df_trends.rolling(5, center=True).mean()
dfs['Weight All'] = df_trends['Weight All']
with sns.axes_style('ticks'):
f, ax = plt.subplots()
dfs[['All', 'Land', 'Ocean']].plot(ax=ax, color=['k', 'C2', 'C0']);
ax.hlines(0, -90, 90, color='k', linestyles='--', linewidth=0.5)
ax.set_ylim(-0.23, 1.1)
ax.set_ylabel('Temperature trend (°C per decade)'); ax.set_xlabel('Latitude');
plt.savefig('ERA5/era5_zonal_land_oceans.png', dpi=150, bbox_inches='tight');
with sns.axes_style('ticks'):
f, ax = plt.subplots()
dfs[['Land', 'Glaciers']].plot(ax=ax, color=['C0', 'C3']);
ax.hlines(0, -90, 90, color='k', linestyles='--', linewidth=0.5)
ax.set_ylim(-0.23, 1.1)
ax.set_ylabel('Temperature trend (°C per decade)'); ax.set_xlabel('Latitude');
plt.savefig('ERA5/era5_zonal_land_glaciers.png', dpi=150, bbox_inches='tight');
with sns.axes_style('ticks'):
f, ax = plt.subplots()
dfs[['Land', 'Glaciers', 'Above 1000']].plot(ax=ax, color=['C0', 'C3', 'k'], style=['-', '-', '--']);
ax.hlines(0, -90, 90, color='k', linestyles='--', linewidth=0.5)
ax.set_ylim(-0.23, 1.1)
ax.set_ylabel('Temperature trend (°C per decade)'); ax.set_xlabel('Latitude');
plt.savefig('ERA5/era5_zonal_land_glaciers_elev.png', dpi=150, bbox_inches='tight');
dfs['Cumulative All'] = df_trends['Weight All'][::-1].cumsum()
dfs['Cumulative Glaciers'] = df_trends['Weight Glaciers'][::-1].cumsum()
with sns.axes_style('ticks'):
f, ax = plt.subplots()
(dfs[['Weight All', 'Weight Glaciers', 'Cumulative All', 'Cumulative Glaciers']] * 100.).plot(ax=ax, secondary_y=['Cumulative All', 'Cumulative Glaciers'],
color=['C0', 'C3'], style=['-', '-', ':', ':']);
ax.set_ylabel('Percentage of total area (%)'); ax.set_xlabel('Latitude');
plt.savefig('ERA5/era5_zonal_weights.png', dpi=150, bbox_inches='tight');
dfs[['Cumulative All', 'Cumulative Glaciers']].loc[dfs['Cumulative Glaciers'] > 0.5].iloc[-1]
Cumulative All 0.945008 Cumulative Glaciers 0.500372 Name: 62.75, dtype: float64
import hvplot.pandas
import holoviews as hv
import geoviews as gv
import geoviews.feature as gf
ds['Trend'] = ds['trends'].astype(np.float32)
ds['StdErr'] = ds['trends_std_err'].astype(np.float32)
dataset = gv.Dataset(ds[['Trend', 'StdErr']], ['longitude', 'latitude'], vdims=[hv.Dimension('Trend', unit='°C per dec'),
hv.Dimension('StdErr', unit='°C per dec')],
crs=ccrs.PlateCarree())
print(dataset)
:Dataset [longitude,latitude] (Trend,StdErr)
from bokeh.models import HoverTool
hover = HoverTool(
tooltips=[
("Trend (°C per decade)", "@image{0[.]00} ± @StdErr{0[.]00}")
]
)
image = dataset.to(gv.Image).redim.range(Trend=(-1.5, 1.5))
plot = image.opts(cmap='RdBu_r', colorbar=True, width=700, height=350,
projection=ccrs.Robinson(), tools=[hover], active_tools=['pan','wheel_zoom']) * gf.coastline
renderer = hv.renderer('bokeh')
renderer.save(plot, 'ERA5/trend-map')