Table of Contents
- 1 Create and Connect to Dask Distributed Cluster
- 2 Load data into xarray from an intake catalog
- 3 Compute Weighted Means
- 4 Get Observations (HadCRUT4; Morice et al. 2012)
- 5 Figure 2: Global surface temperature anomaly (1961-90 base period) for individual ensemble members, and observations
- 6 Compute Linear Trend for Winter Seasons
- 7 Figure 4: Global maps of historical (1979 - 2012) boreal winter (DJF) surface air trends
In [1]:
# TODO: document this notebook
# References: https://journals.ametsoc.org/doi/full/10.1175/BAMS-D-13-00255.1
In [2]:
%matplotlib inline
import warnings
warnings.filterwarnings("ignore")
# Silence dask.distributed logs
import logging
logger = logging.getLogger("distributed.utils_perf")
logger.setLevel(logging.ERROR)
import intake
import numpy as np
import pandas as pd
import xarray as xr
Create and Connect to Dask Distributed Cluster¶
In [4]:
from dask_kubernetes import KubeCluster
from dask.distributed import Client
cluster = KubeCluster()
cluster.adapt(minimum=2, maximum=100, wait_count=60)
client = Client(cluster)
cluster
☝️ Don't forget to click the link above to view the scheduler dashboard!
Load data into xarray from an intake catalog¶
In [5]:
cat = intake.Catalog(
"https://raw.githubusercontent.com/NCAR/cesm-lens-aws/master/intake-catalogs/atmosphere/daily.yaml"
)
ds_20C = cat["reference_height_temperature_20C"].to_dask()
ds_rcp = cat["reference_height_temperature_RCP85"].to_dask()
t_20c = ds_20C["TREFHT"]
t_rcp = ds_rcp["TREFHT"]
In [6]:
t_ref = t_20c.sel(time=slice("1961", "1990"))
t_ref
Out[6]:
In [7]:
areacella = ds_20C.area
total_area = areacella.sum()
areacella
Out[7]:
Compute Weighted Means¶
In [8]:
t_ref_ts = (
(t_ref.resample(time="AS").mean("time") * areacella).sum(dim=("lat", "lon"))
/ total_area
).mean(dim=("time", "member_id"))
t_20c_ts = (
(t_20c.resample(time="AS").mean("time") * areacella).sum(dim=("lat", "lon"))
) / total_area
t_rcp_ts = (
(t_rcp.resample(time="AS").mean("time") * areacella).sum(dim=("lat", "lon"))
) / total_area
In [9]:
%%time
t_ref_mean = t_ref_ts.load()
t_ref_mean
Out[9]:
In [10]:
t_20c_ts_df = t_20c_ts.to_series().unstack().T
t_20c_ts_df.head()
Out[10]:
In [11]:
t_rcp_ts_df = t_rcp_ts.to_series().unstack().T
t_rcp_ts_df.head()
Out[11]:
Get Observations (HadCRUT4; Morice et al. 2012)¶
In [12]:
ds = xr.open_dataset(
"https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/cru/hadcrut4/air.mon.anom.median.nc"
).load()
ds
Out[12]:
- Obs mean: weight by days in each month
In [13]:
def weighted_temporal_mean(ds):
time_bound_diff = ds.time_bnds.diff(dim="nbnds")[:, 0]
wgts = time_bound_diff.groupby("time.year") / time_bound_diff.groupby(
"time.year"
).sum(xr.ALL_DIMS)
np.testing.assert_allclose(wgts.groupby("time.year").sum(xr.ALL_DIMS), 1.0)
obs = ds["air"]
cond = obs.isnull()
ones = xr.where(cond, 0.0, 1.0)
obs_sum = (obs * wgts).resample(time="AS").sum(dim="time")
ones_out = (ones * wgts).resample(time="AS").sum(dim="time")
obs_s = (obs_sum / ones_out).mean(("lat", "lon")).to_series()
return obs_s
In [14]:
obs_s = weighted_temporal_mean(ds)
obs_s.head()
Out[14]:
In [15]:
all_ts_anom = pd.concat([t_20c_ts_df, t_rcp_ts_df]) - t_ref_mean.data
years = [val.year for val in all_ts_anom.index]
- Confirm that after using area weighted average, max temp increase is 5k
In [16]:
np.testing.assert_allclose(all_ts_anom.values.max(), 5.0, rtol=0.02)
Figure 2: Global surface temperature anomaly (1961-90 base period) for individual ensemble members, and observations¶
In [17]:
import matplotlib.pyplot as plt
In [18]:
ax = plt.axes()
ax.tick_params(right=True, top=True, direction="out", length=6, width=2, grid_alpha=0.5)
ax.plot(years, all_ts_anom, color="grey")
ax.plot(years, all_ts_anom[1], color="black")
ax.plot(obs_s.index.year.tolist(), obs_s, color="red")
ax.text(
0.3,
0.4,
"observations",
verticalalignment="bottom",
horizontalalignment="left",
transform=ax.transAxes,
color="red",
fontsize=10,
)
ax.text(
0.3,
0.3,
"members 1-40",
verticalalignment="bottom",
horizontalalignment="left",
transform=ax.transAxes,
color="grey",
fontsize=10,
)
ax.set_xticks([1850, 1920, 1950, 2000, 2050, 2100])
plt.ylim(-1, 5)
plt.xlim(1850, 2100)
plt.ylabel("Global Surface\nTemperature Anomaly (K)")
plt.show()
Compute Linear Trend for Winter Seasons¶
In [19]:
def linear_trend(da, dim="time"):
da_chunk = da.chunk({dim: -1})
trend = xr.apply_ufunc(
calc_slope,
da_chunk,
vectorize=True,
input_core_dims=[[dim]],
output_core_dims=[[]],
output_dtypes=[np.float],
dask="parallelized",
)
return trend
def calc_slope(y):
"""ufunc to be used by linear_trend"""
x = np.arange(len(y))
return np.polyfit(x, y, 1)[0]
In [20]:
t = xr.concat([t_20c, t_rcp], dim="time")
seasons = t.sel(time=slice("1979", "2012")).resample(time="QS-DEC").mean("time")
# Include only full seasons from 1979 and 2012
seasons = seasons.sel(time=slice("1979", "2012")).load()
seasons
Out[20]:
In [21]:
winter_seasons = seasons.sel(
time=seasons.time.where(seasons.time.dt.month == 12, drop=True)
)
winter_trends = linear_trend(
winter_seasons.chunk({"lat": 20, "lon": 20, "time": -1})
).load() * len(winter_seasons.time)
In [22]:
# Make sure that we have 34 seasons
assert len(winter_seasons.time) == 34
Figure 4: Global maps of historical (1979 - 2012) boreal winter (DJF) surface air trends¶
In [23]:
!pip install git+https://github.com/hhuangwx/cmaps.git -q
In [24]:
import cmaps # for NCL colormaps
import cartopy.crs as ccrs
import matplotlib.colors as colors
cmap = cmaps.amwg_blueyellowred
levels = [-7, -6, -5, -4, -3, -2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, 2, 3, 4, 5, 6, 7]
norm = colors.Normalize(vmin=levels[0], vmax=levels[-1], clip=True)
In [25]:
fig = plt.figure(dpi=300)
fg = winter_trends.plot(
col="member_id",
col_wrap=4,
transform=ccrs.PlateCarree(),
subplot_kws={"projection": ccrs.Robinson(central_longitude=180)},
add_colorbar=False,
levels=levels,
cmap=cmap,
norm=norm,
extend="neither",
)
for ax in fg.axes.flat:
ax.coastlines(color="grey")
# TODO: move the subplot title to lower left corners
# TODO: Add obs panel and ensemble mean at the end
fg.add_colorbar(orientation="horizontal")
fg.cbar.set_label("1979-2012 DJF surface air temperature trends (K/34 years)")
fg.cbar.set_ticks(levels)
fg.cbar.set_ticklabels(levels)
In [ ]:
# Gracefully destroy/close our cluster
client.close()
cluster.close()