How to make anomalies and plot timeseries (Example: Nino 3.4 Index), and take Composite mean of anomalies for Elnino and Lanina
Authors:
Nazimur Rashid Chowdhury (rashidnazimur@gmail.com) Abdullah Al Fahad (a.fahad@nasa.gov)
Lets load the data first.¶
In [24]:
# Run this cell when you are in colab and then restart runtime
# Do not worry about this cell, you can always copy paste
!apt-get install libproj-dev proj-data proj-bin
!apt-get install libgeos-dev
!pip install cython
!pip install cartopy
!apt-get -qq install python-cartopy python3-cartopy
!pip uninstall -y shapely # cartopy and shapely aren't friends (early 2020)
!pip install shapely --no-binary shapely
In [25]:
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
In [26]:
# from google.colab import drive
# drive.mount('/content/drive')
Now we will download Reynold's SST using WGET.¶
In [131]:
!wget https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2/sst.mnmean.nc
In [132]:
# Read the data using xarray
ds = xr.open_dataset('/content/sst.mnmean.nc')
#ds = xr.decode_cf(ds, use_cftime=True) # Ignore this line. It is only needed if time coordinates are given in floating data and you need to convert it to cftime format
ds
Out[132]:
In [133]:
# Choose the SST variable store it's average
tos = ds.sst
average_sst = tos.mean(dim='time')
del ds
In [134]:
fig, (ax) = plt.subplots(
ncols=1, nrows=1, figsize=[8, 4], subplot_kw={"projection": ccrs.Robinson()}
)
# plot the data
average_sst.plot(
ax=ax,
x="lon",
y="lat",
transform=ccrs.PlateCarree(),
vmin=-1,
vmax=30,
cmap="coolwarm",
robust=True,
)
ax.coastlines()
ax.set_global()
ax.set_title("Average Sea Surface Temperature (°C)")
Out[134]:
Computing anomalies¶
In this tutorial we will learn about how to deal with time series data. In order to do that, we have chosen a phenomenon called El Nino Southern Oscillation (ENSO). It is a climate pattern in the pacific ocean and can affect weather pattern worldwide. Have a go through this link: https://oceanservice.noaa.gov/facts/ninonina.html
In [135]:
# data has been grouped by month
gb = tos.groupby('time.month')
# calculating anomaly
tos_anom = gb - gb.mean(dim='time')
In [136]:
#ploting 2016 january anomalies which is an elnino year.
tos_anom_at_time = tos_anom.sel(time='2016-01-01')
print(tos_anom_at_time.min())
print(tos_anom_at_time.max())
In [137]:
# Select the temperature anomaly at the chosen time
tos_anom_at_time = tos_anom.sel(time='2016-01-01')
lon = tos_anom_at_time.lon
lat = tos_anom_at_time.lat
fig, (ax) = plt.subplots(
ncols=1, nrows=1, figsize=[8, 4], subplot_kw={"projection": ccrs.Robinson()}
)
# plot the data
tos_anom_at_time.plot(
ax=ax,
x="lon",
y="lat",
transform=ccrs.PlateCarree(),
vmin=-4,
vmax=4,
cmap="coolwarm",
robust=True,
)
ax.coastlines()
ax.set_global()
ax.set_title('Sea Surface Temperature Anomaly at Time: {}'.format(tos_anom_at_time.time.values), pad=20)
Out[137]:
In [138]:
# Now we will calculate the mean global anomaly by averaging over all the lats and lons
unweighted_mean_global_anom = tos_anom.mean(dim=['lat','lon'])
unweighted_mean_global_anom.plot()
Out[138]:
As we move toward the poles, each grid cell covers less area. So we need to account for that when calculating mean global anomaly. This is called weighted average and we will use the grid cell area as weights. If the dataset doesn't include any variable containing area of the cells, you can apply the method below which is using the latitudes to calculate the area of the grid cells.¶
In [139]:
# For a rectangular grid the cosine of the latitude is proportional to the grid cell area.
weights = np.cos(np.deg2rad(tos_anom.lat))
weights.name = "weights"
In [140]:
# Calculating the weighted global anomaly using the 'weights' we have created
weighted_mean_global_anom = tos_anom.weighted(weights).mean(dim=['lat','lon'])
In [141]:
# Plotting weighted and unweighted mean together to see difference
unweighted_mean_global_anom.plot(size=5)
weighted_mean_global_anom.plot()
plt.legend(['unweighted','weighted'])
Out[141]:
Now we will calculate rolling mean. A rolling mean helps to reduce the impact of short-term fluctuations or noise in the data. It provides a clearer representation of the underlying trend or pattern by averaging out these variations over a specified window.¶
In [142]:
# Using rolling method to compute averages in a moving window of 5-months of data:
m_avg = weighted_mean_global_anom.rolling(time=5, center=True).mean()
In [143]:
# Now using the methods we will calculate and plot Nino 3.4 Index
# First we will subset our data (for Nino 3.4 region) using two methods
# Slice
tos_nino34 = tos.sel(lat=slice(-5,5), lon=slice(190,240))
# Masking with 'where'
tos_nino34 = tos.where((tos.lat>-5) & (tos.lat<5) & (tos.lon>190) & (tos.lon<240), drop=True)
In [144]:
# Here we are repeating the computation we learned above for calculating anomalies. This time for the Nino 3.4 area.
gb2 = tos_nino34.groupby('time.month')
nino34_anom = gb2 - gb2.mean(dim='time')
index_nino34 = nino34_anom.mean(dim=['lat','lon'])
index_nino34_rolling_mean = index_nino34.rolling(time=5, center=True).mean()
In [145]:
# Let's plot mean anomaly and 5-month running mean anomaly to compare
index_nino34.plot(size=6)
index_nino34_rolling_mean.plot()
plt.legend(['anomaly', '5-month running mean anomaly'])
plt.title('SST anomaly over the Niño 3.4 region');
In [145]:
In [145]:
In [146]:
## routines for calculating seasonal climatology and seasonal timeseries
import xarray as xr
import numpy as np
dpm = {'noleap': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
'365_day': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
'standard': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
'gregorian': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
'proleptic_gregorian': [0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
'all_leap': [0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
'366_day': [0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31],
'360_day': [0, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30]}
def leap_year(year, calendar='standard'):
"""Determine if year is a leap year
Args:
year (numeric)
"""
leap = False
if ((calendar in ['standard', 'gregorian',
'proleptic_gregorian', 'julian']) and
(year % 4 == 0)):
leap = True
if ((calendar == 'proleptic_gregorian') and
(year % 100 == 0) and
(year % 400 != 0)):
leap = False
elif ((calendar in ['standard', 'gregorian']) and
(year % 100 == 0) and (year % 400 != 0) and
(year < 1583)):
leap = False
return leap
def get_days_per_mon(time, calendar='standard'):
"""
return a array of days per month corresponding to the months provided in `months`
Args: time (CFTimeIndex): ie. ds.time.to_index()
calendar (str): default 'standard'
"""
month_length = np.zeros(len(time), dtype=np.int)
cal_days = dpm[calendar]
for i, (month, year) in enumerate(zip(time.month, time.year)):
month_length[i] = cal_days[month]
if leap_year(year, calendar=calendar):
month_length[i] += 1
return month_length
def season_mean(ds, season = "all", cal = "none"):
""" calculate climatological mean by season
Args: ds (xarray.Dataset): dataset
var (str): variable to use
season (str): "all", 'DJF', "MAM", "JJA", "SON"
cal (str): "none"(default) or calendar used for weighting months by number of days
"""
## no weighting of months:
if cal == "none":
if season == "all":
## calculate mean for all season
smean = ds.groupby('time.season').mean('time')
else :
## calculate mean for specified season
smean = ds.where(ds['time.season'] == season).mean('time')
return smean
## weighted months
else:
## create array of month_length (number of days in each month)
## assign time coords matching original ds
month_length = xr.DataArray(get_days_per_mon(ds.time.to_index(), calendar=cal),
coords=[ds.time], name='month_length')
## Calculate the weights by grouping by 'time.season'
weights = month_length.groupby('time.season') / month_length.groupby('time.season').sum()
if season == "all":
## calculate weighted mean for all season
smean = (ds * weights).groupby('time.season').mean('time')
else :
## calculate weighted mean for specified season
smean = (ds * weights).where(ds['time.season'] == season).mean('time')
return smean[season]
def season_ts(ds, season):
""" calculate timeseries of seasonal averages
Args: ds (xarray.Dataset): dataset
var (str): variable to calculate
season (str): 'DJF', 'MAM', 'JJA', 'SON'
"""
## set months outside of season to nan
ds_season = ds.where(ds['time.season'] == season)
# calculate 3month rolling mean (only middle months of season will have non-nan values)
ds_season = ds_season.rolling(min_periods=3, center=True, time=3).mean()
# reduce to one value per year
ds_season = ds_season.groupby('time.year').mean('time')
# FUTURE: remove first year if it has nan?
return ds_season
In [147]:
del gb
In [126]:
# creating nino34 index sst timeseries for DJF month only
tos_djf_anom_global=season_ts(tos_anom,'DJF')
tos_nino34_djf_anom=tos_djf_anom_global.sel(lat=slice(5,-5),lon=slice(190,240)).mean(dim=['lat','lon'])
In [126]:
In [126]:
In [127]:
# Now let's plot the Nino 3.4 Index
fig = plt.figure(figsize=(12, 6))
tos_nino34_djf_anom.plot(color='black')
plt.axhline(0, color='black', lw=0.5)
plt.axhline(0.5, color='black', linewidth=0.5, linestyle='dotted')
plt.axhline(-0.5, color='black', linewidth=0.5, linestyle='dotted')
plt.title('Niño 3.4 Index');
In [128]:
## find elnino years
elnino_year=tos_nino34_djf_anom.year[tos_nino34_djf_anom.data>0.5]
## find lanina years
lanina_year=tos_nino34_djf_anom.year[tos_nino34_djf_anom.data<-0.5]
print('elnino, ',elnino_year.data)
print('lanina, ',lanina_year.data)
In [129]:
# elnino anomaly compoite mean DJF
elnino_anom_compmean=tos_djf_anom_global.sel(year=elnino_year.data).mean(dim='year')
# lanina anomaly compoite mean DJF
lanina_anom_compmean=tos_djf_anom_global.sel(year=lanina_year.data).mean(dim='year')
In [130]:
fig, (ax) = plt.subplots(
ncols=1, nrows=3, figsize=[10, 12], subplot_kw={"projection": ccrs.Robinson()}
)
# Elnino
elnino_anom_compmean.plot(
ax=ax[0],
x="lon",
y="lat",
transform=ccrs.PlateCarree(),
vmin=-4,
vmax=4,
cmap="coolwarm",
robust=True,
)
ax[0].coastlines()
ax[0].set_global()
ax[0].set_title("Composite mean Elnino SST anomalies (°C)")
# Lanina
lanina_anom_compmean.plot(
ax=ax[1],
x="lon",
y="lat",
transform=ccrs.PlateCarree(),
vmin=-4,
vmax=4,
cmap="coolwarm",
robust=True,
)
ax[1].coastlines()
ax[1].set_global()
ax[1].set_title("Composite mean Lanina SST anomalies (°C)")
# Elnino - Lanina
(elnino_anom_compmean-lanina_anom_compmean).plot(
ax=ax[2],
x="lon",
y="lat",
transform=ccrs.PlateCarree(),
vmin=-4,
vmax=4,
cmap="coolwarm",
robust=True,
)
ax[2].coastlines()
ax[2].set_global()
ax[2].set_title("Composite mean Elnino - Lanina SST anomalies (°C)")
Out[130]:
In [130]:
In [130]:
In [130]: