% 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')

"""
Explanation: Trend and Anomaly Analyses of Long-term Tempro-Spatial Dataset
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.

Data Source
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.



1. Load basic libs
End of explanation
"""


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

"""
Explanation: 2. Read SST data and pick variables
2.1 Read SST
End of explanation
"""


utime = netcdftime.utime(t_unit, calendar = t_cal)
datevar = utime.num2date(nctime)
print(datevar.shape)

datevar[0:5]

"""
Explanation: 2.2 Parse time
End of explanation
"""


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

"""
Explanation: 2.3 Read mask (1=ocen, 0=land)
End of explanation
"""


sst = np.ma.masked_array(sst, mask=lsm)

"""
Explanation: 2.3 Mask out Land
End of explanation
"""


#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')

"""
Explanation: 3. Trend Analysis
3.1 Linear trend calculation
End of explanation
"""


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)

"""
Explanation: 3.2 Visualize SST trend
End of explanation
"""


sst_grd_ym  = sst.reshape((12,nt/12, ngrd), order='F').transpose((1,0,2))
sst_grd_ym.shape

"""
Explanation: 4. Anomaly analysis
4.1 Convert sst data into nyear x12 x lat x lon
End of explanation
"""


sst_grd_clm = np.mean(sst_grd_ym, axis=0)
sst_grd_clm.shape

"""
Explanation: 4.2 Calculate seasonal cycle
End of explanation
"""


sst_grd_anom = (sst_grd_ym - sst_grd_clm).transpose((1,0,2)).reshape((nt, nlat, nlon), order='F')
sst_grd_anom.shape

"""
Explanation: 4.3 Remove seasonal cycle
End of explanation
"""


print(lats[0:12])
print(lons[0:12])

"""
Explanation: 4.4 Calculate area-weights
4.4.1 Make sure lat-lon grid direction
End of explanation
"""


lonx, latx = np.meshgrid(lons, lats)
weights = np.cos(latx * np.pi / 180.)

print(weights.shape)

"""
Explanation: 4.4.2 Calculate area-weights with cos(lats)
End of explanation
"""


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,:])    

"""
Explanation: 4.4.3 Calculate valid grids total eareas for Global, NH and SH
End of explanation
"""


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')

"""
Explanation: 5. Visualize monthly SST anomaly time series
End of explanation
"""


