Area-weighted global temperature anomalies
Calculate standard and area-weighted global temperature average anomalies using NCEP/NCAR
imports¶
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
import xarray as xr
import numpy as np
import pandas as pd
from cartopy import crs as ccrs
In [3]:
from calendar import month_abbr
matplotlib parameters for the figures¶
In [4]:
plt.rcParams['figure.figsize'] = (10,8)
plt.rcParams['font.size'] = 14
plt.rcParams['lines.linewidth'] = 2
load the monthly NCEP (NCEP/NCAR) reanalysis temperature (2 m)¶
In [5]:
url = 'https://psl.noaa.gov/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface/air.mon.mean.nc'
In [6]:
ncep_temp = xr.open_dataset(url)
In [7]:
ncep_temp = ncep_temp.sortby('lat')
In [8]:
ncep_temp
Out[8]:
selects only complete years¶
In [9]:
ncep_temp = ncep_temp.sel(time=slice(None, '2021'))
levels for plotting¶
In [10]:
levels=np.arange(-30, 30+5, 5)
In [11]:
levels
Out[11]:
calculate the time-averaged temperature for the 1950s and the 2010s¶
In [12]:
ncep_dec_1 = ncep_temp.sel(time=slice('1950','1959')).mean('time')
ncep_dec_2 = ncep_temp.sel(time=slice('2010','2019')).mean('time')
parameters for the colorbar¶
In [13]:
cbar_kwargs = dict(orientation='horizontal', shrink=0.7, pad=0.05, label='\u00B0C')
map the average temperature and the difference between the 2 decades¶
In [14]:
f, axes = plt.subplots(nrows=3, figsize=(8, 14), subplot_kw=dict(projection=ccrs.PlateCarree(central_longitude=180)))
ncep_dec_1['air'].plot(ax=axes[0], transform=ccrs.PlateCarree(), levels=levels, cbar_kwargs=cbar_kwargs)
ncep_dec_2['air'].plot(ax=axes[1], transform=ccrs.PlateCarree(), levels=levels, cbar_kwargs=cbar_kwargs)
(ncep_dec_2 - ncep_dec_1)['air'].plot(ax=axes[2], transform=ccrs.PlateCarree(), levels=10, cbar_kwargs=cbar_kwargs)
[ax.set_title(title) for ax, title in zip(axes, ['1950s','2010s','difference'])]
[ax.coastlines() for ax in axes];
zoom in on the Northern Hemisphere high latitudes¶
In [15]:
cbar_kwargs.update(shrink=0.5)
In [16]:
f, axes = plt.subplots(nrows=3, figsize=(8, 14), subplot_kw=dict(projection=ccrs.NorthPolarStereo(central_longitude=0)))
ncep_dec_1['air'].plot(ax=axes[0], transform=ccrs.PlateCarree(), levels=levels, cbar_kwargs=cbar_kwargs)
ncep_dec_2['air'].plot(ax=axes[1], transform=ccrs.PlateCarree(), levels=levels, cbar_kwargs=cbar_kwargs)
(ncep_dec_2 - ncep_dec_1)['air'].plot(ax=axes[2], transform=ccrs.PlateCarree(), levels=10, cbar_kwargs=cbar_kwargs)
[ax.set_title(title) for ax, title in zip(axes, ['1950s','2010s','difference'])]
[ax.coastlines() for ax in axes];
[ax.set_extent([0, 360, 30, 90], crs=ccrs.PlateCarree()) for ax in axes];
what does the time-series of surface temperature looks like for one location¶
In [17]:
Auckland_coordinates = [174.7645, -36.8509] # longitude, latitude
In [18]:
Auckland = ncep_temp.sel(lon=Auckland_coordinates[0], lat=Auckland_coordinates[1], method='nearest')
In [19]:
Auckland['air'].plot()
Out[19]:
In [20]:
f, ax = plt.subplots()
Auckland.groupby(Auckland.time.dt.month).mean('time')['air'].plot(ax=ax, marker='o')
ax.set_xticks(np.arange(12) + 1)
ax.set_xticklabels(month_abbr[1:]);
ax.grid(ls=':')
ax.set_ylabel('\u00B0C')
Out[20]:
we calculate, for each grid point, the anomalies with respect to the mean seasonal cycle (average for each month)¶
1st step: Calculate the climatological 'normal' (now 1991 - 2010)¶
In [21]:
climatology = ncep_temp.sel(time=slice('1991','2010'))
In [22]:
climatology = climatology.groupby(climatology.time.dt.month).mean('time')
plots the climatologhy¶
In [23]:
cbar_kwargs = {'shrink':0.7, 'label':'\u00B0C'}
In [24]:
fg = climatology['air'].plot(figsize=(10,8), col='month', col_wrap=3, transform=ccrs.PlateCarree(), \
subplot_kws={"projection":ccrs.PlateCarree(central_longitude=180)}, cbar_kwargs=cbar_kwargs)
for i, ax in enumerate(fg.axes.flat):
ax.coastlines()
ax.gridlines()
ax.set_title(month_abbr[i+1])
2nd step: Remove this climatology from the raw data¶
In [25]:
ncep_temp_anomalies = ncep_temp.groupby(ncep_temp.time.dt.month) - climatology
plots the anomalies for the December 2021¶
In [26]:
ncep_temp_anomalies.isel(time=-1)['air'].plot()
Out[26]:
In [27]:
f, ax = plt.subplots(subplot_kw={"projection":ccrs.Orthographic(central_longitude=0, central_latitude=90)})
ncep_temp_anomalies.isel(time=-1)['air'].plot(ax=ax, levels=20, transform=ccrs.PlateCarree())
ax.coastlines(resolution='10m')
Out[27]:
give the large difference in the area represented by grid points at low and high latitudes, we will calculate the weighted average global mean temperature anomlaies and compare the results to the standard (unweighted) anomalies¶
In [28]:
latitudes = ncep_temp['lat']
In [29]:
latitudes
Out[29]:
In [30]:
weights = np.cos(np.deg2rad(latitudes))
weights.name = "weights"
weights
Out[30]:
In [31]:
weights.plot()
Out[31]:
Standard mean¶
In [32]:
standard_mean = ncep_temp_anomalies['air'].mean(['lat','lon'])
In [33]:
standard_mean_smooth = standard_mean.rolling(time=12*5, center=True).mean().dropna("time")
Weighted mean¶
In [34]:
weighted_mean = ncep_temp_anomalies['air'].weighted(weights)
In [35]:
weighted_mean = weighted_mean.mean(("lon", "lat"))
In [36]:
weighted_mean_smooth = weighted_mean.rolling(time=12*5, center=True).mean().dropna("time")
In [37]:
f, ax = plt.subplots()
weighted_mean.plot(ax=ax, label='weighted mean', color='coral', lw=0.5, alpha=0.5)
standard_mean.plot(ax=ax, label='standard mean', color='steelblue', lw=0.5, alpha=0.5)
weighted_mean_smooth.plot(ax=ax, label='weighted mean (5 years rolling ave.)', color='r')
standard_mean_smooth.plot(ax=ax, label='standard mean (5 years rolling ave.)', color='b')
ax.legend(fontsize=12, frameon=False)
ax.grid(ls=':')
ax.axhline(0, color='k', zorder=0)
ax.set_title('NCEP/NCAR global temperature average')
ax.set_ylabel('\u00B0C' + ' anomaly')
Out[37]:
exporting the time-series in csv¶
In [38]:
standard_mean
Out[38]:
In [39]:
weighted_mean
Out[39]:
In [40]:
weighted_mean.name = 'weighted mean'
In [41]:
standard_mean.name = 'standard mean'
In [42]:
global_means = xr.merge([weighted_mean, standard_mean])
In [43]:
global_means = global_means.to_pandas()
In [44]:
global_means.head()
Out[44]:
In [45]:
global_means = global_means.drop('month', axis=1)
In [46]:
global_means.head()
Out[46]:
In [47]:
global_means.plot()
Out[47]:
In [48]:
global_means_smooth = global_means.rolling(window=12*5, min_periods=12*5, center=True).mean()
In [49]:
global_means_smooth.columns = ['weighted mean (5 years rolling ave.)', \
'standard mean (5 years rolling ave.)']
In [50]:
global_means_smooth.plot()
Out[50]:
In [51]:
all_means = pd.concat([global_means, global_means_smooth], axis=1)
In [52]:
all_means
Out[52]:
In [53]:
all_means.to_csv('../data/NCEP_global_means.csv')
In [54]:
data = pd.read_csv('../data/NCEP_global_means.csv', index_col=0, parse_dates=True)
In [55]:
data.head()
Out[55]:
In [56]:
data.plot()
Out[56]:
Comments