<jupyter_start><jupyter_text>Sea Surface Altimetry Data Analysis<img src="http://marine.copernicus.eu/documents/IMG/SEALEVEL_GLO_SLA_MAP_L4_REP_OBSERVATIONS_008_027.png"      width="15%"      align=left     alt="Globe">For this example we will use gridded [sea-surface altimetry data from The Copernicus Marine Environment](http://marine.copernicus.eu/services-portfolio/access-to-products/?option=com_csw&view=details&product_id=SEALEVEL_GLO_PHY_L4_REP_OBSERVATIONS_008_047):This is a widely used dataset in physical oceanography and climate.The dataset has already been extracted from copernicus and stored in google cloud storage in [xarray-zarr](http://xarray.pydata.org/en/latest/io.htmlzarr) format.<jupyter_code>%matplotlib inline

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import dask.array as dsa

plt.rcParams['figure.figsize'] = (15,10)<jupyter_output><empty_output><jupyter_text>Initialize DatasetHere we load the dataset from the zarr store. Note that this very large dataset initializes nearly instantly, and we can see the full list of variables and coordinates.<jupyter_code>import intake
cat = intake.Catalog("https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/ocean.yaml")
ds = cat["sea_surface_height"].to_dask()
ds<jupyter_output><empty_output><jupyter_text>Examine MetadataFor those unfamiliar with this dataset, the variable metadata is very helpful for understanding what the variables actually represent<jupyter_code>for v in ds.data_vars:
    print('{:>10}: {}'.format(v, ds[v].attrs['long_name']))<jupyter_output><empty_output><jupyter_text>Visually Examine Some of the DataLet's do a sanity check that the data looks reasonable:<jupyter_code>plt.rcParams['figure.figsize'] = (15, 8)
ds.sla.sel(time='2000-01-01', method='nearest').plot()<jupyter_output><empty_output><jupyter_text>Same thing using interactive graphics<jupyter_code>import holoviews as hv
from holoviews.operation.datashader import regrid
hv.extension('bokeh')
dataset = hv.Dataset(ds.sla)
hv_im = (dataset.to(hv.Image, ['longitude', 'latitude'], dynamic=True)
                .redim.range(sla=(-0.5, 0.5))
                .options(cmap='RdBu_r', width=800, height=450, colorbar=True))

%output holomap='scrubber' fps=2
regrid(hv_im, precompute=True)<jupyter_output><empty_output><jupyter_text>Create and Connect to Dask Distributed Cluster<jupyter_code>from dask.distributed import Client, progress

from dask_kubernetes import KubeCluster
cluster = KubeCluster()
cluster.adapt(minimum=4, maximum=40)
cluster<jupyter_output><empty_output><jupyter_text>** ☝️ Don't forget to click the link above to view the scheduler dashboard! **<jupyter_code>client = Client(cluster)
client<jupyter_output><empty_output><jupyter_text>Timeseries of Global Mean Sea LevelHere we make a simple yet fundamental calculation: the rate of increase of global mean sea level over the observational period.<jupyter_code># the number of GB involved in the reduction
ds.sla.nbytes/1e9
# the computationally intensive step
sla_timeseries = ds.sla.mean(dim=('latitude', 'longitude')).load()
sla_timeseries.plot(label='full data')
sla_timeseries.rolling(time=365, center=True).mean().plot(label='rolling annual mean')
plt.ylabel('Sea Level Anomaly [m]')
plt.title('Global Mean Sea Level')
plt.legend()
plt.grid()<jupyter_output><empty_output><jupyter_text>In order to understand how the sea level rise is distributed in latitude, we can make a sort of [Hovmöller diagram](https://en.wikipedia.org/wiki/Hovm%C3%B6ller_diagram).<jupyter_code>sla_hov = ds.sla.mean(dim='longitude').load()
fig, ax = plt.subplots(figsize=(12, 4))
sla_hov.name = 'Sea Level Anomaly [m]'
sla_hov.transpose().plot(vmax=0.2, ax=ax)<jupyter_output><empty_output><jupyter_text>We can see that most sea level rise is actually in the Southern Hemisphere.  Sea Level VariabilityWe can examine the natural variability in sea level by looking at its standard deviation in time.<jupyter_code>sla_std = ds.sla.std(dim='time').load()
sla_std.name = 'Sea Level Variability [m]'
ax = sla_std.plot()
_ = plt.title('Sea Level Variability')<jupyter_output><empty_output><jupyter_text>Spectral AnalysisThis is an advanced, research-grade example. Here we perform wavenumber-frequency spectral analysis of the SSH signal using methods similar to those described in [Abernathey & Wortham (2015)](https://journals.ametsoc.org/doi/10.1175/JPO-D-14-0160.1). Step 1: Extract a sector in the Eastern PacificThis sector is chosen because it has very little land.<jupyter_code>sector = ds.sla.sel(longitude=slice(180, 230), latitude=slice(-70, 55, 4))
sector_anom = (sector - sector.mean(dim='longitude'))
sector_anom
sector_anom[0].plot()<jupyter_output><empty_output><jupyter_text>Step 2: Rechunk, reshape, and window the data for efficient to prepare for FFT calculation<jupyter_code># reshape data into arrays 365 days long and rechunk
nsegments = 24
segment_len = 365
sector_reshape = (sector_anom.isel(time=slice(0, nsegments*segment_len))
                             .transpose('latitude', 'time', 'longitude')
                             .chunk({'time': segment_len}))
sector_reshape
# now get the raw dask array
data = sector_reshape.data

arrays = [data[:, n*segment_len:(n + 1)*segment_len][np.newaxis]
          for n in range(nsegments)]
stacked = dsa.concatenate(arrays)
stacked
# apply windows
data_windowed = (stacked
                 * np.hanning(stacked.shape[-1])[None, None, None, :]
                 * np.hanning(stacked.shape[-2])[None, None, :, None])<jupyter_output><empty_output><jupyter_text>Step 3: Actually calculate the Fourier transform and power spectral density<jupyter_code># take FFT
data_fft = dsa.fft.fftn(data_windowed, axes=(-2, -1))

# take power spec and average over segments
power_spec = np.real(data_fft * np.conj(data_fft)).mean(axis=0)
power_spec
# do the computation and load results into memory
power_spec_shift = np.fft.fftshift(power_spec.compute(), axes=(-2, -1))<jupyter_output><empty_output><jupyter_text>Step 4: Define spectral coordinates and put everything back together in an DataArray<jupyter_code>freq = np.fft.fftshift(np.fft.fftfreq(segment_len))

# wavelength is a bit trickier because it depends on latitude
R = 6.37e6
# in km
dx = np.deg2rad(0.25) * R * np.cos(np.deg2rad(sector.latitude)) / 1000
inv_wavelength = np.vstack([np.fft.fftshift(np.fft.fftfreq(len(sector.longitude), d))
                            for d in dx.values])

ps_da = xr.DataArray(power_spec_shift, dims=('latitude', 'freq', 'wavenumber'),
                     coords={'latitude': sector.latitude,
                             'freq': ('freq', -freq, {'units': r'days$^{-1}$'}),
                              'inverse_wavelength': (('latitude', 'wavenumber'),
                                                     inv_wavelength, {'units': r'km$^{-1}$'})},
                     name='SSH_power_spectral_density')
ps_da<jupyter_output><empty_output><jupyter_text>Step 5: Plot wavenumber-frequency power spectra at different latitudes<jupyter_code>from matplotlib.colors import LogNorm

for lat in range(-55, 55, 10):
    plt.figure()
    (ps_da.sel(latitude=lat, method='nearest')
          .swap_dims({'wavenumber': 'inverse_wavelength'})
          .transpose().plot(norm=LogNorm()))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11844978165938864, 'fraction_numerical': 0.017876637554585153, 'mean_word_length': 5.112593828190159, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 63, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '37489', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2558, 'n_tokens_neox': 2409, 'n_words': 583, 'path': 'notebooks/ocean_ssh_example.ipynb', 'repo_name': 'valentina-s/Oceans19-data-science-tutorial', 'size': '7328'}