California august temperature anomaly¶
How anomalous was the August 2020 average temperature?
In this first section, we load required packages and modules
[1]:
##This is so variables get printed within jupyter
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
[2]:
##import packages
import os
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import matplotlib.ticker as mticker
import cdsapi
#for rank calculation
# import bottleneck
[3]:
os.chdir(os.path.abspath('../../'))
Load ERA5¶
We have retrieved netcdf files of global monthly 2m temperature and 2m dewpoint temperature for each year over 1979-2020.
We load all files with xarray open_mfdataset.
[4]:
ERA5 = xr.open_mfdataset('../California_example/ERA5/ERA5_????.nc',combine='by_coords') ## open the data
ERA5#
[4]:
xarray.Dataset
- latitude: 51
- longitude: 61
- time: 42
- longitude(longitude)float32-130.0 -129.0 ... -71.0 -70.0
- units :
- degrees_east
- long_name :
- longitude
array([-130., -129., -128., -127., -126., -125., -124., -123., -122., -121., -120., -119., -118., -117., -116., -115., -114., -113., -112., -111., -110., -109., -108., -107., -106., -105., -104., -103., -102., -101., -100., -99., -98., -97., -96., -95., -94., -93., -92., -91., -90., -89., -88., -87., -86., -85., -84., -83., -82., -81., -80., -79., -78., -77., -76., -75., -74., -73., -72., -71., -70.], dtype=float32) - latitude(latitude)float3270.0 69.0 68.0 ... 22.0 21.0 20.0
- units :
- degrees_north
- long_name :
- latitude
array([70., 69., 68., 67., 66., 65., 64., 63., 62., 61., 60., 59., 58., 57., 56., 55., 54., 53., 52., 51., 50., 49., 48., 47., 46., 45., 44., 43., 42., 41., 40., 39., 38., 37., 36., 35., 34., 33., 32., 31., 30., 29., 28., 27., 26., 25., 24., 23., 22., 21., 20.], dtype=float32) - time(time)datetime64[ns]1979-08-01 ... 2020-08-01
- long_name :
- time
array(['1979-08-01T00:00:00.000000000', '1980-08-01T00:00:00.000000000', '1981-08-01T00:00:00.000000000', '1982-08-01T00:00:00.000000000', '1983-08-01T00:00:00.000000000', '1984-08-01T00:00:00.000000000', '1985-08-01T00:00:00.000000000', '1986-08-01T00:00:00.000000000', '1987-08-01T00:00:00.000000000', '1988-08-01T00:00:00.000000000', '1989-08-01T00:00:00.000000000', '1990-08-01T00:00:00.000000000', '1991-08-01T00:00:00.000000000', '1992-08-01T00:00:00.000000000', '1993-08-01T00:00:00.000000000', '1994-08-01T00:00:00.000000000', '1995-08-01T00:00:00.000000000', '1996-08-01T00:00:00.000000000', '1997-08-01T00:00:00.000000000', '1998-08-01T00:00:00.000000000', '1999-08-01T00:00:00.000000000', '2000-08-01T00:00:00.000000000', '2001-08-01T00:00:00.000000000', '2002-08-01T00:00:00.000000000', '2003-08-01T00:00:00.000000000', '2004-08-01T00:00:00.000000000', '2005-08-01T00:00:00.000000000', '2006-08-01T00:00:00.000000000', '2007-08-01T00:00:00.000000000', '2008-08-01T00:00:00.000000000', '2009-08-01T00:00:00.000000000', '2010-08-01T00:00:00.000000000', '2011-08-01T00:00:00.000000000', '2012-08-01T00:00:00.000000000', '2013-08-01T00:00:00.000000000', '2014-08-01T00:00:00.000000000', '2015-08-01T00:00:00.000000000', '2016-08-01T00:00:00.000000000', '2017-08-01T00:00:00.000000000', '2018-08-01T00:00:00.000000000', '2019-08-01T00:00:00.000000000', '2020-08-01T00:00:00.000000000'], dtype='datetime64[ns]')
- t2m(time, latitude, longitude)float32dask.array<chunksize=(1, 51, 61), meta=np.ndarray>
- units :
- K
- long_name :
- 2 metre temperature
Array Chunk Bytes 522.65 kB 12.44 kB Shape (42, 51, 61) (1, 51, 61) Count 126 Tasks 42 Chunks Type float32 numpy.ndarray - d2m(time, latitude, longitude)float32dask.array<chunksize=(1, 51, 61), meta=np.ndarray>
- units :
- K
- long_name :
- 2 metre dewpoint temperature
Array Chunk Bytes 522.65 kB 12.44 kB Shape (42, 51, 61) (1, 51, 61) Count 126 Tasks 42 Chunks Type float32 numpy.ndarray
- Conventions :
- CF-1.6
- history :
- 2021-02-09 15:47:30 GMT by grib_to_netcdf-2.16.0: /opt/ecmwf/eccodes/bin/grib_to_netcdf -S param -o /cache/data7/adaptor.mars.internal-1612885647.5164566-26296-16-d75758d7-9566-4339-b834-0e54e545a802.nc /cache/tmp/d75758d7-9566-4339-b834-0e54e545a802-adaptor.mars.internal-1612885647.5171535-26296-6-tmp.grib
Retrieve the land-sea mask¶
We retrieve the land-sea mask for ERA5 from CDS.
[5]:
c = cdsapi.Client()
c.retrieve(
'reanalysis-era5-single-levels-monthly-means',
{
'format': 'netcdf',
'product_type': 'monthly_averaged_reanalysis',
'variable': 'land_sea_mask',
'grid': [1.0, 1.0],
'year': '1979',
'month': '01',
'time': '00:00',
},
'../California_example/ERA_landsea_mask.nc')
2021-08-12 10:14:27,929 INFO Welcome to the CDS
2021-08-12 10:14:27,931 INFO Sending request to https://cds.climate.copernicus.eu/api/v2/resources/reanalysis-era5-single-levels-monthly-means
2021-08-12 10:14:28,432 INFO Request is completed
2021-08-12 10:14:28,433 INFO Downloading https://download-0011.copernicus-climate.eu/cache-compute-0011/cache/data5/adaptor.mars.internal-1628605177.1251855-440-12-4553cc0b-4a6e-4405-85fa-b9d5396ac1f1.nc to ../California_example/ERA_landsea_mask.nc (130.5K)
2021-08-12 10:14:28,662 INFO Download rate 574.8K/s
[5]:
Result(content_length=133596,content_type=application/x-netcdf,location=https://download-0011.copernicus-climate.eu/cache-compute-0011/cache/data5/adaptor.mars.internal-1628605177.1251855-440-12-4553cc0b-4a6e-4405-85fa-b9d5396ac1f1.nc)
And here we open the dataset. It contains dimensionless values from 0 to 1. From CDS: “Grid boxes where this parameter has a value above 0.5 can be comprised of a mixture of land and inland water but not ocean. Grid boxes with a value of 0.5 and below can only be comprised of a water surface.”
[5]:
LSMask = xr.open_dataset('../California_example/ERA_landsea_mask.nc')
LSMask.load()
[5]:
xarray.Dataset
- latitude: 181
- longitude: 360
- time: 1
- longitude(longitude)float320.0 1.0 2.0 ... 357.0 358.0 359.0
- units :
- degrees_east
- long_name :
- longitude
array([ 0., 1., 2., ..., 357., 358., 359.], dtype=float32)
- latitude(latitude)float3290.0 89.0 88.0 ... -89.0 -90.0
- units :
- degrees_north
- long_name :
- latitude
array([ 90., 89., 88., 87., 86., 85., 84., 83., 82., 81., 80., 79., 78., 77., 76., 75., 74., 73., 72., 71., 70., 69., 68., 67., 66., 65., 64., 63., 62., 61., 60., 59., 58., 57., 56., 55., 54., 53., 52., 51., 50., 49., 48., 47., 46., 45., 44., 43., 42., 41., 40., 39., 38., 37., 36., 35., 34., 33., 32., 31., 30., 29., 28., 27., 26., 25., 24., 23., 22., 21., 20., 19., 18., 17., 16., 15., 14., 13., 12., 11., 10., 9., 8., 7., 6., 5., 4., 3., 2., 1., 0., -1., -2., -3., -4., -5., -6., -7., -8., -9., -10., -11., -12., -13., -14., -15., -16., -17., -18., -19., -20., -21., -22., -23., -24., -25., -26., -27., -28., -29., -30., -31., -32., -33., -34., -35., -36., -37., -38., -39., -40., -41., -42., -43., -44., -45., -46., -47., -48., -49., -50., -51., -52., -53., -54., -55., -56., -57., -58., -59., -60., -61., -62., -63., -64., -65., -66., -67., -68., -69., -70., -71., -72., -73., -74., -75., -76., -77., -78., -79., -80., -81., -82., -83., -84., -85., -86., -87., -88., -89., -90.], dtype=float32) - time(time)datetime64[ns]1979-01-01
- long_name :
- time
array(['1979-01-01T00:00:00.000000000'], dtype='datetime64[ns]')
- lsm(time, latitude, longitude)float320.0 0.0 0.0 0.0 ... 1.0 1.0 1.0 1.0
- units :
- (0 - 1)
- long_name :
- Land-sea mask
- standard_name :
- land_binary_mask
array([[[0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], [0., 0., 0., ..., 0., 0., 0.], ..., [1., 1., 1., ..., 1., 1., 1.], [1., 1., 1., ..., 1., 1., 1.], [1., 1., 1., ..., 1., 1., 1.]]], dtype=float32)
- Conventions :
- CF-1.6
- history :
- 2021-08-10 14:19:37 GMT by grib_to_netcdf-2.20.0: /opt/ecmwf/mars-client/bin/grib_to_netcdf -S param -o /cache/data5/adaptor.mars.internal-1628605177.1251855-440-12-4553cc0b-4a6e-4405-85fa-b9d5396ac1f1.nc /cache/tmp/4553cc0b-4a6e-4405-85fa-b9d5396ac1f1-adaptor.mars.internal-1628605176.529759-440-21-tmp.grib
[6]:
LSMask['lsm'].plot()
[6]:
<matplotlib.collections.QuadMesh at 0x7fbcb9f389a0>
We can select all gridcelss where the land-sea mask values are > 0.5 in order to remove ocean gridcells:
[7]:
LSMask['lsm'].where(LSMask['lsm'] > 0.5).plot()
[7]:
<matplotlib.collections.QuadMesh at 0x7fbcb9d95d00>
The longitude values in this dataset run from 0:360. We need to convert this to -180:180:
[8]:
# convert the longitude from 0:360 to -180:180
LSMask['longitude'] = (((LSMask['longitude'] + 180) % 360) - 180)
Calculating the anomaly¶
We want to show how anomalous the recorded monthly average temperature for 2020 is compared to the 1979-2010 average. We first calculate the temperature anomaly from the 1979-2010 mean and then calculate the standardardized anomaly by dividing the anomaly by the standard deviation:
[9]:
ERA5_anomaly = ERA5['t2m'] - ERA5['t2m'].sel(time=slice('1979','2010')).mean('time')
ERA5_anomaly.attrs = {
'long_name': 'August temperature anomaly',
'units': 'C'
}
ERA5_sd_anomaly = ERA5_anomaly / ERA5['t2m'].sel(time=slice('1979','2010')).std('time')
ERA5_sd_anomaly.attrs = {
'long_name': 'August temperature standardized anomaly',
'units': '-'
}
Plotting the 2020 temperature anomaly¶
We define a function to plot the data on a global map:
[10]:
def plot_California(ERA5_input):
extent = [-120, -80, 20, 50]
central_lon = np.mean(extent[:2])
central_lat = np.mean(extent[2:])
plt.figure(figsize=(12, 6))
ax = plt.axes(projection=ccrs.AlbersEqualArea(central_lon, central_lat))
ax.set_extent(extent)
ERA5_input.plot(
ax=ax,
transform=ccrs.PlateCarree(),
extend='both')
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.coastlines(
resolution='110m') #Currently can be one of “110m”, “50m”, and “10m”.
ax.set_title('')
gl = ax.gridlines(crs=ccrs.PlateCarree(),
draw_labels=True,
linewidth=1,
color='gray',
alpha=0.5,
linestyle='--')
gl.top_labels = False
gl.right_labels = False
And plot the ERA5 2020 temperature anomaly
[11]:
plot_California(ERA5_anomaly.sel(time = '2020'))
# plt.savefig('graphs/California_anomaly.png')
Plot the standardized anomaly
[12]:
plot_California(ERA5_sd_anomaly.sel(time = '2020'))
Selecting a domain for further analysis¶
We define the domain as a contiguous, land-only region within the box 125:100W, 20:45N with temperature anomalies above 2 standard deviation.
[17]:
ERA5_sd_anomaly_masked = (ERA5_sd_anomaly.
sel(longitude = slice(-125,-100), #select the domain
latitude = slice(45,20)). #select the domain
where(ERA5_sd_anomaly.sel(time = '2020').squeeze('time')>2). #Select the region where sd>2
where(LSMask['lsm'].sel(time = '1979').squeeze('time') > 0.5) #Select land-only gridcells
)
ERA5_sd_anomaly_masked.load()
[17]:
xarray.DataArray
't2m'
- time: 42
- latitude: 26
- longitude: 26
- nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan
array([[[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], [[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], [[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], ..., [[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], [[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]], [[nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], ..., [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan], [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32) - latitude(latitude)float6445.0 44.0 43.0 ... 22.0 21.0 20.0
- units :
- degrees_north
- long_name :
- latitude
array([45., 44., 43., 42., 41., 40., 39., 38., 37., 36., 35., 34., 33., 32., 31., 30., 29., 28., 27., 26., 25., 24., 23., 22., 21., 20.]) - longitude(longitude)float64-125.0 -124.0 ... -101.0 -100.0
- units :
- degrees_east
- long_name :
- longitude
array([-125., -124., -123., -122., -121., -120., -119., -118., -117., -116., -115., -114., -113., -112., -111., -110., -109., -108., -107., -106., -105., -104., -103., -102., -101., -100.]) - time(time)datetime64[ns]1979-08-01 ... 2020-08-01
- long_name :
- time
array(['1979-08-01T00:00:00.000000000', '1980-08-01T00:00:00.000000000', '1981-08-01T00:00:00.000000000', '1982-08-01T00:00:00.000000000', '1983-08-01T00:00:00.000000000', '1984-08-01T00:00:00.000000000', '1985-08-01T00:00:00.000000000', '1986-08-01T00:00:00.000000000', '1987-08-01T00:00:00.000000000', '1988-08-01T00:00:00.000000000', '1989-08-01T00:00:00.000000000', '1990-08-01T00:00:00.000000000', '1991-08-01T00:00:00.000000000', '1992-08-01T00:00:00.000000000', '1993-08-01T00:00:00.000000000', '1994-08-01T00:00:00.000000000', '1995-08-01T00:00:00.000000000', '1996-08-01T00:00:00.000000000', '1997-08-01T00:00:00.000000000', '1998-08-01T00:00:00.000000000', '1999-08-01T00:00:00.000000000', '2000-08-01T00:00:00.000000000', '2001-08-01T00:00:00.000000000', '2002-08-01T00:00:00.000000000', '2003-08-01T00:00:00.000000000', '2004-08-01T00:00:00.000000000', '2005-08-01T00:00:00.000000000', '2006-08-01T00:00:00.000000000', '2007-08-01T00:00:00.000000000', '2008-08-01T00:00:00.000000000', '2009-08-01T00:00:00.000000000', '2010-08-01T00:00:00.000000000', '2011-08-01T00:00:00.000000000', '2012-08-01T00:00:00.000000000', '2013-08-01T00:00:00.000000000', '2014-08-01T00:00:00.000000000', '2015-08-01T00:00:00.000000000', '2016-08-01T00:00:00.000000000', '2017-08-01T00:00:00.000000000', '2018-08-01T00:00:00.000000000', '2019-08-01T00:00:00.000000000', '2020-08-01T00:00:00.000000000'], dtype='datetime64[ns]')
- long_name :
- August temperature standardized anomaly
- units :
- -
The resulting domain looks as follows:
[18]:
plot_California(ERA5_sd_anomaly_masked.sel(time = '2020'))
Let’s make a nicer plot, using the previously defined function but adding an outline for the domain
[21]:
def plot_California(ERA5_input, ERA5_masked):
extent = [-120, -80, 20, 50]
central_lon = np.mean(extent[:2])
central_lat = np.mean(extent[2:])
plt.figure(frameon=False, figsize=(90 / 25.4, 60 / 25.4))
ax = plt.axes(projection=ccrs.AlbersEqualArea(central_lon, central_lat))
ax.set_extent(extent)
ERA5_input.plot(
ax=ax,
transform=ccrs.PlateCarree(),
extend='both')
(ERA5_masked.fillna(0).sel(time = '2020').
squeeze('time').
plot.contour(levels = [0],
colors = 'black',
transform=ccrs.PlateCarree(),
ax = ax)
)
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.coastlines(
resolution='110m') #Currently can be one of “110m”, “50m”, and “10m”.
ax.set_title('')
gl = ax.gridlines(crs=ccrs.PlateCarree(),
draw_labels=True,
linewidth=1,
color='gray',
alpha=0.5,
linestyle='--')
gl.top_labels = False
gl.right_labels = False
[20]:
plot_California(ERA5_sd_anomaly.sel(time = '2020'),ERA5_sd_anomaly_masked)
/soge-home/users/cenv0732/.conda/envs/UNSEEN-open/lib/python3.8/site-packages/cartopy/mpl/geoaxes.py:1478: UserWarning: No contour levels were found within the data range.
result = matplotlib.axes.Axes.contour(self, *args, **kwargs)
Some functions that were used to create the figure for publication: - Set figure size as 90 by 60mm and remove frame: plt.figure(frameon=False, figsize=(90 / 25.4, 60 / 25.4)) - Set the font: plt.rcParams["font.family"] = "sans-serif" - Set the font size: plt.rcParams['font.size'] = 10 - Set the font type so editing software (e.g. inkscape) can recognize text for svg graphics: plt.rcParams['svg.fonttype'] = 'none' - Same but for pdf files: plt.rcParams['pdf.fonttype'] = 42
[22]:
plt.rcParams["font.family"] = "sans-serif" ##change font
plt.rcParams['font.size'] = 10 ## change font size
# plt.rcParams['svg.fonttype'] = 'none' ## so inkscape recognized texts in svg file
plt.rcParams['pdf.fonttype'] = 42 ## so illustrator can recognize text
plot_California(ERA5_sd_anomaly.sel(time = '2020'),ERA5_sd_anomaly_masked)
plt.savefig('graphs/California_sd_anomaly_contour.pdf')
/soge-home/users/cenv0732/.conda/envs/UNSEEN-open/lib/python3.8/site-packages/cartopy/mpl/geoaxes.py:1478: UserWarning: No contour levels were found within the data range.
result = matplotlib.axes.Axes.contour(self, *args, **kwargs)
Timeseries for the selected domain¶
Here we take the areal average over the selected domain and plot the resulting timeseries. Gridcell need to be weighted with their cell area when taking the spatial mean. We first calculate the area weights and then use these to average.
[23]:
area_weights = np.cos(np.deg2rad(ERA5_sd_anomaly.latitude))
ERA5_std_anomaly_timeseries = ERA5_sd_anomaly_masked.weighted(area_weights).mean(['longitude','latitude'])
[28]:
plt.figure(frameon=False, figsize=(90 / 25.4, 60 / 25.4))
ERA5_std_anomaly_timeseries.plot()
plt.ylabel('August standardized temperature anomaly')
plt.savefig('graphs/California_anomaly_timeseries.pdf')
[28]:
<Figure size 255.118x170.079 with 0 Axes>
[28]:
[<matplotlib.lines.Line2D at 0x7fbcb0643f40>]
[28]:
Text(0, 0.5, 'August standardized temperature anomaly')
Another option would be to use the Californian domain using regionmask: states = regionmask.defined_regions.natural_earth.us_states_50