Wave measurements from SAR

This notebook is from Salman Khan and Emilio Echevarria, there is a great presentation about how they use synthetic aperture radar here.

About Sentinel-1 Ocean Surface Wave Measurements

Note

Sentinel-1 A and B satellites are part of Europe’s Copernicus programme and were launched in April 2014 and 2016, respectively.

The satellites fly in the same polar orbital plane, but 180° out of orbital phase, and offer a 6-day repeat cycle at the equator and more frequently at higher latitudes.

Both satellites are equipped with identical C-band SAR instruments, which can measure the directional wind-wave spectra for wavelengths longer than ~150m in azimuth direction. In the open ocean, they operate in wave mode (WV), measuring alternately in near-range with an incidence angle of ~23.5° (WV1) and far-range with an incidence angle of ~36.5° (WV2), and acquire 20 km x 20 km vignettes every 100 km along-track. Since they alternate between WV1 and WV2, consecutive vignettes for each incidence angle are 200 km apart.

In this notebook, we will use Sentinel-1 A and B SAR measurements over the Australasian region. The data is accesible through the AODN portal, or Thredds server.

The amount of data for each sattelite (A and B) for different wave modes (WV1 and WV2) and orbit’s direction (ascending or descending passes), allocated in 2.5°x2.5° bins looks like this:

Data acquisition instruments

Imported Python packages

  • cartopy = 0.16.0

  • matplotlib = 2.2.2

  • numpy = 1.15.0

  • xarray = 0.11.3

  • wavespectra = 3.5.3

  • utility functions from utilfuncs.py

import os
import sys
import warnings

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
from wavespectra import SpecArray
import xarray as xr

from utilfuncs import k2f, ek2f, s1_compute_efth, cmap_wavespectra, plot_part

warnings.filterwarnings('ignore')

%matplotlib inline
%config InlineBackend.figure_format = 'svg'
plt.rcParams['mathtext.fontset'] = 'cm'

Daily along-track measurements

We will first read a daily along-track netCDF file and visualise selected measurements off Albany, WA on a map:

# Read and select measurements near South Australia
url = 'http://thredds.aodn.org.au/thredds/dodsC/IMOS/SRS/Surface-Waves/SAR/DELAYED/SENTINEL-1A/2019/08/IMOS_SRS-Surface-Waves_W_20190801T040752Z_Sentinel-1A_FV01_DM00_K2_END-20190801T220024Z.nc'+'#fillmismatch'

dset_s1 = xr.open_dataset(url, decode_times =  True, decode_cf = True, mask_and_scale =  True)

# Selecting the measurement close to Albany
dset_WA = dset_s1.isel(TIME = slice(384, 404))

# Plot measurement locations
fig = plt.figure(figsize = (10, 5), dpi = 100)
ax = plt.axes(projection = ccrs.PlateCarree(central_longitude = 195))
ax.coastlines()
ax.plot(dset_s1.LONGITUDE, dset_s1.LATITUDE, '.', 
        transform = ccrs.PlateCarree(), markersize = 1)

ax.plot(dset_WA.LONGITUDE, dset_WA.LATITUDE,
        'r.', transform = ccrs.PlateCarree(), markersize = 3)
plt.show()
../_images/sar_waves_7_0.svg

Auxialliary variables

We can plot some auxialliary variables such as the satellite measures along-track:

  • satellite heading HEADING

  • incidence angle INC_ANGLE

  • distance to coast DIST2COAST, and

  • depth BOT_DEPTH

To check the variables available we can do:

dset_WA
<xarray.Dataset>
Dimensions:               (DIRECTION: 72, PARTITION: 5, TIME: 20, WAVNUM: 60)
Coordinates:
  * DIRECTION             (DIRECTION) float32 0.0 5.0 10.0 ... 345.0 350.0 355.0
  * PARTITION             (PARTITION) int8 0 1 2 3 4
  * WAVNUM                (WAVNUM) float32 0.005236 0.005574 ... 0.1967 0.2094
    LATITUDE              (TIME) float32 -53.57 -52.14 -51.87 ... -38.07 -36.74
    LONGITUDE             (TIME) float32 124.2 126.3 123.3 ... 119.7 117.3 119.1
  * TIME                  (TIME) datetime64[ns] 2019-08-01T10:45:42.000000256...
Data variables: (12/17)
    AZ_CUTOFF             (TIME) float32 344.0 378.0 304.0 ... 220.0 210.0 236.0
    HEADING               (TIME) float32 341.0 339.1 341.7 ... 343.5 345.1 343.9
    INC_ANGLE             (TIME) float32 22.89 36.0 22.94 ... 36.36 23.29 36.41
    NRCS                  (TIME) float32 -1.507 -10.63 -2.39 ... -7.262 -16.84
    SNR                   (TIME) float32 30.75 27.38 29.15 ... 25.23 31.75 25.9
    WDIR_ECMWF            (TIME) float32 288.8 307.4 307.0 ... 358.5 1.01 40.07
    ...                    ...
    DP_PART               (TIME, PARTITION) float32 93.5 315.6 nan ... nan nan
    EKTH                  (TIME, DIRECTION, WAVNUM) float32 362.2 419.6 ... 0.0
    EKTH_PART             (TIME, DIRECTION, WAVNUM) float32 1.0 1.0 ... nan nan
    EKTH_quality_control  (TIME, DIRECTION, WAVNUM) float32 2.0 2.0 ... nan nan
    HS_PART               (TIME, PARTITION) float32 3.398 2.14 nan ... nan nan
    WP_PART               (TIME, PARTITION) float32 269.2 561.9 nan ... nan nan
Attributes: (12/30)
    Conventions:                     CF-1.6,IMOS-1.4
    abstract:                        Sentinel-1 A and B satellites are part o...
    acknowledgement:                 Any users of IMOS data are required to c...
    author:                          Khan, Salman Saeed
    author_email:                    salmansaeed.khan@csiro.au
    citation:                        The citation in a list of references is:...
    ...                              ...
    source:                          Copernicus Australasia
    standard_name_vocabulary:        NetCDF Climate and Forecast (CF) Metadat...
    time_coverage_end:               2019-08-01T22:00:24Z
    time_coverage_start:             2019-08-01T04:07:52Z
    title:                           Sentinel-1 A C-SAR surface wave spectra
    DODS_EXTRA.Unlimited_Dimension:  TIME

Let’s plot the aforementionned variables

fig = plt.figure(figsize = (7.5, 3))
ax = plt.subplot(111)
dset_WA.HEADING.plot(ax = ax, lw=2, c='r')
plt.title('Heading') 
plt.grid(linestyle = 'dashed')

fig = plt.figure(figsize = (7.5, 3))
ax = plt.subplot(111)
dset_WA.INC_ANGLE.plot(ax = ax, lw=2, c='b')
plt.title('Incidence Angle') 
plt.grid(linestyle = 'dashed')

fig = plt.figure(figsize = (7.5, 3))
ax = plt.subplot(111)
dset_WA.DIST2COAST.plot(ax = ax, lw=2, c='k')
plt.title('Distance to coast')
plt.grid(linestyle = 'dashed')

fig = plt.figure(figsize = (7.5, 3))
ax = plt.subplot(111)
dset_WA.BOT_DEPTH.plot(ax = ax, lw=2, c='g')
plt.title('Depth')
plt.grid(linestyle = 'dashed')
../_images/sar_waves_11_0.svg../_images/sar_waves_11_1.svg../_images/sar_waves_11_2.svg../_images/sar_waves_11_3.svg

Frequency spectra conversion

We will now convert directional wavenumber spectra to frequency spectra using deep water assumptions: efth - spectral variance density - units m2.s.deg-1

# compute efth
dset_WA_ws = dset_WA.assign(efth = s1_compute_efth(dset_WA))

Plot some wave statistics derived from directional frequency spectra:

  • swell significant height

  • mean period - Tm02

  • mean wave direction

Note

A directional wave spectrum describes how the wave energy (more specifically, the variance spectral density) is distributed over frequencies (or wavenumbers) and directions. For a detailed description of wave spectra, its uses and applications, and the transformations from wavenumber to frequency spectra, see: Holthuijsen, L. (2007). Waves in Oceanic and Coastal Waters. Cambridge: Cambridge University Press. https://doi.org/10.1017/CBO9780511618536

It is very common to describe the wave climate with integrated (also called bulk) wave parameters, instead of with a directional wave spectra. These parameters represent average conditions of the wave climate.

# plot stats
# hs
fig = plt.figure(figsize = (7.5, 3))
ax = plt.subplot(111)
dset_WA_ws.spec.hs(tail = True).plot(ax = ax)
plt.title('Swell significant wave height')
plt.grid(linestyle = 'dashed')

# tm02
fig = plt.figure(figsize = (7.5, 3))
ax = plt.subplot(111)
dset_WA_ws.spec.tm02().plot(ax = ax)
plt.title('Mean period - Tm02')
plt.grid(linestyle = 'dashed')

# dm
fig = plt.figure(figsize = (7.5, 3))
ax = plt.subplot(111)
dset_WA_ws.spec.dm().plot(ax = ax)
plt.ylabel('Mean wave to direction $(^\circ)$')
plt.title('Mean wave direction')
plt.grid(linestyle = 'dashed')
plt.show()
../_images/sar_waves_15_0.svg../_images/sar_waves_15_1.svg../_images/sar_waves_15_2.svg

Visualise directional frequency spectra:

  • 1st directional spectra

  • all directional spectra

  • all omni-directional spectra

# plot 1st directional spectra
cmap = cmap_wavespectra(n = 64)
qcs = (dset_WA_ws
       .isel(time = 0)
       .spec
       .plot
       .contourf(cmap = cmap,
                 as_period = True, 
                 as_log10 = False, 
                 show_direction_label = True,
                 levels = cmap.N,
                 figsize = (7.5, 7.5),
                 robuts = True))
angle_labels = ['N', 'NE', 'E', 'SE', 'S', 'SW', '','NW']
qcs.ax.set_rlabel_position(180 + 45)
qcs.ax.set_ylim(0, 30)#set_rmax(30)
qcs.ax.set_thetagrids(angles = range(0, 360, 45),
                      labels = angle_labels)
(<a list of 16 Line2D ticklines objects>,
 [Text(0.0, 0, 'N'),
  Text(0.7853981633974483, 0, 'NE'),
  Text(1.5707963267948966, 0, 'E'),
  Text(2.356194490192345, 0, 'SE'),
  Text(3.141592653589793, 0, 'S'),
  Text(3.9269908169872414, 0, 'SW'),
  Text(4.71238898038469, 0, ''),
  Text(5.497787143782138, 0, 'NW')])
../_images/sar_waves_17_1.svg
## plot all directional spectra
ncols = 4
fg = (dset_WA_ws
      .efth
      .spec
      .plot
      .contourf(col = "time", col_wrap = ncols, cmap = cmap,
                as_period = True, as_log10 = False, 
                show_direction_label = True, levels = cmap.N,
                robuts = True, size = 3, aspect = 1))
angle_labels = ['N', 'NE', 'E', 'SE', 'S', 'SW', 'W','NW']        
for i, ax in enumerate(fg.axes.flat):
    ax.set_rlabel_position(180 + 45)
    ax.set_ylim(0, 30)#set_rmax(30)
    ax.set_thetagrids(angles = range(0, 360, 45),
                      labels = angle_labels)
fg.set_titles(pad = 15)
fg.set_ylabels('Wave period [s]', labelpad = 20)
fg.set_xlabels('Wave to direction [$^\circ$]', labelpad = 0)
fg.cbar.set_label('Variance density [$m^2 s deg{-1}$]', labelpad = 10)
Error in callback <function flush_figures at 0x7f7ff6b0c8c8> (for post_execute):
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/ipykernel/pylab/backend_inline.py in flush_figures()
    118         # ignore the tracking, just draw and close all figures
    119         try:
--> 120             return show(True)
    121         except Exception as e:
    122             # safely show traceback if in IPython, else raise

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/ipykernel/pylab/backend_inline.py in show(close, block)
     40             display(
     41                 figure_manager.canvas.figure,
---> 42                 metadata=_fetch_figure_metadata(figure_manager.canvas.figure)
     43             )
     44     finally:

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/IPython/core/display.py in display(include, exclude, metadata, transient, display_id, *objs, **kwargs)
    311             publish_display_data(data=obj, metadata=metadata, **kwargs)
    312         else:
--> 313             format_dict, md_dict = format(obj, include=include, exclude=exclude)
    314             if not format_dict:
    315                 # nothing to display (e.g. _ipython_display_ took over)

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/IPython/core/formatters.py in format(self, obj, include, exclude)
    178             md = None
    179             try:
--> 180                 data = formatter(obj)
    181             except:
    182                 # FIXME: log the exception

<decorator-gen-9> in __call__(self, obj)

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/IPython/core/formatters.py in catch_format_error(method, self, *args, **kwargs)
    222     """show traceback on failed format call"""
    223     try:
--> 224         r = method(self, *args, **kwargs)
    225     except NotImplementedError:
    226         # don't warn on NotImplementedErrors

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/IPython/core/formatters.py in __call__(self, obj)
    339                 pass
    340             else:
--> 341                 return printer(obj)
    342             # Finally look for special method names
    343             method = get_real_method(obj, self.print_method)

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/IPython/core/pylabtools.py in <lambda>(fig)
    252         jpg_formatter.for_type(Figure, lambda fig: print_figure(fig, 'jpg', **kwargs))
    253     if 'svg' in formats:
--> 254         svg_formatter.for_type(Figure, lambda fig: print_figure(fig, 'svg', **kwargs))
    255     if 'pdf' in formats:
    256         pdf_formatter.for_type(Figure, lambda fig: print_figure(fig, 'pdf', **kwargs))

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, **kwargs)
    130         FigureCanvasBase(fig)
    131 
--> 132     fig.canvas.print_figure(bytes_io, **kw)
    133     data = bytes_io.getvalue()
    134     if fmt == 'svg':

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2194 
   2195                     bbox_inches = self.figure.get_tightbbox(
-> 2196                         renderer, bbox_extra_artists=bbox_extra_artists)
   2197                     if pad_inches is None:
   2198                         pad_inches = rcParams['savefig.pad_inches']

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/matplotlib/figure.py in get_tightbbox(self, renderer, bbox_extra_artists)
   2504 
   2505         for a in artists:
-> 2506             bbox = a.get_tightbbox(renderer)
   2507             if bbox is not None and (bbox.width != 0 or bbox.height != 0):
   2508                 bb.append(bbox)

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/matplotlib/artist.py in get_tightbbox(self, renderer)
    284             if clip_path is not None and bbox is not None:
    285                 clip_path = clip_path.get_fully_transformed_path()
--> 286                 bbox = Bbox.intersection(bbox, clip_path.get_extents())
    287         return bbox
    288 

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/matplotlib/path.py in get_extents(self, transform, **kwargs)
    597             for curve, code in self.iter_bezier(**kwargs):
    598                 # places where the derivative is zero can be extrema
--> 599                 _, dzeros = curve.axis_aligned_extrema()
    600                 # as can the ends of the curve
    601                 xys.append(curve([0, *dzeros, 1]))

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/matplotlib/bezier.py in axis_aligned_extrema(self)
    296         roots = []
    297         for i, pi in enumerate(dCj.T):
--> 298             r = np.roots(pi[::-1])
    299             roots.append(r)
    300             dims.append(np.full_like(r, i))

<__array_function__ internals> in roots(*args, **kwargs)

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/lib/polynomial.py in roots(p)
    241         A = diag(NX.ones((N-2,), p.dtype), -1)
    242         A[0,:] = -p[1:] / p[0]
--> 243         roots = eigvals(A)
    244     else:
    245         roots = NX.array([])

<__array_function__ internals> in eigvals(*args, **kwargs)

/usr/share/miniconda/envs/coast/lib/python3.6/site-packages/numpy/linalg/linalg.py in eigvals(a)
   1067         _raise_linalgerror_eigenvalues_nonconvergence)
   1068     signature = 'D->D' if isComplexType(t) else 'd->D'
-> 1069     w = _umath_linalg.eigvals(a, signature=signature, extobj=extobj)
   1070 
   1071     if not isComplexType(t):

KeyboardInterrupt: 
# plot all omni-directional spectra
ncols = 4
fg = dset_WA_ws.spec.oned().plot.line(col = 'time', col_wrap = ncols)
fg.set_axis_labels(y_var = 'var. spec. density ($m^2s/deg$)')
for i, ax in enumerate(fg.axes.flat):
    ax.grid(True)
../_images/sar_waves_19_0.svg

Visualise directional frequency spectra swell partitions:

  • paritions of first directional spectrum

  • partitions of all directional spectra

# plot partitions for 1st 2d spectrum
fig = plot_part(dset_WA_ws.EKTH_PART.isel(TIME = 0))
../_images/sar_waves_21_0.svg
# plot partitions for all 2d spectra
fig = plot_part(dset_WA_ws.EKTH_PART)
../_images/sar_waves_22_0.svg

Plot some primary swell wave statistics:

# plot primary swell stats
# swell height
fig = plt.figure(figsize = (7.5, 3))
ax = plt.subplot(111)
dset_WA_ws.HS_PART.isel(PARTITION = 0).plot()
plt.title('Primary swell significant height') 
plt.grid(linestyle = 'dashed')

# direction at spectral peak
fig = plt.figure(figsize = (7.5, 3))
ax = plt.subplot(111)
dset_WA_ws.DP_PART.isel(PARTITION = 0).plot()
plt.title('Primary swell peak direction') 
plt.grid(linestyle = 'dashed')

# wavelength at spectral peak
fig = plt.figure(figsize = (7.5, 3))
ax = plt.subplot(111)
dset_WA_ws.WP_PART.isel(PARTITION = 0).plot()
plt.title('Primary swell peak wavelength') 
plt.grid(linestyle = 'dashed')
../_images/sar_waves_24_0.svg../_images/sar_waves_24_1.svg../_images/sar_waves_24_2.svg