An example is presented here of a numerical analysis of Sea Surface Temperature (SST) where the global change rate per decade has been calculated from 1982 to 2015 using the daily NCEI AVHRR_OI Level 4 v2 dataset (https://podaac.jpl.nasa.gov/dataset/AVH ... -GLOB-v2.0). A linear trend has been fitted to each data location for each day in the SST time series and the anomaly has been calculated as described in the post https://podaac.jpl.nasa.gov/forum/viewtopic.php?f=23&t=346.
The python code is shown here:
Trend Analysis
- Code: Select all
#! /usr/bin/env python2
#
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
import os
import glob
import calendar
#*** Setup Global Variables
nlon = 1440
nlat = 720
startY = 1982
endY = 2015
#*** Setup Data File Directory Name
dirname = './'
#*** Calculate how many days of data files
nt = 0
for i in range(startY,endY+1):
if calendar.isleap(i):
nt = nt + 366
else:
nt = nt + 365
sst_all=np.empty((nt,nlat/2,nlon/2))
# Pre-define latitude and longitude grid
lats = [-89.875+0.5*x for x in range(0, nlat/2)]
lons = [-179.875+0.5*x for x in range(0, nlon/2)]
# Read SST using the netCDF4 module.
idx = 0
for i in range(startY, endY+1):
for j in range(1, 366):
for filename in glob.glob(dirname+str(i)+'/'+"{0:0>3}".format(j)+'/*.nc'):
ncin = Dataset(filename, 'r')
sst = ncin.variables['analysed_sst'][:]
sst_all[idx,:,:] = sst[0,0:nlat:2,0:nlon:2]
ncin.close()
# Linear trend calculation
sst=np.empty((idx))
x = range(idx)
sst_rate=np.empty((nlat/2,nlon/2))
for i in range(nlon/2):
for j in range(nlat/2):
sst[:] = sst_all[0:idx,j,i]
z = np.polyfit(x, sst, 1)
sst_rate[j,i] = z[0]*3650.0
#*** Plot the result
plt.figure()
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 = plt.colorbar(cs, orientation='horizontal')
cb.set_label('SST Changing Rate (deg/decade)', fontsize=12)
plt.title('SST Changing Rate (degC/decade)', fontsize=16)
plt.show()
Anomaly Analysis
- Code: Select all
#! /usr/bin/env python2
#
from netCDF4 import Dataset
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
import os
import glob
import calendar
#*** Setup Global Variables
nlon = 1440
nlat = 720
startY = 1982
endY = 2015
#*** Setup Data File Directory Name
dirname = './'
#*** Calculate how many days of data files
nt = 0
for i in range(startY,endY+1):
if calendar.isleap(i):
nt = nt + 366
else:
nt = nt + 365
#*** Array defination
sst_glb_avg=np.empty(nt)
sst_glb_anomaly=np.empty(nt)
sst_glb_std=np.empty(nt)
sst_year_day=np.empty((endY-startY+1,366))
sst_season=np.empty((366))
sst_season_long=np.empty((nt))
# Read SST using the netCDF4 module.
idx = 0
for i in range(startY, endY+1):
for j in range(1, 367):
for filename in glob.glob(dirname+str(i)+'/'+"{0:0>3}".format(j)+'/*.nc'):
ncin = Dataset(filename, 'r')
sst = ncin.variables['analysed_sst'][:]
ncin.close()
asst = sst[sst > 0]
sst_glb_avg[idx] = np.average(asst)
sst_glb_std[idx] = np.std(asst)
sst_year_day[i-startY,j-1] = sst_glb_avg[idx]
idx = idx + 1
sst_glb_avg_new = sst_glb_avg[0:idx]
sst_glb_std_new = sst_glb_std[0:idx]
#*** Seasonal Cycle
for j in range(0, 366):
sst_season[j] = 0.0
np = 0
for i in range(0, endY-startY+1):
if sst_year_day[i,j] > 0:
sst_season[j] = sst_season[j] + sst_year_day[i,j]
np = np + 1
sst_season[j] = sst_season[j]/float(np)
np = 0
for i in range(startY,endY+1):
if calendar.isleap(i):
np = np + 366
sst_season_long[np-366:np] = sst_season[0:366]
else:
np = np + 365
sst_season_long[np-365:np] = sst_season[0:365]
#*** Calculate anomaly by removing seasonal cycle
sst_glb_anomaly[:] = sst_glb_avg[:] - sst_season_long[:]
#*** PLot the results
plt.figure()
# Plot the time series
plt.figure()
days = range(0, idx)
plt.plot(sst_glb_avg, color='b', linewidth=2)
plt.axhline(0, color='k')
plt.title('SST Anomaly Time Series')
plt.xlabel('Day #')
plt.ylabel('Normalized Units')
plt.ylim(min(sst_glb_avg), max(sst_glb_avg))
plt.show()
User could also download the Docker image for this recipe from dockerhub, and run the recipe using jupyter nootbook:
https://hub.docker.com/r/podaacdatarecipes/trend_analysis/