In [31]:
%matplotlib inline
import os
import numpy as np
import matplotlib.pyplot as mplt
import xarray as xr
import xarray.plot as xplt
import pandas as pd
import bottleneck as bn
from ipywidgets import interact
ecv_base = 'D:\\EOData\\CCI-TBX'
In [32]:
oc_pattern = os.path.join(ecv_base, 'occci-v2.0/data/geographic/netcdf/monthly/chlor_a/2010/*.nc')
st_pattern = os.path.join(ecv_base, 'sst/data/lt/Analysis/L4/v01.1/2010/*.nc')
oc_ds = xr.open_mfdataset(oc_pattern)
sst_ds = xr.open_mfdataset(st_pattern)
In [33]:
chla = oc_ds.chlor_a
In [34]:
chla
Out[34]:
In [35]:
xplt.imshow(chla.isel(time=0), vmax=0.3)
Out[35]:
In [36]:
sst = sst_ds.analysed_sst
xarray has a number of issues with SST data, see XXX
In [37]:
sst
Out[37]:
In [38]:
xr.plot.imshow(sst_ds.analysed_sst.isel(time=2))
Out[38]:
In [39]:
import scipy.ndimage
lat_factor = 3600./4320.
lon_factor = 7200./8640.
def resample(x):
y = scipy.ndimage.zoom(x, [lat_factor, lon_factor])
return xr.DataArray(y)
# Help! This is soooo slow... few minutes on Norman's PC
chla_temp = chla.groupby('time').apply(resample)
lon = scipy.ndimage.zoom(chla.lon, [lon_factor])
lat = scipy.ndimage.zoom(chla.lat, [lat_factor])
chla = xr.DataArray(chla_temp, dims=['time', 'lat', 'lon'],
coords=dict(time=chla.time, lat=lat, lon=lon), attrs=chla.attrs)
chla
Out[39]:
Note of OC data we have to supply lat slices in reverse order!
In [43]:
chla = chla.sel(lat=slice(45., 30.), lon=slice(-60., -45.))
chla
Out[43]:
In [44]:
sst = sst.sel(lat=slice(30., 45.), lon=slice(-60., -45.))
sst
Out[44]:
In [46]:
def cov(a, b):
#print(a.shape)
x = (a - bn.nanmean(a)) * (b - bn.nanmean(b))
n = np.count_nonzero(np.isnan(x))
return bn.nansum(x) / n if n > 0 else np.nan
def autocov(x, lag):
#print(x, type(x))
a = x[:-lag]
b = x[lag:]
return cov(a, b)
def myfunc(x, axis, **kwargs):
#print(x.shape, axis)
res = np.apply_along_axis(autocov, axis, x, 1)
#print(res.shape)
return res
chla_auto_corr = chla.reduce(myfunc, dim='time', keep_attrs=True)
In [13]:
plot.imshow(chla_auto_corr)
In [ ]: