## Trend and Anomaly Analyses of Long Time Series Datasets

### Trend and Anomaly Analyses of Long Time Series Datasets

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 Datasetimport matplotlib.pyplot as pltfrom mpl_toolkits.basemap import Basemapimport numpy as npimport osimport globimport calendar#*** Setup Global Variablesnlon = 1440nlat = 720startY = 1982endY = 2015#*** Setup Data File Directory Namedirname = './'#*** Calculate how many days of data filesnt = 0for i in range(startY,endY+1):   if calendar.isleap(i):     nt = nt + 366   else:     nt = nt + 365sst_all=np.empty((nt,nlat/2,nlon/2))# Pre-define latitude and longitude gridlats = [-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 = 0for 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 calculationsst=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 resultplt.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 Datasetimport matplotlib.pyplot as pltfrom mpl_toolkits.basemap import Basemapimport numpy as npimport osimport globimport calendar#*** Setup Global Variablesnlon = 1440nlat = 720startY = 1982endY = 2015#*** Setup Data File Directory Namedirname = './'#*** Calculate how many days of data filesnt = 0for i in range(startY,endY+1):   if calendar.isleap(i):     nt = nt + 366   else:     nt = nt + 365#*** Array definationsst_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 = 0for 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 + 1sst_glb_avg_new = sst_glb_avg[0:idx]sst_glb_std_new = sst_glb_std[0:idx]#*** Seasonal Cyclefor 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 = 0for 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 cyclesst_glb_anomaly[:] = sst_glb_avg[:] - sst_season_long[:]#*** PLot the resultsplt.figure()# Plot the time seriesplt.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/
yiboj

Posts: 130
Joined: Mon Mar 30, 2015 11:22 am

### Re: Trend and Anomaly Analyses of Long Time Series Datasets

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.
esherman

Posts: 1
Joined: Thu Aug 10, 2017 9:46 am

### Re: Trend and Anomaly Analyses of Long Time Series Datasets

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.
yiboj

Posts: 130
Joined: Mon Mar 30, 2015 11:22 am