Ray tracing in the Agulhas current

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import time

import mantaray

from src import *

Introduction

This notebook test the Ray Tracing Package in a realistic analysis. We focus on the wave-trapping process occuring in the Agulhas current along the South Africa Coast. We use a realistic current field and a constant bathymetry field. Both fields are already constructed but can be re-constructed by the use of the present notebook. Before running this example, it is recommended that users run idealized_fields.ipynb first

Initialization:

  • The initial conditions (wave wavelength/period/wavenumber, and direction) for the present example are inspired from the Halsne et al. 2023. (current data are from the globcurrent product: Rio, M. H., Mulet, S., & Picot, N. (2014). Beyond GOCE for the ocean circulation estimate: Synergetic use of altimetry, gravimetry, and in situ data provides new insight into geostrophic and Ekman currents. Geophysical Research Letters, 41(24), 8918-8925., https://data.marine.copernicus.eu/product/MULTIOBS_GLO_PHY_MYNRT_015_003/description)

  • num_rays: is the number of rays shot. All rays have the same initial parameters

  • x\(_{0}\), y\(_{0}\): are the locations where the rays are shot

  • plot_current: a flag to vizualize the current field.

plot_current = True # flag to plot current

T0 = 10 # Period [s]
theta0_deg = 60 # Direction [deg]

theta0 = theta0_deg * np.pi/180
g = 9.81 # Acceleration of Gravity

# Convert period to wavenumber magnitude
f = (1/T0)
k0 = (2*np.pi*f)**2/g

n_rays = 500 # number of rays to shot
depth0 = 4000 # constant depth for this case

Load data and convert lon/lat to meter axes. Define a sub domain that will be cropped from the global current dataset (Globcurrent product Rio et al. 2014)

##########
# -- Define the subdomain
##########
lon_min = 16
lon_max = 29.5

lat_min = -42
lat_max = -32

# Create a rectangle
lons_domain = [lon_min, lon_max, lon_max, lon_min, lon_min]
lats_domain = [lat_min, lat_min, lat_max, lat_max, lat_min]

##########
# -- Define the starting line for ray tracing computation
##########
x_line = [10_000, 400_000]
y_line = [200_000, 100_000]

##########
# -- Load current data and crop it
##########
ds_current = xr.open_dataset('../data/globcurrent.nc')
ds_current['current_speed'] = (ds_current.eastward_geostrophic_current_velocity**2 + ds_current.northward_geostrophic_current_velocity**2)**(1/2)

sub_ds_current = ds_current.sel(lon = slice(lon_min, lon_max), lat = slice(lat_min, lat_max))


# --- Convert lon/lat to meter/meter
Lon_meter = ((sub_ds_current.lon.values-np.mean(sub_ds_current.lon.values))*1_852*60*np.cos(np.mean(sub_ds_current.lat.values)*np.pi/180))
Lat_meter = ((sub_ds_current.lat.values-np.mean(sub_ds_current.lat.values))*1_852*60)

Lon_meter += Lon_meter.max()
Lat_meter += Lat_meter.max()

Plot

if plot_current:
    fig, axes = plt.subplots(ncols = 2, figsize = (10, 5))
    ax = axes[0]
    ds_current.current_speed.squeeze().plot(ax =ax, vmin = 0, vmax = 2, add_colorbar = False)
    ax.fill_between(lons_domain, lats_domain, facecolor = 'r', alpha = .6)
    ax.set_aspect('equal', 'box')

    ax = axes[1]
    p2 = sub_ds_current.current_speed.squeeze().plot(ax = ax, vmin = 0, vmax = 2, add_colorbar = False)
    cax = fig.add_axes([.1, .02, .8, 0.02])
    cbar = plt.colorbar(p2, cax = cax, orientation = 'horizontal')
    ax.set_aspect('equal', 'box')
    ax.set_title(' ')
    cbar.ax.set_xlabel('Current Speed [m/sec]')
_images/agulhas_example_11_0.png

Define the initial location of rays. Instead of shotting rays from the edges of the domain, we create a dummy line that will be used to defined starting points (\(x_{0}\), \(y_{0}\))

x0, y0 = start_line_ray_tracing(x_line, y_line,  Lon_meter, Lat_meter, n_rays) # starting locations

Create forcings for ray tracing (Not mandatory because already created)

ucur = sub_ds_current.eastward_geostrophic_current_velocity.values.squeeze()
vcur = sub_ds_current.northward_geostrophic_current_velocity.values.squeeze()

output_file_cur = "../data/current_agulhas.nc"
output_file_depth =  "../data/bathy_agulhas.nc"
create_current_forcing(Lon_meter, Lat_meter, ucur, vcur, output_file_cur, output_file_depth, depth =  depth0)

Perform the ray tracing

cg = group_velocity(k0, f, depth0)

# Calculate wavenumber components
kx0 = k0*np.cos(theta0)
ky0 = k0*np.sin(theta0)

# Number of rays
# Initialize wavenumber for all rays
Kx0 = kx0*np.ones(n_rays)
Ky0 = ky0*np.ones(n_rays)

# Current and bathymetry file path
current = output_file_cur
bathymetry = output_file_depth

# Read x and y from file to get domain size
ds = xr.open_dataset(current)

# Estimates CFL
# Computes grid smallest spacing
dd = np.min([np.diff(ds.x).mean(), np.diff(ds.y).mean()])
# Computes group velocity
cg = group_velocity(k0, f, depth0)
# Computes CFL
cfl = dd/cg

duration = round(ds.x.max().values/cg)
step_size = cfl

t = time.time() # - How long this ray tracing process takes? (1/2)
bando = mantaray.ray_tracing(x0, y0, Kx0, Ky0, duration, step_size, bathymetry, current)
elapsed = time.time() - t # - How long this ray tracing process takes? (2/2)
print(f'{n_rays} rays computed in {elapsed} sec')

500 rays computed in 0.10207271575927734 sec

Remove 10 outliers rays (certainly due to discontinuity in current field)

id_rays1 = np.array(np.arange(0, 370, 1), dtype = int)
id_rays2 = np.array(np.arange(380, 500, 1), dtype = int)
id_ray_final = np.concatenate([id_rays1, id_rays2])

Final Plot In the Agulhas current

fig, ax = plt.subplots()

ax.plot(x_line, y_line, color = 'w')

cs = ax.pcolormesh(Lon_meter, Lat_meter, (ds.u**2 + ds.v**2)**(1/2), vmin = 0, vmax = 2)
for i in id_ray_final:
    ray = bando.isel(ray=i)
    ax.plot(ray.x, ray.y, 'r', lw=.78)
cbar = plt.colorbar(cs, orientation = 'vertical')
cbar.ax.set_ylabel('Current Speed [m/sec]')
ax.set_xlabel('Distance [km]')
ax.set_ylabel('Distance [km]')

ax.set_yticks([0, 250_000, 500_000, 750_000, 1_000_000])
ax.set_xticks([0, 250_000, 500_000, 750_000, 1_000_000])

ax.set_yticklabels([0, 250, 500, 750, 1_000])
ax.set_xticklabels([0, 250, 500, 750, 1_000])

ax.set_aspect('equal', 'box')
_images/agulhas_example_23_0.png