In [1]:
%load_ext load_style
%load_style talk.css
Trend and anomaly analyses are widely used in atmospheric and oceanographic research for detecting long term change.
An example is presented in this notebook of a numerical analysis of Sea Surface Temperature (SST) where the global change rate per decade has been calculated from 1982 to 2016. Moreover, its area-weighted global monthly SST anomaly time series is presented, too. In addition, all of calculating processes is list step by step.
NOAA Optimum Interpolation (OI) Sea Surface Temperature (SST) V2 is downloaded from https://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.html.
Spatial Coverage:
- 1.0 degree latitude x 1.0 degree longitude global grid (180x360).
- 89.5N - 89.5S, 0.5E - 359.5E.
Because oisst is an interpolated data, so it covers ocean and land. As a result, have to use land-ocean mask data at the same time, which is also available from the webstie.
We select SST from 1982 to 2016.
In [2]:
% matplotlib inline
from pylab import *
import numpy as np
import datetime
from netCDF4 import netcdftime
from netCDF4 import Dataset as netcdf # netcdf4-python module
from netcdftime import utime
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import matplotlib.dates as mdates
from matplotlib.dates import MonthLocator, WeekdayLocator, DateFormatter
import matplotlib.ticker as ticker
from matplotlib.pylab import rcParams
rcParams['figure.figsize'] = 15, 6
import warnings
warnings.simplefilter('ignore')
In [3]:
ncset= netcdf(r'data/sst.mnmean.nc')
lons = ncset['lon'][:]
lats = ncset['lat'][:]
sst = ncset['sst'][1:421,:,:] # 1982-2016 to make it divisible by 12
nctime = ncset['time'][1:421]
t_unit = ncset['time'].units
try :
t_cal =ncset['time'].calendar
except AttributeError : # Attribute doesn't exist
t_cal = u"gregorian" # or standard
nt, nlat, nlon = sst.shape
ngrd = nlon*nlat
In [4]:
utime = netcdftime.utime(t_unit, calendar = t_cal)
datevar = utime.num2date(nctime)
print(datevar.shape)
datevar[0:5]
Out[4]:
In [5]:
lmfile = 'data\lsmask.nc'
lmset = netcdf(lmfile)
lsmask = lmset['mask'][0,:,:]
lsmask = lsmask-1
num_repeats = nt
lsm = np.stack([lsmask]*num_repeats,axis=-1).transpose((2,0,1))
lsm.shape
Out[5]:
In [6]:
sst = np.ma.masked_array(sst, mask=lsm)
In [7]:
#import scipy.stats as stats
sst_grd = sst.reshape((nt, ngrd), order='F')
x = np.linspace(1,nt,nt)#.reshape((nt,1))
sst_rate = np.empty((ngrd,1))
sst_rate[:,:] = np.nan
for i in range(ngrd):
y = sst_grd[:,i]
if(not np.ma.is_masked(y)):
z = np.polyfit(x, y, 1)
sst_rate[i,0] = z[0]*120.0
#slope, intercept, r_value, p_value, std_err = stats.linregress(x, sst_grd[:,i])
#sst_rate[i,0] = slope*120.0
sst_rate = sst_rate.reshape((nlat,nlon), order='F')
In [8]:
m = Basemap(projection='cyl', llcrnrlon=min(lons), llcrnrlat=min(lats),
urcrnrlon=max(lons), urcrnrlat=max(lats))
x, y = m(*np.meshgrid(lons, lats))
clevs = np.linspace(-0.5, 0.5, 21)
cs = m.contourf(x, y, sst_rate.squeeze(), clevs, cmap=plt.cm.RdBu_r)
m.drawcoastlines()
#m.fillcontinents(color='#000000',lake_color='#99ffff')
cb = m.colorbar(cs)
cb.set_label('SST Changing Rate ($^oC$/decade)', fontsize=12)
plt.title('SST Changing Rate ($^oC$/decade)', fontsize=16)
Out[8]:
In [9]:
sst_grd_ym = sst.reshape((12,nt/12, ngrd), order='F').transpose((1,0,2))
sst_grd_ym.shape
Out[9]:
In [10]:
sst_grd_clm = np.mean(sst_grd_ym, axis=0)
sst_grd_clm.shape
Out[10]:
In [11]:
sst_grd_anom = (sst_grd_ym - sst_grd_clm).transpose((1,0,2)).reshape((nt, nlat, nlon), order='F')
sst_grd_anom.shape
Out[11]:
In [12]:
print(lats[0:12])
print(lons[0:12])
In [13]:
lonx, latx = np.meshgrid(lons, lats)
weights = np.cos(latx * np.pi / 180.)
print(weights.shape)
In [14]:
sst_glb_avg = np.zeros(nt)
sst_nh_avg = np.zeros(nt)
sst_sh_avg = np.zeros(nt)
for it in np.arange(nt):
sst_glb_avg[it] = np.ma.average(sst_grd_anom[it, :], weights=weights)
sst_nh_avg[it] = np.ma.average(sst_grd_anom[it,0:nlat/2,:], weights=weights[0:nlat/2,:])
sst_sh_avg[it] = np.ma.average(sst_grd_anom[it,nlat/2:nlat,:], weights=weights[nlat/2:nlat,:])
In [15]:
fig, ax = plt.subplots(1, 1 , figsize=(15,5))
ax.plot(datevar, sst_glb_avg, color='b', linewidth=2, label='GLB')
ax.plot(datevar, sst_nh_avg, color='r', linewidth=2, label='NH')
ax.plot(datevar, sst_sh_avg, color='g', linewidth=2, label='SH')
ax.axhline(0, linewidth=1, color='k')
ax.legend()
ax.set_title('Monthly SST Anomaly Time Series (1982 - 2016)', fontsize=16)
ax.set_xlabel('Month/Year #', fontsize=12)
ax.set_ylabel('$^oC$', fontsize=12)
ax.set_ylim(-0.6, 0.6)
fig.set_figheight(9)
# rotate and align the tick labels so they look better
fig.autofmt_xdate()
# use a more precise date string for the x axis locations in the toolbar
ax.fmt_xdata = mdates.DateFormatter('%Y')
http://unidata.github.io/netcdf4-python/
Kalnay et al.,The NCEP/NCAR 40-year reanalysis project, Bull. Amer. Meteor. Soc., 77, 437-470, 1996.
Matplotlib: A 2D Graphics Environment by J. D. Hunter In Computing in Science & Engineering, Vol. 9, No. 3. (2007), pp. 90-95
Reynolds, R.W., N.A. Rayner, T.M. Smith, D.C. Stokes, and W. Wang, 2002: An improved in situ and satellite SST analysis for climate. J. Climate, 15, 1609-1625.
In [ ]: