In [2]:
import xarray as xr
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import timeit 
import dask
# Note install catopy and shapely from Gohlke's wheels using pip
import cartopy.crs as ccrs
import bokeh.plotting as bk
# plotly dies not seem to have a Python 3.5 combatible version, unfortunately.

Loading dataset


In [3]:
fn_chl = 'C:/Users/gunnar/src/CABLAB/xarray_test/data/chl_a/*.nc'
chl =  xr.open_mfdataset(fn_chl)

Subsetting a region: Baltic Sea


In [4]:
# For the sake of the example, we selct two arrays from the same data set, chorophyll a and the log10 bias.
chla  = chl.chlor_a.loc[dict(lat=slice(90.,50.),lon=slice(7.,30.))]
chl_rms = chl.chlor_a_log10_bias.loc[dict(lat=slice(90.,50.),lon=slice(7.,30.))]

In [5]:
chla.load()
chl_rms.load()


Out[5]:
<xarray.DataArray 'chlor_a_log10_bias' (time: 12, lat: 960, lon: 552)>
array([[[ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        ..., 
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan]],

       [[ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        ..., 
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan]],

       [[ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        ..., 
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan]],

       ..., 
       [[ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        ..., 
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan]],

       [[ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        ..., 
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan]],

       [[ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        ..., 
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan],
        [ nan,  nan,  nan, ...,  nan,  nan,  nan]]])
Coordinates:
  * lon      (lon) float32 7.02083 7.0625 7.10417 7.14583 7.1875 7.22917 ...
  * lat      (lat) float32 89.9792 89.9375 89.8958 89.8542 89.8125 89.7708 ...
  * time     (time) datetime64[ns] 2010-01-01 2010-02-01 2010-03-01 ...
Attributes:
    comment: Uncertainty lookups derived from file: /data/datasets/CCI/v2.0-production/stage09b-uncertainty_tables/chlor_a/cci_chla_bias.dat
    grid_mapping: crs
    ref: http://www.uncertweb.org
    rel: uncertainty
    long_name: Bias of log10-transformed chlorophyll-a concentration in seawater.

Correlation Analysis


In [6]:
# defining method to calculate cross-correlation between two data-sets along a given axis.
def mycorr(x,y,dim=None):
    return (((x-x.mean(dim=dim))*(y-y.mean(dim=dim))).sum(dim=dim)/x.count(dim=dim))/(x.std(dim=dim)*y.std(dim=dim))

In [7]:
# calling the function
corr_chla = mycorr(chla,chl_rms,dim="time")

In [8]:
# plotting the result as a 2D overlay on a map with coastlines. The normalized correlation in time varies between -1 and 1

map_extent = [0, 35, 40, 80]
dtick = 10.
crs = ccrs.NorthPolarStereo()
ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()
ax.coastlines('50m')
ax.set_extent(map_extent, ccrs.PlateCarree())  
ax.set_xticks(np.arange(map_extent[0],map_extent[1],dtick), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(map_extent[2],map_extent[3],dtick), crs=ccrs.PlateCarree())
corr_chla.plot.contourf(ax=ax, transform=ccrs.PlateCarree(),levels = np.linspace(-1.,1.,20.),cmap="hot")


Out[8]:
<matplotlib.contour.QuadContourSet at 0x225014f0128>

In [10]:
# plotting the result as a 2D overlay on a map with coastlines. The normalized correlation in time varies between -1 and 1
map_extent = [1., 360, 50, 90]
dtick = 10.
ax = plt.axes(projection=ccrs.NorthPolarStereo())

ax.set_extent(map_extent, crs=ccrs.PlateCarree())
#ax.set_xticks(np.arange(map_extent[0],map_extent[1],dtick), crs=ccrs.PlateCarree())
#ax.set_yticks(np.arange(map_extent[2],map_extent[3],dtick), crs=ccrs.PlateCarree())
#ax.plot(x,y,transform=ccrs.Geodetic())
corr_chla.plot.contourf(ax=ax, transform=ccrs.PlateCarree(),levels = np.linspace(-1.,1.,20.),cmap="hot")
ax.stock_img()
#ax.coastlines()


Out[10]:
<matplotlib.image.AxesImage at 0x2250f3e98d0>

In [11]:
# monthly plot of chla only, standard xarray plotting
chla.plot(x='lon', y='lat', col='time', col_wrap=3, vmax = 7)


Out[11]:
<xarray.plot.facetgrid.FacetGrid at 0x22503596940>

In [12]:
def plot_TimeSeries(data,lat,lon,tdim = "time", latdim="lat",londim="lon", xlabel = None, ylabel =None, title=None):
    xx = data.sel(**{latdim : lat, londim :lon}, method='nearest')
    y = xx.values
    t = xx[tdim].values
    y_mean = data.mean(dim=[latdim,londim]).values
    # output to static HTML file:
    bk.output_notebook()

    # create a new plot
    TOOLS="resize,crosshair,pan,wheel_zoom,box_zoom,reset,box_select,lasso_select"

    p = bk.figure(
       tools=TOOLS, title=title,
       x_axis_label=xlabel, y_axis_label=ylabel, x_axis_type="datetime"
    )

    # add some renderers
    p.line(t, y, legend='Time Series at Latitude '+str(lat)+', Longitude '+str(lon))
    p.line(t, y_mean, legend='Mean in the entire region', color="red")
    # show the results
    bk.show(p)

In [13]:
llat = 55.
llon = 11.
plot_TimeSeries(chla,llat,llon, title="MyFigure", xlabel = "Time", ylabel = "Chl a")


Loading BokehJS ...

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]: