PODAAC

Parameter trend and anomaly analyses are widely used in atmospheric and oceanographic research for detecting long term change.

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/

Line 48 of the Anomaly Analysis code takes a global average (sst_glb_avg) of sea surface temperature, but does not weight grid cell values as a function of cosine(degrees latitude). Why are these values not weighted?

Lines 59-64 computes the average SST everyday for the given time period (sst_season), creating a climatology of everyday of the year. The global average values (sst_glb_avg) computed in line 48 are used to calculate the climatology. As in, the global average values for January 1st is averaged over the time period to create the January 1st climatology. These lines of code are computing the average of an average, which may not be the same value if you computed the climatology using the raw daily sst data. If the number of observations for each value of sst_glb_avg differs, I suspect taking the average of the global averages (sst_glb_avg) is not computing the climatology you want.

Hi,
Thanks for the comments, it is very helpful for us to improve the algorithm later.

As you suggested, the value should be weighted by cos(latitude), and the climatology should be recalculated. We will update the code accordingly.

Thanks again for your valuable inputs.