Ray tracing in the Agulhas current ================================== .. code:: ipython3 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\ :math:`_{0}`, y\ :math:`_{0}`: are the locations where the rays are shot - plot_current: a flag to vizualize the current field. .. code:: ipython3 plot_current = True # flag to plot current .. code:: ipython3 T0 = 10 # Period [s] theta0_deg = 60 # Direction [deg] .. code:: ipython3 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 .. code:: ipython3 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) ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' .. code:: ipython3 ########## # -- 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 '''' .. code:: ipython3 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]') .. image:: ./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 (:math:`x_{0}`, :math:`y_{0}`) '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' .. code:: ipython3 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) ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' .. code:: ipython3 ucur = sub_ds_current.eastward_geostrophic_current_velocity.values.squeeze() vcur = sub_ds_current.northward_geostrophic_current_velocity.values.squeeze() .. code:: ipython3 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 ''''''''''''''''''''''' .. code:: ipython3 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 .. code:: ipython3 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') .. parsed-literal:: 500 rays computed in 0.10207271575927734 sec Remove 10 outliers rays (certainly due to discontinuity in current field) ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' .. code:: ipython3 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 ''''''''''''''''''''''''''''''''' .. code:: ipython3 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') .. image:: ./agulhas_example_23_0.png