

In [1]:

    
import xarray as xr
import oocgcm.filtering.linearfilters as tf
import pylab as plt
import numpy as np
%matplotlib inline



In [2]:

    
sigdir = "../examples/signals/"
test_file = 'ORCA025_SSH_SWPAC.nc'

Here is an example of simulated sea surface height from the NEMO model run at 1/4°, that is represented by an xarray.DataArray object. Note that dask array can be used by precising chunks.



In [24]:

    
signal_xyt = xr.open_dataset(sigdir + test_file, decode_times=False)['sossheig'].chunk(chunks={'time_counter': 50})
print signal_xyt
signal_xyt.isel(time_counter=0).plot(vmin=-0.07, vmax=0.07, cmap='seismic')









    



<xarray.DataArray 'sossheig' (time_counter: 300, y: 184, x: 224)>
dask.array<xarray-..., shape=(300, 184, 224), dtype=float64, chunksize=(50, 184, 224)>
Coordinates:
  * time_counter  (time_counter) float64 73.04 nan nan nan nan nan nan nan ...
  * y             (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * x             (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...






    Out[24]:





<matplotlib.collections.QuadMesh at 0x7fe7d6d9b1d0>

A window object that is linked to the previous data is created by the following line:



In [4]:

    
win1D = signal_xyt.win

Temporal filtering: Boxcar window

A boxcar window object that will be applied along the time dimension is simply be defined by setting its different properties:



In [5]:

    
win1D.set(window_name='boxcar', n=[5], dims=['time_counter'])
print win1D._depth.values()
print win1D









    



[5]
Window [name->boxcar,dims->['time_counter'],order->{'time_counter': 5}]

The weight distribution and the frequency response of the window may be plotted for one dimensional windows using the plot method:



In [6]:

    
win1D.plot()

This window can now be applied on the data by the associated procedure:



In [7]:

    
signal_LF_box = win1D.apply(compute=True)









    



[########################################] | 100% Completed |  0.6s






    



/home/serazin/anaconda2/envs/test_env/lib/python2.7/site-packages/xarray/core/variable.py:1046: RuntimeWarning: invalid value encountered in divide
  if not reflexive

By default the filtering is computed when using the apply function. If compute is set to False, the filtering will be computed when the output data is required. This allows the definition of several dask object before the global computation is performed.

Temporal filtering: Lanczos window

Setting now different properties to use a Lanczos window:



In [8]:

    
win1D.set(window_name='lanczos', n=[5], dims=['time_counter'], fc=0.1)
print win1D
win1D.plot()









    



Window [name->lanczos,dims->['time_counter'],order->{'time_counter': 5}]






    



/home/serazin/anaconda2/envs/test_env/lib/python2.7/site-packages/oocgcm-0.0.1-py2.7.egg/oocgcm/filtering/linearfilters.py:56: RuntimeWarning: invalid value encountered in divide



In [9]:

    
signal_LF_lcz = win1D.apply(compute=True)









    



[########################################] | 100% Completed |  0.6s

The following figure gives the comparison between the boxcar and the lanczos window applied for the low-pass filtering at one grid point of the dataset.



In [10]:

    
signal_xyt.isel(x=50, y=50).plot()
signal_LF_box.isel(x=50, y=50).plot(color='g')
signal_LF_lcz.isel(x=50, y=50).plot(color='r')
plt.legend(['raw', 'boxcar 10yr', 'lanczos 10yr'])









    Out[10]:





<matplotlib.legend.Legend at 0x7fe7d9a7de10>

Spatial filtering

The window object extend to multidimensional filtering such bidimensional spatial filtering. The filtering method used is able to deal with missing data or coastline by reweighting the filters weights. For example, a 2D Lanczos window may thus be associated to a dataset.



In [11]:

    
signal_xyt = xr.open_dataset(sigdir + test_file, decode_times=False)['sossheig'].chunk(chunks={'x': 40, 'y':40})
win_box2D = signal_xyt.win
win_box2D.set(window_name='lanczos', n=[24, 24], dims=['x', 'y'], fc=0.0416)
print signal_xyt
print win_box2D









    



<xarray.DataArray 'sossheig' (time_counter: 300, y: 184, x: 224)>
dask.array<xarray-..., shape=(300, 184, 224), dtype=float64, chunksize=(300, 40, 40)>
Coordinates:
  * time_counter  (time_counter) float64 73.04 nan nan nan nan nan nan nan ...
  * y             (y) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
  * x             (x) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
Window [name->lanczos,dims->['x', 'y'],order->{'y': 24, 'x': 24}]



In [12]:

    
win_box2D.plot()









    



{1: 24, 2: 24}
(49, 49)

Original dataset

Temporal standard deviation of the raw timeseries:



In [13]:

    
bw = win_box2D.boundary_weights()



In [14]:

    
bw.isel(time_counter=1).plot(vmin=0, vmax=1, cmap='spectral')









    Out[14]:





<matplotlib.collections.QuadMesh at 0x7fe7d9b6ef50>



In [15]:

    
signal_LS = win_box2D.apply(weights=bw)









    



[########################################] | 100% Completed |  3min 22.0s



In [16]:

    
# Original dataset
signal_xyt.std(dim='time_counter').plot(cmap='jet', vmin=0, vmax=0.06)









    



/home/serazin/anaconda2/envs/test_env/lib/python2.7/site-packages/dask/array/reductions.py:265: RuntimeWarning: invalid value encountered in true_divide
  u = total/n
/home/serazin/anaconda2/envs/test_env/lib/python2.7/site-packages/dask/array/numpy_compat.py:44: RuntimeWarning: invalid value encountered in divide
  x = np.divide(x1, x2, out)






    Out[16]:





<matplotlib.collections.QuadMesh at 0x7fe7d98b6e90>

Large Scales (>6°)

Temporal standard deviation of the large field:



In [17]:

    
# Large-scale (>6°) dataset
signal_LS.std(dim='time_counter').plot(cmap='jet', vmin=0, vmax=0.06)









    Out[17]:





<matplotlib.collections.QuadMesh at 0x7fe7d977b350>

Small Scales (<6°)

Temporal standard deviation of the large field:



In [18]:

    
# Small-scale (<6°) dataset
(signal_xyt-signal_LS).std(dim='time_counter').plot(cmap='jet', vmin=0, vmax=0.06)









    Out[18]:





<matplotlib.collections.QuadMesh at 0x7fe7d6e99790>



In [ ]: