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.


    

In [3]:

    

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


    

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.

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>