Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Study wave propagation from the cross spectra between Sentinel2 optical channels

As part of the ILA9 session of the OTC23 online training course.

This Jupyter notebook will use the SEAScope viewer and SEAScope Python bindings to retrieve inputs data

  • author: dr.fab, OceanDataLab

  • date: 2022-10-31

Study wave propagation from the cross spectra between Sentinel2 optical channels

Sentinel-2 data are great to estimate wave properties in cloud free conditions and in the presence of sun or sky glitter.

The different optical bands of Sentinel2 are sensing the sea surface at slightly different times since the different detectors for each band are placed at a slighly different place in the MSI camera focal plane, creating an optical parallax.

This notebook will teach you how to:

  • Import data from SEAScope

  • Investigate cross spectra between optical Sentinel2 optical channels

This notebook uses the following data collection:

  • directory: sentinel2, label: Ocean Color sentinel2 RGB

Check that you have them in your SEAScope catalogue before continuing with this notebook
These collections can be downloaded at https://ftp.odl.bzh/events/otc23/oct31.zip

To learn how to add data to SEAScope, you can watch one of the following YouTube videos:

import os
import numpy
from numpy.fft import fft2, fftshift, ifft2, fftfreq
import math
import matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from scipy import signal
import scipy.io as sp
from scipy.ndimage.filters import gaussian_filter
from scipy.interpolate import interp1d
from scipy.ndimage import rotate 
from numpy.lib.recfunctions import stack_arrays
from datetime import datetime,timedelta

Extract Sentinel2 data from SEAScope

locate the Sentinel2 data in the SEAScope viewer

  1. In the “Catalogue” (the panel on the right), select :

    • collection Ocean Color Sentinel2 RGB, variable rgb

  2. Select the date where you have data available (highlighted in white in the timeline)

    • The 12th of August 2018 should be good

  3. Locate the data on the globe

    • It should be in the Agulhas region

Steps for data extraction from SEAScope

  1. Zoom-in to properly see the waves (the resolution of the extraction depends on the zoom level)

  2. Draw a polygon over the region you want to extract:

    1. Start drawing by clicking on the 3rd button to the right of the pen icon in the toolbar or by hitting the P key

    2. Place the corners of the polygon by left-clicking on the globe

    3. Once all the corners placed, hit the Enter or Return key to commit your polygon

  3. Select the polygon you just created by clicking on it

  4. Extract the data intersecting that polygon by clicking on the extract button in the “Selected annotation” panel

# Load data directly from viewer memory
from SEAScope.lib import get_extracted_data
extractions = get_extracted_data()
for k, data in enumerate(extractions.keys()):
    print('{} - {}'.format(k, os.path.basename(data)))
    print('\n'.join(['\t{}'.format(x) for x in extractions[data]['data']]))

Play with the data

** Compute Cross spectra and compare the phase to the linear dispertion relation:
**

granule_uri =  next( v for i, v in enumerate(extractions.keys()) if i == 0)
extraction = extractions[granule_uri]
start = extraction['meta']['start']
print(extraction['meta']['fields'])

** Selection of the two channels you want to compare:
** B02_TOA_reflectance is the blue channel
B04_TOA_reflectance is the red channel

### Set name_var_b4 and name_var_b6
name_var_b4 = 'B04_TOA_reflectance'
name_var_b2 = 'B02_TOA_reflectance'
# Extract data
extraction = extractions[granule_uri]   
b4 = extraction['data'][name_var_b4]
b2 = extraction['data'][name_var_b2]

** Plot first channel **

fig = plt.figure(figsize=(16,10))
plt.imshow(numpy.flipud(b2),interpolation='bicubic',cmap='gray')
cbar=plt.colorbar()

Crop the image borders to avoid missing data in FFT calculations

** Define crop number (by default 40) to remove borders **

# define the numbers of pixels to crop from the border (will be applied to crop left, right, up, down)
crop = 40
b4_crop = b4[crop:-crop, crop:-crop]
b2_crop = b2[crop:-crop, crop:-crop]
fig = plt.figure(figsize=(16,10))
plt.imshow(numpy.flipud(b4_crop),interpolation='bicubic',cmap='gray')
cbar=plt.colorbar()
# Reconstruct coordinates and compute ground_spacing
from SEAScope.lib.utils import get_lonlat
lon2D, lat2D = get_lonlat(extraction, numpy.shape(b4))
ground_spacing = (lat2D[int(b4.shape[0]/2+1), int(b4.shape[1]/2)]
                  - lat2D[int(b4.shape[0]/2),int(b4.shape[1]/2)]) * 100000
print(ground_spacing)
def rebin(a, shape):
    sh = shape[0],a.shape[0]//shape[0],shape[1],a.shape[1]//shape[1]
    return a.reshape(sh).mean(-1).mean(1)
def looks2xspec(im1,im2,periodo_size):
    imshape = numpy.array(im1.shape, dtype='int32')

    ###########################################################################
    # Set periodograms/looks/specs sizes and positions
    ###########################################################################

    aziwindow = numpy.hanning(periodo_size+2)[1:-1]
    ranwindow = numpy.hanning(periodo_size+2)[1:-1]
    window = numpy.sqrt(numpy.outer(aziwindow, ranwindow))
    ###########################################################################
    # Compute periodograms and compute co/cross spectra
    ###########################################################################

    count = numpy.ceil(imshape[0]/periodo_size).astype('int32')
    perpos = (numpy.floor(numpy.linspace(0, imshape[0]-periodo_size, num=count)+0.5).astype('int32'),
              numpy.floor(numpy.linspace(0, imshape[1]-periodo_size, num=count)+0.5).astype('int32'))
    specshape = numpy.array((periodo_size, periodo_size), dtype='int32')
    specs = numpy.zeros(specshape, dtype='complex64')

    for appos in iter(perpos[0]):
        for rppos in iter(perpos[1]):
            sub1 = im1[appos:appos+periodo_size, rppos:rppos+periodo_size]
            sub1 = sub1 - numpy.mean(sub1)
            per1 = fftshift(fft2(sub1*window))/periodo_size
            sub2 = im2[appos:appos+periodo_size, rppos:rppos+periodo_size]
            sub2 = sub2 - numpy.mean(sub2)
            per2 = fftshift(fft2(sub2*window))/periodo_size
            specs += per1 * numpy.conj(per2)
            
    specs[int(periodo_size/2-2):int(periodo_size/2+3), int(periodo_size/2-2):int(periodo_size/2+3)] = 0
    specs[int(periodo_size/2-3):int(periodo_size/2+4), int(periodo_size/2)] = 0
    specs[int(periodo_size/2), int(periodo_size/2-3):int(periodo_size/2+4)] = 0

    return specs
# Define period to compute spectrum
periodo_size = 64
spec = looks2xspec(b4_crop, b2_crop,periodo_size)
kran = (numpy.arange(periodo_size)-periodo_size/2.)/periodo_size*2*numpy.pi/ground_spacing
kazi = (numpy.arange(periodo_size)-periodo_size/2.)/periodo_size*2*numpy.pi/ground_spacing
fig = plt.figure(figsize = (14,10))
plt.subplot(2, 2, 1)
plt.imshow(numpy.flipud(numpy.abs(spec)*1e5), extent=[kran[0], kran[-1],
            kazi[-1], kazi[1]], aspect='auto',cmap='jet')
plt.xlabel('wavenumber [rad/m]')
plt.ylabel('wavenumber [rad/m]')
plt.title('Absolute value')
plt.colorbar()
plt.subplot(2, 2, 2)
plt.imshow(numpy.flipud(numpy.real(spec)*1e5), extent=[kran[0], kran[-1],
            kazi[-1], kazi[1]],aspect='auto',cmap='jet')
plt.xlabel('wavenumber [rad/m]')
plt.ylabel('wavenumber [rad/m]')
plt.title('Real Part')
plt.colorbar()
plt.subplot(2, 2, 3)
plt.imshow(numpy.flipud(numpy.imag(spec)*1e5), extent=[kran[0], kran[-1],
            kazi[-1], kazi[1]],aspect='auto',cmap='jet')
plt.xlabel('wavenumber [rad/m]')
plt.ylabel('wavenumber [rad/m]')
plt.title('Imaginary part')
plt.colorbar()
plt.subplot(2, 2, 4)
plt.imshow(numpy.flipud(numpy.angle(spec)), extent=[kran[0], kran[-1],
           kazi[-1], kazi[1]],aspect='auto',cmap='jet')
plt.xlabel('wavenumber [rad/m]')
plt.ylabel('wavenumber [rad/m]')
plt.title('Phase')
plt.colorbar()

Rotate the spectrum to have the maximum energy in the horizontal direction

** Adapt the rotation angle until the spectum maximum energy is aligned with horizontal direction **

rotation_angle = 68#! parameter that may be changed

rot_spec_imag = rotate(numpy.imag(spec),rotation_angle,reshape=False)
rot_spec_real = rotate(numpy.real(spec),rotation_angle,reshape=False)
rot_spec = rot_spec_real + numpy.complex(0,1)*rot_spec_imag
size = numpy.shape(rot_spec)
fig = plt.figure(figsize = (10,10))
plt.imshow(numpy.flipud(numpy.abs(rot_spec)*1e5), extent=[kran[0], kran[-1],
            kazi[-1], kazi[1]], aspect='auto',cmap='jet')

estimate phase velocity

** Define time shift between the two bands (variable shift_time_bands).
If you are comparing red and blue channels shift_time_bands is around 0.8. More information can be found in the Sentinel 2 documentation.**

fig = plt.figure(figsize = (10,10))
shift_time_bands = 0.5

mini = int(periodo_size/2+4)
spec = looks2xspec(b4_crop, b2_crop,periodo_size)
vphase = numpy.angle(numpy.mean(rot_spec[:,mini:-5],axis=0)) / kran[mini:-5] / shift_time_bands
sp = numpy.mean(numpy.real(rot_spec[:,mini:-5]),axis=0)/numpy.mean(numpy.real(spec[:,mini:-5]))
plt.plot(kran[mini:-5],vphase, 'b', label='observed phase velocity')
plt.plot(kran[mini:-5],sp, 'g', label='observed spectrum')
plt.plot(kran[mini:-5],numpy.sqrt(9.81/kran[mini:-5]), 'r', label='theoretical phase velocity')
plt.plot(kran[mini:-5],-numpy.sqrt(9.81/kran[mini:-5]), 'r')

plt.ylim(-25,25)
plt.xlim(0.01,0.08)
plt.legend()
plt.xlabel('wavenumber')
plt.ylabel('m/s')
OVL Manual
Computing Sentinel-1 radial velocities
OVL Manual
Swell refraction