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 current interaction in the Agulhas current region

As part of the ILB9 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 current interaction in the Agulhas current region

Wave current interactions are well known to bend wave propagation rays.

We will simulate the wave current interactions and compare it to atimeter and SAR wave observations

This notebook will teach you how to:

  • Import surface current data from SEAScope

  • Simulate wave current interactions

  • Export simulated wave rays in SEAScope

This notebook uses the following data collection:

  • directory: woc-bfnqg-3h, label: WOC BFN QG current

  • directory: sentinel1_sar_roughness, label: Sentinel1 SAR

  • directory: jason2_swh, label: Jason2 GDR

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 matplotlib
%matplotlib inline
import matplotlib.pyplot as plt
import numpy
import os
import datetime
import scipy
from numpy.random import randn
import pyproj


import numpy as np
from netCDF4 import Dataset, chartostring, date2num, num2date, default_fillvals
import pyproj

from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))

Load and select data from SEAScope

locate the WOC BFN QG current data in the SEAScope viewer

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

    • collection WOC BFN QG current, variable current

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

    • The 17th of December 2015 should be good

  3. Locate the data on the globe

    • It should be in the Agulhas current region

Steps for data extraction from SEAScope

  1. Zoom-in to properly see the current (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()

Print selected granules

#from SEAScope.lib import load_pyo
#data_path = ''
#extractions = load_pyo(data_path)
## Print a list of path of extracted granules 
print('\n'.join(extractions.keys()))

Play with the data

** Compute trajectory of fictive swell ray at initial position (x0, y0) using the dispertion relation in presence of surface current :
ωc=ω+k⃗.u⃗ \omega_c = \omega+ \vec{k}.\vec{u}
dxdt=dωcdk \frac{dx}{dt}=\frac{d\omega_c}{dk}
dkdt=−dωcdx \frac{dk}{dt}=-\frac{d\omega_c}{dx}
**

** Fill variable name for your selected velocity **

# Get longitudes, latitudes from gcps and time from first granule
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'])
name_var_u = 'ug'
name_var_v = 'vg'

extraction = extractions[granule_uri]   
u = extraction['data'][name_var_u]
v = extraction['data'][name_var_v]

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

velocity = (u**2+v**2)**(1./2)
print(velocity.shape)

Verify the extracted current by plotting it

fig = plt.figure(figsize = (16,10))
plt.imshow(numpy.flipud(velocity),interpolation='bicubic',cmap='jet', vmin=0 , vmax=2)
cbar=plt.colorbar()
X, Y = numpy.meshgrid(numpy.arange(0, u.shape[1]), numpy.arange(0, u.shape[0]))
Q = plt.quiver(X[::15, ::15], Y[::15, ::15],numpy.clip(numpy.flipud(u[::15,::15]),-2,2), numpy.clip(numpy.flipud(v[::15,::15]),-2,2))
qk = plt.quiverkey(Q, 0.5, 0.98, 2, r'$2 \frac{m}{s}$', labelpos='W',
                   fontproperties={'weight': 'bold'})
#plt.xlim(40, 150)
#plt.ylim(180, 120)
res = (lat2D[1,0] - lat2D[0,0])*100000
res = (lat2D[int(u.shape[0]/2+1), int(u.shape[1]/2)]
       - lat2D[int(u.shape[0]/2), int(u.shape[1]/2)]) * 100000
print(res)

Simulate wave current interactions

** Play with nstep (number of time steps) and nseed (number of rays) **

def ev(var, x, y):
    if numpy.isnan(x) or numpy.isnan(y):
        return numpy.nan
    nx, ny = numpy.shape(var)
    i = int(x)
    j = int(y)
    if (i > nx - 2) or (j > ny - 2):
        return numpy.nan
    dx = x - i
    dy = y - j
    uix = dx * var[i + 1, j] + (1 - dx) * var[i, j]
    uiy = dx * var[i + 1, j + 1] + (1 - dx) * var[i, j + 1]
    ui = dy * uiy + (1 - dy) * uix
    return ui


# Initialize wave ray position

granule_id = None
# Time step for advection (in sec)
s=1
dt = 3600/3/s
# Number of seeds
nseed = 120*2 #! Change this
# Wave number

# 17 dec 2015
wl = 250 
tetad = 35
tetad = (30-numpy.arange(nseed)/nseed*20)*0+35
c0 = 10 # Starting column
# Starting line
l0 = u.shape[0]-130
# Total time of advection 
nstep = int(130*s) #! Change this

k = 2*numpy.pi/wl
# Wave direction in radians
teta0 = tetad*numpy.pi/180
# Starting column

tmpu = extraction['data'][name_var_u]
tmpv = extraction['data'][name_var_v]
res = (lat2D[int(u.shape[0]/2+1), int(u.shape[1]/2)]
       -lat2D[int(u.shape[0]/2), int(u.shape[1]/2)]) * 100000
x = numpy.full((nstep+1, nseed), numpy.nan)
y = numpy.full((nstep+1, nseed), numpy.nan)
lon = numpy.full((nstep+1, nseed), numpy.nan)
lat = numpy.full((nstep+1, nseed), numpy.nan)
for c in range(nseed):
    x[0,c] = -l0 - 1.5*c/3
    y[0,c] = c0 + 6*c/3
    lon[0,c] = ev(lon2D, x[0,c], y[0,c])
    lat[0,c] = ev(lat2D, x[0,c], y[0,c])
    teta = teta0[c]
    kx = k*numpy.sin(teta)
    ky = k*numpy.cos(teta)
    for i in range(nstep):
        u = tmpu
        v = tmpv
        if x[i, c] is numpy.ma.masked or y[i, c] is numpy.ma.masked:
            continue
        omega = numpy.sqrt(9.81*(kx**2+ky**2)) + (ky*ev(u,x[i,c],y[i,c])
                            +kx*ev(v,x[i,c],y[i,c]))
        dx = dt*(ev(v,x[i,c],y[i,c])+numpy.sqrt(9.81/(numpy.sqrt(kx**2+ky**2)))
                 *numpy.sin(teta)/2)/res
        dy = dt*(ev(u,x[i,c],y[i,c])+numpy.sqrt(9.81/(numpy.sqrt(kx**2+ky**2)))
                 *numpy.cos(teta))/2/res
        teta = numpy.arctan2(kx,ky)
        if dx is numpy.ma.masked or dy is numpy.ma.masked:
            continue
        if (x[i, c]+dx<(1-u.shape[0])) or (y[i, c]+dy>(u.shape[1]-1)):
            continue
        if (x[i, c]+dx>0) or (y[i, c]+dy<0):
            continue
        domegax = omega - (numpy.sqrt(9.81*(kx**2+ky**2))
                           + (ky*ev(u,x[i,c]+dx,y[i,c])
                           + kx*ev(v,x[i,c]+dx,y[i,c])))
        domegay = omega - (numpy.sqrt(9.81*(kx**2+ky**2))
                           + (ky*ev(u,x[i,c],y[i,c]+dy)
                           + kx*ev(v,x[i,c],y[i,c]+dy)))

        x[i+1,c] = x[i,c] + dx
        y[i+1,c] = y[i,c] + dy
        lon[i+1,c] = ev(lon2D, x[i+1,c], y[i+1,c])
        lat[i+1,c] = ev(lat2D, x[i+1,c], y[i+1,c])
        dkx = dt*domegax/dx/res
        dky = dt*domegay/dy/res
        kx = kx + dkx
        ky = ky + dky

plot the simulated wave rays

fig = plt.figure(figsize = (18,10))
plt.imshow(numpy.flipud(velocity),interpolation='bicubic',cmap='jet', vmin=0 , vmax=2)
cbar=plt.colorbar()
Q = plt.quiver(X[::15, ::15], Y[::15, ::15],
               numpy.clip(numpy.flipud(u[::15,::15]),-2,2),
               numpy.clip(numpy.flipud(v[::15,::15]),-2,2))
qk = plt.quiverkey(Q, 0.5, 0.98, 2, r'$2 \frac{m}{s}$', labelpos='W',
                   fontproperties={'weight': 'bold'})
for c in range(nseed):
    plt.plot(y[:,c],-x[:,c],color='w',linewidth=1.5)

Export wave rays to SEAScope

Define Collection:
** A granule has to belong to a collection **

from SEAScope.lib.utils import create_collection, init_ids
# Make sure to keep the collection identifier, we will need it later
# collection_id, collection = create_collection('Simulated wave rays 11sec 25°- 10 aout 2019 - mercator GLOBAL OCEAN 1/12°')
# collection_id, collection = create_collection('Simulated wave rays 15sec 37°- 20 fev 2019 - BFN')
collection_id, collection = create_collection('Simulated wave rays - BFNQG ')
# collection_id, collection = create_collection('Simulated wave rays 15sec 37°- 20 fev 2019 - DUACS')
# collection_id, collection = create_collection('Simulated wave rays 15sec 35°- 20 fev 2019 - DUACS')
print(collection_id)

# IP address and port used to reach the SEAScope standalone application
host = '127.0.0.1'
port = 11155

import SEAScope.upload
with SEAScope.upload.connect(host, port) as link:
# Send create collection job
    SEAScope.upload.collection(link, collection)

Upload granule to SEAScope:
Define shell of granule, create granule
The granule will contain each matrix of modified data as what we call a “field”
 ** Enter a name for each field (field_name_time) **
Upload collection and granule to SEAScope

from SEAScope.lib.utils import create_granule
# Shell of Granule
# Each position of the trajectory will be used as GCP
with SEAScope.upload.connect(host, port) as link:   
    for c in range(nseed):
        gcps = []
        count = 0
        for i in range(0, nstep+1):
            gcp = {'lon': lon[i,c], 'lat': lat[i,c], 'i': i, 'j': 0}
            if lat[i,c] != 0:
                gcps.append(gcp)
                count=count+1

        # Create granule
        stop_dt = start + datetime.timedelta(days=7)
        granule_id, granule = create_granule(collection_id, gcps, start, stop_dt)

        # Set variables
        # Enter a name for each variable field
        from SEAScope.lib.utils import set_field
        field_name_time ='time rays'
        set_field(granule, field_name_time, numpy.arange(count)*dt/3600)
        
        # Send create granule job
        SEAScope.upload.granule(link, granule)

**Telling SEAScope what to display:
The next step is variables definition. **
Enter a name for each variable you want to display (variable_name_time, ...)

from SEAScope.lib.utils import create_variable

# Create a variable for each parameter
variable_name_time = 'time_rays'
var1 = create_variable(collection, variable_name_time, [field_name_time], dims=1)

# Upload to SEAScope
with SEAScope.upload.connect(host, port) as link:
    # Send create variable jobs
    SEAScope.upload.variable(link, var1)

** Telling SEAScope how to display variables:
The rendering options are defined in the following cell.
You can play with min, max values, the colormap ... and check how the rendering changes on SEAScope **

# Rendering configuration for the time
rcfg1 = var1['rendering']
rcfg1['min'] = 1
rcfg1['max'] = 40
rcfg1['colormap'] = ''
rcfg1['color'] = [255,255,255]
rcfg1['opacity'] = 0.25
rcfg1['zindex'] = .923


# Upload to SEAScope
with SEAScope.upload.connect(host, port) as link:
    # Send create rendering config jobs
    SEAScope.upload.rendering_config(link, rcfg1)
OVL Manual
Sentinel-2 swell
OVL Manual
SEAShot