ROMS - Forward tracking

PyLag includes support for reading in Regional Ocean Modelling System (ROMS) outputs and using them to track the movement of particles in the simulated flow field. Here we demonstrate this functionality using outputs from a multi-year numerical simulation of flow and water properties over the Texas-Louisiana contintal shelf. ROMS outputs for the region have been made available via a public thredds server that has been set up by the Physical Oceanography Numerical Group (PONG) from Texas A&M University. For convenience, a small subset of the full dataset has been pulled and uploaded to google drive with the other input datasets used when generating PyLag’s documentation. If you would like to run the code in this notebook interactively, download the data into a directory of your chooising. By default, the notebook will look for these files in the directory ${HOME}/data/pylag_doc. To change this, simply update the data_dir path below.

[1]:

import warnings
import matplotlib
import os
import pathlib

warnings.filterwarnings('ignore')

# Ensure inline plotting
%matplotlib inline

# Root directory for PyLag example input files
home_dir = os.environ['HOME']
data_dir = f'{home_dir}/data/pylag_doc'

# Keep a copy of the cwd
cwd = os.getcwd()

# Create run directory
simulation_dir = f'{cwd}/simulations/roms_texas'
pathlib.Path(simulation_dir).mkdir(parents=True, exist_ok=True)

Background

The Regional Ocean Modelling System (ROMS) is a free-surface, terrain-following, primitive equation ocean model. In the horizontal, the model employs a structured, curvilinear grid. The grid is a staggered, Arakawa C-grid. In the vertical, the model employs a generalized vertical, terrain-following coordinate system.

On the Arakawa C-grid, u-, v- and w- velocity components are evaluated on cell faces while all active and passive tracers are evaluated at cell centres. While ignoring those variables defined at cell corners, this effectively yields three separate horizontal grids for the U-, V- and rho-points; and two separate vertical grids for rho- and w-points. To handle grids of this type, PyLag contructs three separate horizontal unstructured grids and two separate vertical grids for the purpose of locating particles and interpolating field data at a given particle’s location.

The adopted approach works best for ROMS configurations in which the model grid is rectilinear, with u and v velocity components defined along lines of constant longitude and latitude. In this case, the model essentially reflects the approach taken in the LTRANS particle tracking model (North et al, 2011). Rudimentary support for setups involving curvilinear grids with u and v velocity components defined on the native model grid is also included. In such cases, PyLag first moves the u and v velocity components to rho points, then rotates the velocity vector onto a lat-lon geographic grid. It is anticipated this approach is likely sufficient for the majority of applications, especially those focussed on simulating particle dynamics and marine connectivity. However, there are two negative consequences of the adopted approach. First, the remapped velocity field is no longer identical to that generated by ROMS, which may negatively impact the representation of the flow field, especially around land boundaries; and b) model run times are increased, since rotating the velocity field, which is done on the fly, is a CPU intensive operation that invokes multiple calls to transcendental functions. In the future, to negate these problems, an approach similar to that used in TRACMASS (Döös et al, 2017) or Parcels (Delandmeter et al, 2019) could be implemented, although this would also necessitate extending the model to include more general support for tracking particles on quadrilateral grids.

The example uses outputs from a high resolution hindcast simulation of the Texas-Louisianna continental shelf region, based on a ROMS setup that uses a curvilinear grid. The hindcast has been run from 1993 through to the present day. More information on the model system along with links to scientific papers can be found here. In the example, forward in time particle tracks are created for a short (two day) simulation.

A map of the model domain with the bathymetry shaded is shown below.

[2]:

from netCDF4 import Dataset
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from pylag.processing.plot import create_figure, colourmap
from pylag.processing.plot import ArakawaCPlotter

# The name of the grid metrics file
grid_metrics_file_name = f'{data_dir}/txla_hindcast_grid_metrics.nc'

# Read the bathymetry
ds = Dataset(grid_metrics_file_name, 'r')
bathy = -ds.variables['h'][:]
ds.close()

# Create figure
font_size = 15
cmap = colourmap('h_r')
fig, ax = create_figure(figure_size=(26., 26.), projection=ccrs.PlateCarree(),
                        font_size=font_size, bg_color='white')

# Configure plotter
plotter = ArakawaCPlotter(grid_metrics_file_name,
                          geographic_coords=True,
                          font_size=font_size)

# Plot bathymetry
extents = None
ax, plot = plotter.plot_field(ax, 'grid_rho', bathy, extents=extents, add_colour_bar=True,
                              cb_label='Depth (m)', vmin=None, vmax=0., cmap=cmap)

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '10m',
                                            edgecolor='face',
                                            facecolor=cfeature.COLORS['land_alt1']));
../_images/examples_roms_txla_3_0.png

When making the above plot, grid information was read in from a grid metrics file, which has been pre-generated for the purpose of running the example.

Setting particle initial positions

We will create an initial position file in which to record particle initial positions. This will be read in by PyLag and used to create the initial particle seed. Particles will be released from a single release zone, which is located south of Marsh Island. The release zone will have a radius of 10000 m, and will contain 100 particles.

[3]:

from pylag.processing.coordinate import utm_from_lonlat, lonlat_from_utm
from pylag.processing.release_zone import create_release_zone
from pylag.processing.input import create_initial_positions_file_single_group

# The group ID of this particle set
group_id = 1

# Lat and lon coordiantes for the centre of the release zone
lat = 29.2
lon = -91.8

# Release zone radius (m)
radius = 10000.0

# Target number of particles to be released. Only a target,
# since we are evenly distributing particles in the release
# zone, which has no unique solution.
n_particles_target = 100

# Release depths
depth_below_surface = 0.0

# Create the release zone
surface_release_zone = create_release_zone(group_id = group_id,
                                           radius = radius,
                                           centre = [lon, lat],
                                           coordinate_system = 'geographic',
                                           n_particles = n_particles_target,
                                           depth = depth_below_surface,
                                           random = False)

# Get the actual number of particles
n_particles = surface_release_zone.get_number_of_particles()

# Get coordinates
lons, lats, depths = surface_release_zone.get_coordinates()

# Create input sub-directory
input_dir = f'{simulation_dir}/input'
pathlib.Path(input_dir).mkdir(parents=True, exist_ok=True)

# Output filename
file_name = f'{input_dir}/initial_positions.dat'

# Write data to file
create_initial_positions_file_single_group(file_name,
                                           n_particles,
                                           group_id,
                                           lons,
                                           lats,
                                           depths)

To see the initial positions of particles, we can simply plot them. Note the ordered scattering of particles within the circular release zone, which is clearly visibile after zooming in on the release zone.

[4]:

import numpy as np

# Create figure
font_size = 15
cmap = colourmap('h_r')
fig, ax = create_figure(figure_size=(26., 26.), projection=ccrs.PlateCarree(),
                        font_size=font_size, bg_color='gray')

# Configure plotter
plotter = ArakawaCPlotter(grid_metrics_file_name,
                          geographic_coords=True,
                          font_size=font_size)

# Plot bathymetry
extents = np.array([-92.5, -91.1, 28.6, 29.8], dtype=float)
ax, plot = plotter.plot_field(ax, 'grid_rho', bathy, extents=extents, add_colour_bar=True,
                              cb_label='Depth (m)', vmin=-60., vmax=0., cmap=cmap)

# Plot particle initial positions
plotter.scatter(ax, 'grid_rho', lons, lats, configure=True, extents=extents, s=16,
                color='#e50000', edgecolors='none');

# Add land
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '10m',
                                            edgecolor='face',
                                            facecolor=cfeature.COLORS['land_alt1']));
../_images/examples_roms_txla_8_0.png

Creating a ROMS grid metrics file

When running with a new input dataset, it is necessary to generate a grid metrics file. The grid metrics file holds information about the grid on which data are defined. PyLag reads this file in during model start up. The data is pre-generated, so that the grid information can be reused in later simulations using the same input dataset. This avoids incurring repeated costs associated with generating the data. As stated above, a grid metrics file for this example has been pre-generated; indeed, it was used to create the above plot of the domain’s bathymetry. However, for the sake of completeness, we regenerate the grid metrics file so it is clear how this is done in an end-to-end example. To generate a new ROMS grid metrics file, we need an ROMS output file for the domain in question. At a minimum, this should conain the list of dimensions and variables given below when running with a curvilinear grid. If the grid is rectilinear, the angle argument can be omitted. Further assistance on the creation of ROMS grid metrics file can be found in the help pages for the module pylag.grid_metrics.

[5]:

from pylag.grid_metrics import create_roms_grid_metrics_file

# A ROMS output file
roms_file_name = f'{data_dir}/txla_hindcast_agg_0001.nc'

# The name of the output file
grid_metrics_file_name = f'{input_dir}/grid_metrics.nc'.format(input_dir)

# Optional variable names
angles_var_name='angle'
vtransform=2

# Generate the file
create_roms_grid_metrics_file(roms_file_name,
                              angles_var_name=angles_var_name,
                              grid_metrics_file_name=grid_metrics_file_name,
                              vtransform=vtransform)

Reading variable `h` ... done
Sorting axes for variable `h`... done
Reading variable `angle` ... done
Sorting axes for variable `angle`... done

Processing grid_u variables ...
Reading variable `lon_u` ... done
Reading variable `lat_u` ... done
Sorting axes for variable `lon_u`... done
Sorting axes for variable `lat_u`... done
Calculating lons and lats at element centres ... done

Processing grid_v variables ...
Reading variable `lon_v` ... done
Reading variable `lat_v` ... done
Sorting axes for variable `lon_v`... done
Sorting axes for variable `lat_v`... done
Calculating lons and lats at element centres ... done

Processing grid_rho variables ...
Reading variable `lon_rho` ... done
Reading variable `lat_rho` ... done
Sorting axes for variable `lon_rho`... done
Sorting axes for variable `lat_rho`... done
Calculating lons and lats at element centres ... done

Processing grid_psi variables ...
Reading variable `lon_psi` ... done
Reading variable `lat_psi` ... done
Sorting axes for variable `lon_psi`... done
Sorting axes for variable `lat_psi`... done
Calculating lons and lats at element centres ... done

Calculating land sea mask ...
Reading variable `mask_rho` ... done
Sorting axes for variable `mask_rho`... done

Reading vertical grid vars
Reading variable `s_rho` ... done
Reading variable `Cs_r` ... done
Reading variable `s_w` ... done
Reading variable `Cs_w` ... done
Reading variable `hc` ... done

Creating grid metrics file /users/modellers/jcl/code/git/PyLag/PyLag/doc/source/examples/simulations/roms_texas/input/grid_metrics.nc

Creating the run configuration file

As in the other examples, a template run configuration file has been provided. We will use configparser to look at some settings specific to the example. First, we read in the file, then we print out some of the key options for the run:

[6]:

import configparser

config_file_name = './configs/roms_txla_template.cfg'

cf = configparser.ConfigParser()
cf.read(config_file_name)

# Start time
print(f"Start time: {cf.get('SIMULATION', 'start_datetime')}")

# End time
print(f"End time: {cf.get('SIMULATION', 'end_datetime')}")

# Specify that this is a forward tracking experiment
print(f"Time direction: {cf.get('SIMULATION', 'time_direction')}")

# We will do a single run, rather than an ensemble run
print(f"Number of particle releases: {cf.get('SIMULATION', 'number_of_particle_releases')}")

# Use depth restoring, and restore particle depths to the ocean surface
print(f"Use depth restoring: {cf.get('SIMULATION', 'depth_restoring')}")
print(f"Restore particles to a depth of: {cf.get('SIMULATION', 'fixed_depth')} m")

# Specify that we are working with Arakawa A-grid in geographic coordinates
print(f"Model name: {cf.get('OCEAN_DATA', 'name')}")
print(f"Coordinate system: {cf.get('SIMULATION', 'coordinate_system')}")

# Set the location of the grid metrics and input files
print(f"Data directory: {cf.get('OCEAN_DATA', 'data_dir')}")
print(f"Path to grid metrics file: {cf.get('OCEAN_DATA', 'grid_metrics_file')}")
print(f"File name stem of input files: {cf.get('OCEAN_DATA', 'data_file_stem')}")

Start time: 2018-06-29 00:00:00
End time: 2018-06-30 23:00:00
Time direction: forward
Number of particle releases: 1
Use depth restoring: True
Restore particles to a depth of: 0.0 m
Model name: ROMS
Coordinate system: geographic
Data directory:
Path to grid metrics file:
File name stem of input files: txla_hindcast_agg

From the above options, it can be seen we hvae not yet specified a directory within which the model can find ROMS output files, or the location of the grid metrics file. We will set these using data_dir, and save the new config file in the simulation directory. All other paths are relative, and can be left as they are.

[7]:

cf.set('OCEAN_DATA', 'data_dir', data_dir)
cf.set('OCEAN_DATA', 'grid_metrics_file', grid_metrics_file_name)

# Save a copy in the simulation directory
with open(f"{simulation_dir}/pylag.cfg", 'w') as config:
    cf.write(config)

Running the model

With the run configuration file saved, we can now run the example. While PyLag can be used interactivatly, it is most commonly launched from the command line like so:

[8]:

# Change to the run directory
os.chdir(f'{simulation_dir}')

# Run the model
!{"python -m pylag.main -c pylag.cfg"}

# Return to the cwd
os.chdir(cwd)


Starting ensemble member 1 ...
Progress:
100% |###########################################|

Visualising the result

With the model having run, the final step is to visulise the result. First, we produce a animation of the velocity field over the course of the simulation and near to where the particles were released. We do this by first moving the u and v velocity components to rho-points, then rotating the velocity field onto a geographic lat-lon grid. For this, we leverage the seapy library. This is the same procedure that PyLag (currently) uses internally when working with ROMS outputs on a curvilinear grid.

[10]:

from netCDF4 import num2date
from matplotlib import pyplot as plt
from matplotlib import animation
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from IPython.display import HTML
import seapy

from pylag.processing.ncview import Viewer

# Open the dataset
viewer = Viewer(f'{data_dir}/txla_hindcast_agg_0001.nc'.format(data_dir))

# Get dates
time = viewer('ocean_time')
date = num2date(time[:], units=time.units)

# Plot surface velocity field on model grid at rho points
lon_rho = viewer('lon_rho')[:]
lat_rho = viewer('lat_rho')[:]
u = viewer('u')[:, -1, :, :]
v = viewer('v')[:, -1, :, :]
angle = viewer('angle')[:]

# Move to rho points
u = seapy.model.u2rho(u)
v = seapy.model.v2rho(v)

# Rotate
u_rot, v_rot = seapy.rotate(u, v, angle)

# Plot configuration
extents = np.array([-92.5, -91.1, 28.6, 29.8], dtype=float)
step = 5

# Set up the plot
fig, ax = create_figure(figure_size=(26, 26), projection=ccrs.PlateCarree())

gl = ax.gridlines(linewidth=1, draw_labels=True, linestyle='--', color='k')

gl.xlabel_style = {'fontsize': 10}
gl.ylabel_style = {'fontsize': 10}

gl.xlabels_top=False
gl.ylabels_right=False
gl.xlabels_bottom=True
gl.ylabels_left=True

gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER

# Make the animation
quiver = ax.quiver(lon_rho[::step, ::step], lat_rho[::step, ::step],
                   u_rot[0, ::step, ::step], v_rot[0, ::step, ::step],
                   units='inches', scale_units='inches', scale=0.5)
plt.quiverkey(quiver, 0.9, 0.9, 0.5, '{} '.format(0.5) + r'$\mathrm{ms^{-1}}$',
              coordinates='axes')

ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m',
                                            edgecolor='face',
                                            facecolor=cfeature.COLORS['land_alt1']))

if extents is not None:
    ax.set_extent(extents, ccrs.PlateCarree())

# For animating the field ...

def init():
    quiver.set_UVC(u_rot[0, ::step, ::step], v_rot[0, ::step, ::step])
    return (quiver,)

def animate(i):
    ax.set_title('{}'.format(date[i]))
    quiver.set_UVC(u_rot[i, ::step, ::step], v_rot[i, ::step, ::step])
    return (quiver,)

# Prevent the basic figure from being plotted too
plt.close(fig)

anim = animation.FuncAnimation(fig, animate, init_func=init, frames=50, interval=50, blit=True)

HTML(anim.to_html5_video())

[10]:

Next we plot particle trajectories over the course of the simulation.

[10]:

from datetime import timedelta

# Dataset holding PyLag ouputs
file_name = '{}/output/pylag_1.nc'.format(simulation_dir)
pylag_viewer = Viewer(file_name, time_rounding=900)

# Create figure
font_size = 15
cmap = colourmap('h_r')
fig, ax = create_figure(figure_size=(26., 26.), projection=ccrs.PlateCarree(),
                        font_size=font_size, bg_color='gray')

# Configure plotter
plotter = ArakawaCPlotter(grid_metrics_file_name,
                          geographic_coords=True,
                          font_size=font_size)

# Plot bathymetry
extents = np.array([-92.5, -91.1, 28.6, 29.8], dtype=float)
ax, plot = plotter.plot_field(ax, 'grid_rho', bathy, extents=extents,
                              add_colour_bar=True, cb_label='Depth (m)',
                              vmin=-60., vmax=0., cmap=cmap)

# Add land
ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m',
                                            edgecolor='face',
                                            facecolor=cfeature.COLORS['land_alt1']))

# Time of flight
time_of_flight = timedelta(hours=47)

# Date index for plotting
pylag_date = pylag_viewer.date[0] + time_of_flight
time_index = pylag_viewer.date.tolist().index(pylag_date)

# Plot particle final positions
ax, scatter = plotter.scatter(ax, 'grid_rho', pylag_viewer('longitude')[time_index, :].squeeze(),
                             pylag_viewer('latitude')[time_index, :].squeeze(),
                             s=8, color='#e50000', edgecolors='none')

# Add path lines
ax, lines = plotter.plot_lines(ax, pylag_viewer('longitude')[:time_index, :],
                               pylag_viewer('latitude')[:time_index, :],
                               linewidth=0.15, alpha=1,
                               color='#0504aa')
../_images/examples_roms_txla_20_0.png

It is by no means a verification that the particles are moving as they should be, but at a basic level it can be seen the particle tracks are following the velocity vectors, as they should.

References

  1. North, E. W., E. E. Adams, S. Schlag, C. R. Sherwood, R. He, S. Socolofsky. 2011. Simulating oil droplet dispersal from the Deepwater Horizon spill with a Lagrangian approach. AGU Book Series: Monitoring and Modeling the Deepwater Horizon Oil Spill: A Record Breaking Enterprise.

  2. Döös, K., Jönsson, B., and Kjellsson, J.: Evaluation of oceanic and atmospheric trajectory schemes in the TRACMASS trajectory model v6.0, Geosci. Model Dev., 10, 1733–1749, https://doi.org/10.5194/gmd-10-1733-2017, 2017

  3. Delandmeter, P. and van Sebille, E.: The Parcels v2.0 Lagrangian framework: new field interpolation schemes, Geosci. Model Dev., 12, 3571–3584, https://doi.org/10.5194/gmd-12-3571-2019, 2019.