In [1]:
%matplotlib inline
from matplotlib import pyplot as plt
In [2]:
import os, sys
import numpy as np
from numpy import ma
import pandas as pd
import xray
from scipy.signal import detrend
In [3]:
dset = xray.open_dataset('/Users/nicolasf/data/SST/ERSST/ersst.realtime.nc')
In [4]:
dset
Out[4]:
In [5]:
dset['lon']
Out[5]:
In [6]:
dset['lat']
Out[6]:
In [7]:
dset = dset.sel(zlev=0, time=slice('1948','2014'))
In [8]:
lat = dset['lat'].data
lon = dset['lon'].data
In [9]:
lon
Out[9]:
In [10]:
def get_mask(var):
mat = var.data
X = np.reshape(mat, (mat.shape[0], mat.shape[1] * mat.shape[2]), order='F')
X = ma.masked_array(X, np.isnan(X))
mask = X.sum(0).mask
mask = mask.reshape((mat.shape[1], mat.shape[2]), order='F')
return mask
In [17]:
mask = get_mask(dset['sst'])
In [18]:
mask
Out[18]:
In [15]:
#mask = np.logical_not(mask)
In [19]:
plt.imshow(mask); plt.colorbar()
Out[19]:
In [27]:
mask = mask.astype(np.int)
In [29]:
plt.imshow(mask); plt.colorbar()
Out[29]:
In [28]:
mask
Out[28]:
In [38]:
d = {}
d['time'] = ('time',dset['time'])
d['latitudes'] = ('latitudes',lat)
d['longitudes'] = ('longitudes', lon)
d['mask'] = (['latitudes','longitudes'], mask)
d['sst'] = (['time','latitudes','longitudes'], dset['sst'].data)
In [39]:
dset_out = xray.Dataset(d)
In [40]:
dset_out['mask'].plot(yincrease=True)
Out[40]:
In [33]:
plt.imshow(dset_out['sst'][0,:,:]); plt.colorbar()
Out[33]:
In [41]:
dset_out.to_netcdf('../data/ERSST_monthly_SST_1948_2014.nc')
In [ ]:
def demean(x):
return x - x.sel(time=slice('1981-1-1','2010-12-1')).mean('time')
In [ ]:
sst_anomalies = dset['sst'].groupby('time.month').apply(demean)
In [ ]:
def detrend_var(var):
"""
var must be a xray dataarray
"""
mat = var.data
X = np.reshape(mat, (mat.shape[0], mat.shape[1] * mat.shape[2]), order='F')
X = ma.masked_array(X, np.isnan(X))
land = X.sum(0).mask
ocean = -land
Xocean = X[:,ocean]
Xocean = detrend(Xocean, axis=0) + Xocean.mean(0)
Xdetrend = np.ones(X.shape) * -999.
Xdetrend[:,ocean] = Xocean
Xdetrend = np.reshape(Xdetrend, mat.shape, order='F')
Xdetrend = ma.masked_values(Xdetrend, -999.)
mask = land.reshape((mat.shape[1], mat.shape[2]), order='F')
return mask, Xdetrend
In [ ]:
mask, sst_detrend = detrend_var(dset['sst'])
In [ ]:
mask, sst_anomalies_detrend = detrend_var(sst_anomalies)
In [ ]:
plt.imshow(sst_anomalies_detrend[0,:,:])
In [ ]:
plt.imshow(sst_detrend[-2,::-1,:], vmin=0, vmax=30)
In [ ]:
plt.plot(sst_detrend[:,50,100])
plt.plot(dset['sst'][:,50,100])
In [ ]:
plt.plot()
In [ ]:
detrend?
In [ ]:
X = np.reshape(mat, (mat.shape[0], len(lat) * len(lon)), order='F')
In [ ]:
X = ma.masked_array(X, np.isnan(X))
In [ ]:
land = X.sum(0).mask
In [ ]:
ocean = -land
In [ ]:
Xocean = X[:,ocean]
In [ ]:
Xocean.shape
In [ ]:
Xocean = detrend(Xocean, axis=0)
In [ ]:
Xdetrend = np.ones(X.shape) * -999.
In [ ]:
Xdetrend[:,ocean] = Xocean
In [ ]:
Xdetrend.shape
In [ ]:
X.shape
In [ ]:
Xdetrend = np.reshape(Xdetrend, mat.shape, order='F')
In [ ]:
Xdetrend = ma.masked_values(Xdetrend, -999.)
In [ ]:
plt.imshow(Xdetrend[0,::-1,:])
In [ ]:
mask = land.reshape((len(lat), len(lon)), order='F')
In [ ]:
In [ ]:
d = {}
d['time'] = ('time',dset['time'])
d['lat'] = ('latitudes',lat)
d['lon'] = ('longitudes', lon)
d['mask'] = (['latitudes','longitudes'], mask)
d['sst'] = (['time','latitudes','longitudes'], dset['sst'].data)
d['sst_detrend'] = (['time','latitudes','longitudes'], dset['sst'].data)
d['sst_anomalies'] = (['time','latitudes','longitudes'], sst_anoms)
d['sst_anomalies_detrend'] = (['time','latitudes','longitudes'], sstd)
In [ ]:
dset_detrend = xray.Dataset(d)
In [ ]:
dset_detrend
In [ ]:
dset_detrend.to_netcdf('../data/ERSST_anoms_Detrend_1979_2010.nc')
In [ ]:
sst_seas
In [ ]:
import bottleneck as bn
In [ ]:
bn.move_mean?
In [ ]:
m3 = bn.move_mean(dset_detrend['sstd'], 3, axis=0)
m6 = bn.move_mean(dset_detrend['sstd'], 6, axis=0)
In [ ]:
m3.shape
In [ ]:
dset_detrend
In [ ]:
plt.imshow(m6[5,::-1,:])
In [ ]: