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

author: Fabrice COLLARD, OceanDataLab
date: 2024-09-10

To install SEAScope, follow the instructions on the SEAScope website and demonstrated in our “SEAScope How To” YouTube playlist.

To install the required Python packages, retrieve the file requirement.txt. Create an environment:

python -m venv seascopenv

For Linux / MacOS users:

source seascopenv/bin/activate

pip install -r requirements.txt

For Windows:

source seascopenv\bin\activate.csh

pip install -r requirements.txt

compute radial velocities from Sentinel-1 Doppler Shift¶

The radial velocities are the ocean surface motions detected via the Doppler centroid shift. This Doppler shift is composed of several components:

The electronic misspointing induced Doppler
The plateform attitude induced Doppler
The wave induced Doppler (CDOP)
The current induced Doppler, or horizontal radial velocities Vh

First the three first components of the Doppler shift must be estimated and compensated before accessiong the surface current induced Doppler.

The relation between current induced Doppler shift $\Delta$ F and horizontal radial velocities $Vh$ is :

Vh = -\frac{\pi \Delta F}{ke sin(\theta)} = -\frac{\lambda \Delta F}{2 sin(\theta)}

(1)

with $ke$ the radar wavenumber, and $\theta$ the incidence angle.

In this notebook, we will use the Sentinel1 Level2 calibrated RVL product call RVC to estimate the horizontal radial velocities.

This notebook will teach you how to:

Import data from SEAScope
Export data to SEAScope
Derive horizontal radial velocity from Sentinel1 Doppler shift

This notebook uses the following data collection:

directory: sentinel1_l2_rvc, label: Sentinel1 L2 RVC

Check that you have it in your SEAScope catalogue before continuing with this notebook.
This collection was included in the data package you downloaded (the data_space_week_2024.zip).

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

How to add data to SEAScope on Windows
- or the “Add data to SEAScope” chapter of How to add data to SEAScope on Windows
How to add data to SEAScope on macOS
How to add data to SEAScope on Linux

# Import necessary modules
%matplotlib inline
import matplotlib.pyplot as plt
import numpy
import numpy as np
import os

# Create once the collection for this notebook
collection_id = None
import SEAScope.lib.utils
SEAScope.lib.utils.init_ids(20, 2000)

def cdop_func(x):
    """
    """
    return 1./(1.+np.exp(-x))


def cdop(u10, phi, inc, pol):
    '''
        Input
        -----------
        u10 : float, numpy.array
              wind speed in m/s
        phi : float, numpy.array
              wind direction relative to SAR look direction
        inc : float, numpy.array
              SAR incidence angle
        pol : string
              SAR polarization (VV or HH)
        Return
        -----------
        dop : numpy.array
              Estimated Doppler shift [Hz]
    '''
    # Check inputs
    sizes = np.array([np.size(inc), np.size(u10), np.size(phi)])
    size = sizes.max()
    if ((sizes != size) & (sizes != 1)).any():
        raise Exception('Inputs sizes do not agree.')
    if pol.upper() not in ['VV', 'HH']:
        raise Exception('Unknown polarisation : '+pol)
    # NN coefficients (W=weights and B=biases)
    # (coefficient names in mouche2012 are given)
    if pol.upper() == 'VV':
        # lambda[0:2,1]
        B1 = np.array([-0.343935744939, 0.108823529412, 0.15],
                      dtype='float32')
        # lambda[0:2,0]
        W1 = np.array([0.028213254683, 0.0411764705882, .00388888888889],
                      dtype='float32')
        # omega[i,0]
        B2 = np.array([14.5077150927, -11.4312028555, 1.28692747109,
                       -1.19498666071, 1.778908726, 11.8880215573,
                       1.70176062351, 24.7941267067, -8.18756617111,
                       1.32555779345, -9.06560116738],
                      dtype='float32')
        # omega[i,[3,2,1]]
        W2 = np.array([[19.7873046673, 22.2237414308, 1.27887019276],
                       [2.910815875, -3.63395681095, 16.4242081101],
                       [1.03269004609, 0.403986575614, 0.325018607578],
                       [3.17100261168, 4.47461213024, 0.969975702316],
                       [-3.80611082432, -6.91334859293, -0.0162650756459],
                       [4.09854466913, -1.64290475596, -13.4031862615],
                       [0.484338480824, -1.30503436654, -6.04613303002],
                       [-11.1000239122, 15.993470129, 23.2186869807],
                       [-0.577883159569, 0.801977535733, 6.13874672206],
                       [0.61008842868, -0.5009830671, -4.42736737765],
                       [-1.94654022702, 1.31351068862, 8.94943709074]],
                      dtype='float32')
        # gamma[0]
        B3 = np.array(4.07777876994, dtype='float32')
        # gamma[1:11]
        W3 = np.array([7.34881153553, 0.487879873912, -22.167664703,
                       7.01176085914, 3.57021820094, -7.05653415486,
                       -8.82147148713, 5.35079872715, 93.627037987,
                       13.9420969201, -34.4032326496],
                      dtype='float32')
        # beta
        B4 = np.array(-52.2644487109, dtype='float32')
        # alpha
        W4 = np.array(111.528184073, dtype='float32')
    elif pol.upper() == 'HH':
        # lambda[0:2,1]
        B1 = np.array([-0.342097701547, 0.118181818182, 0.15],
                      dtype='float32')
        # lambda[0:2,0]
        W1 = np.array([0.0281843837385, 0.0318181818182, 0.00388888888889],
                      dtype='float32')
        # omega[i,0]
        B2 = np.array([1.30653883096, -2.77086154074, 10.6792861882,
                       -4.0429666906, -0.172201666743, 20.4895916824,
                       28.2856865516, -3.60143441597, -3.53935574111,
                       -2.11695768022, -2.57805898849],
                      dtype='float32')
        # omega[i,[3,2,1]]
        W2 = np.array([[-2.61087309812, -0.973599180956, -9.07176856257],
                       [-0.246776181361, 0.586523978839, -0.594867645776],
                       [17.9261562541, 12.9439063319, 16.9815377306],
                       [0.595882115891, 6.20098098757, -9.20238868219],
                       [-0.993509213443, 0.301856868548, -4.12397246171],
                       [15.0224985357, 17.643307099, 8.57886720397],
                       [13.1833641617, 20.6983195925, -15.1439734434],
                       [0.656338134446, 5.79854593024, -9.9811757434],
                       [0.122736690257, -5.67640781126, 11.9861607453],
                       [0.691577162612, 5.95289490539, -16.0530462],
                       [1.2664066483, 0.151056851685, 7.93435940581]],
                      dtype='float32')
        # gamma[0]
        B3 = np.array(2.68352095337, dtype='float32')
        # gamma[1:11]
        W3 = np.array([-8.21498722494, -94.9645431048, -17.7727420108,
                       -63.3536337981, 39.2450482271, -6.15275352542,
                       16.5337543167, 90.1967379935, -1.11346786284,
                       -17.57689699, 8.20219395141],
                      dtype='float32')
        # beta
        B4 = np.array(-66.9554922921, dtype='float32')
        # alpha
        W4 = np.array(136.216953823, dtype='float32')
    # Make inputs as a single matrix (and clip phi in [0,180])
    inputs = np.zeros((3, size), dtype='float32')
    for ivar, var in enumerate((inc, u10, phi)):
        if sizes[ivar] == 1:
            inputs[ivar, :] = np.repeat(var, size)
        else:
            inputs[ivar, :] = np.ravel(var)
        if ivar == 2:
            inputs[ivar, :] = np.abs(((inputs[ivar, :]+180) % 360)-180)
        inputs[ivar, :] *= W1[ivar]
        inputs[ivar, :] += B1[ivar]
    # Compute CDOP
    B2 = np.tile(B2.reshape((11, 1)), (1, size))
    dop = W4*cdop_func(np.dot(W3, cdop_func(np.dot(W2, inputs) + B2)) + B3) + B4
    # Reshape output
    # (using the shape of input which have the maximum ndim)
    ndims = np.array([np.ndim(inc), np.ndim(u10), np.ndim(phi)])
    tmp = np.where(sizes == size)[0]
    ivar = tmp[ndims[tmp].argmax()]
    shp = np.shape((inc, u10, phi)[ivar])
    dop = dop.reshape(shp)
    return dop

1. Find the data of interest in SEAScope¶

Select and locate the data in the SEAScope viewer¶

In the “Catalogue” (the panel on the right), select collection Sentinel1 l2 RVC, variables:
- doppler_centroid_anomaly
- incidence_angle
- radial_direction
- wind_speed
- wind_direction
Pick a date where you have data available (highlighted in green in the timeline)
- The 25^th of July 2023 should be good
Zoom on the dataset

2. Get data from the SEAScope viewer using the python bindings¶

Extract the data from SEAScope¶

Zoom-in to properly see the structures in the data (the resolution of the extraction depends on the zoom level)
Draw a polygon over the region you want to extract:
1. Start drawing by clicking on the 3^rd 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
Select the polygon you just created by clicking on it
Extract the data intersecting that polygon by clicking on the extract button in the “Selected annotation” panel

You can also watch the How to export data from SEAScope video on YouTube.

Load the extracted data¶

#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']]))

Identify Doppler anomaly, incidence angle, radial direction and stokes drift in extracted granules¶

kt = 0 # ⚠️ You can adapt if you have several granules imported
granule_uri =  next( v for i, v in enumerate(extractions.keys()) if i == kt)
extraction = extractions[granule_uri]
start = extraction['meta']['start']
print(extraction['meta']['fields'])

Select the granules to analyze¶

# Extract data
dca = extraction['data']['doppler_centroid_anomaly']
inc = extraction['data']['incidence_angle']
radir = extraction['data']['radial_direction']
u10 = extraction['data']['wind_speed']
wdir = extraction['data']['wind_direction']

# Spatial coverage will be the mesh of the polygon used to perform the extraction
# It is available in the "gcps" metadata for each extraction
gcps = extraction['meta']['gcps']

# Make sure the result is reset
granule_id = None

3. Get coordinates and plot Doppler anomaly¶

# Get longitudes and latitudes from gcps
lon = [x['lon'] for x in gcps]
lat = [x['lat'] for x in gcps]

# Mask invalid dca values
dca[abs(dca) >100] = numpy.nan

from SEAScope.lib.utils import get_lonlat
lon2D, lat2D = get_lonlat(extraction, numpy.shape(dca))

plt.pcolormesh(lon2D, lat2D, dca, cmap='jet')
plt.xlabel('longitude (degree)')
plt.ylabel('latitude (degree)')
plt.colorbar() 
plt.title('Sentinel-1 Doppler centroid anomaly')

Compute Sea state Doppler¶

ssd = cdop(u10,wdir-radir,inc,'VV')

plt.pcolormesh(lon2D, lat2D, ssd, cmap='jet')
plt.xlabel('longitude (degree)')
plt.ylabel('latitude (degree)')
plt.colorbar() 
plt.title('Sea State Doppler')

plt.pcolormesh(lon2D, lat2D, dca-ssd, cmap='jet',vmin=-40,vmax=40)
plt.xlabel('longitude (degree)')
plt.ylabel('latitude (degree)')
plt.colorbar() 
plt.title('Sentinel-1 surface current Doppler (Hz)')

Compute Radial surface current¶

lambda_e = 0.05546576  # Sentinel-1 radar wavelength
vrad = -lambda_e*(dca-ssd)/2/numpy.sin(numpy.deg2rad(inc))

plt.pcolormesh(lon2D, lat2D, vrad, cmap='jet',vmin=-2,vmax=2)
plt.xlabel('longitude (degree)')
plt.ylabel('latitude (degree)')
plt.colorbar() 
plt.title('Sentinel-1 radial surface current (m/s)')

5. Export data back to SEAScope¶

import datetime

from SEAScope.lib.utils import create_collection, create_variable
from SEAScope.lib.utils import create_granule, set_field

Define Collection:
A granule has to belong to a collection
You may choose a name for the collection collection_name

collection_name = 'User - S-1 Radial surface current' # ⚠️ You can choose your collection name
# Create collection
collection_id, collection = create_collection(collection_name)

Define the granule that will hold the results. Granules are defined by their spatial and temporal coverage

# Reuse the temporal coverage of input granule
# - Start time of granule 
dtstart = extraction['meta']['start']
# - End time of granule 
dtstop = extraction['meta']['stop']

# Create granule
if granule_id is None:
    granule_id, granule = create_granule(collection_id, gcps, dtstart, dtstop)

The granule will contain each result matrix as what we call a “field”
You may choose a name for field_name

# Attach data fields to the granule
field_name = 'radial surface current'
set_field(granule, field_name, numpy.abs(vrad))

Next you need to tell SEAScope what to display: this is done by defining a “variable” which can be composed of one, two or three fields.
When a variable has two fields, these are considered to be the meridional and zonal components of a vector field.
You may choose a name for : variable_name

# Create a variable for the horizontal component
variable_name = f'swath {kt+1}'  # ⚠️ You can choose your label name
variable = create_variable(collection,variable_name, [field_name])

It is also possible to tell SEAScope how to display variables: the following cell shows how to set some rendering options

# Rendering configuration for the horizontal component
rcfg = variable['rendering']
rcfg['min'] = 0
rcfg['max'] = 2
rcfg['colormap'] = 'jet'
rcfg['zindex'] = 0.3
rcfg['filterMode'] = 'BILINEAR'

Finally, send the collection, granule and variable to SEAScope

hostname = '127.0.0.1'
port = 11155
with SEAScope.upload.connect(hostname, port) as link:
    # Upload the collection first (remember: a granule cannot exist before its collection has been uploaded!)
    SEAScope.upload.collection(link, collection)

    # Then the granule
    SEAScope.upload.granule(link, granule)

    # Now SEAscope has access to the data, all that remains is telling the application what and how it should display them, by uploading the variable definition
    SEAScope.upload.variable(link, variable)
    
    SEAScope.upload.rendering_config(link, rcfg)

6. Visualize your radial surface current in the SEAScope viewer¶

Your new collection should have appeared in the “Catalogue”

Your collection is the collection_name you have defined previously
Your variables are the variable_name you have defined previously
Select the CMEMS sea level NRT observations collection, variable geostrophic speed from MADT (streamlines)
Select the data you have just exported

You can still change the rendering configuration using the notebook, just remember to send the modified configuration to SEAScope
You can play with min, max values, the colormap ... and check how the rendering changes on SEAScope