Example

In this example, we use the CoherentStructures.jl package to work with some data derived from satelite altimetry. As is usual for Julia, this will take much longer than normal to run the first time as things are being compiled.

For the sake of simplicity we will restrict ourselves to the case of only a single worker, but parallelising is straightforward.

Getting data

Download velocity fields from Copernicus CMES. After having registered an account click on "Download", then "Download options" and then use the "FTP Access" to download the files you want.

In our case, we have stored the files in /foo/ww_ocean_data, an example filename is dt_global_allsat_phy_l4_20070414_20190101.nc.

Loading the data

Load the packages used:

using OceanTools, CoherentStructures, Dates, StaticArrays

The space-time window we are interested in:

start_date = DateTime("2006-11-25T00:00:00")
end_date = DateTime("2006-12-30T00:00:00")
LL_space = [-50.0, -50.0]
UR_space = [50.0, 50.0]

Read in velocity files (from "ugos" and "vgos" attributes in the NetCDF files)

schema = r"^dt_global_allsat_phy_l4_([0-9][0-9][0-9][0-9])([0-9][0-9])([0-9][0-9])_.*.nc$"
ww_ocean_data = "/foo/ww_ocean_data/"

p, ust1, (Lon, Lat, times)  = read_ocean_velocities(ww_ocean_data,start_date, end_date; schema=schema, LL_space=LL_space, UR_space=UR_space)

Plot a heatmap of the longitudinal component of a tricubic interpolation of the velocity field.

using Plots
xs = range(-20, stop=15, length=400)
ys = range(-10, stop=5, length=400)
t = times[10]
Plots.heatmap(xs, ys, (x,y) -> uv_tricubic((@SVector [x,y]), p, t)[1], title="Longitudinal velocity [deg/day]", color=:viridis, aspect_ratio=1.0)

Notice how areas where the velocity field is missing (i.e. on land), the velocity is zero. Because knowing land areas is sometimes useful, the ust1 variable contains a single time slice of the velocity where missing values are nan.

The velocity fields used are derived from sea surface height measurements. Load the sea surface height mesurements directly:

p2, _  = read_ocean_scalars(ww_ocean_data, start_date, end_date; schema=schema, LL_space=LL_space, UR_space=UR_space, scalar_field_name="sla")
Plots.heatmap(xs, ys, (x,y) -> scalar_tricubic((@SVector [x,y]), p2, t), title="Sea surface height anomaly [m]", color=:viridis, aspect_ratio=1.0)

It is straightforward to plot trajectories of single particles, we make an animation (download this script into "animations.jl" for the animatemp4 function):

include("animations.jl")
x0 = [0.0,0.0]
trajectory = flow(uv_tricubic,x0,times; p=p)

frames = []
for (i,t) in enumerate(times)
    fig = Plots.heatmap(xs,ys,(x,y) -> scalar_tricubic((@SVector [x,y]),p2,t ),title="Sea surface height anomaly [m]",color=:viridis,aspect_ratio=1.0,clim=(-0.1,0.1))
    Plots.scatter!(fig,(trajectory[i][1],trajectory[i][2]),label="simulated drifter position")
    push!(frames,fig)
end
animatemp4(frames)

The flow function used above came from the CoherentStructures.jl package (which in turn calls DifferentialEquations.jl). CoherentStructures.jl is also capable of fancier analyis:

vortices, singularities, bg = materialbarriers(
       uv_tricubic, range(-5,7.5,length=300), range(-40,stop=-28,length=300), range(times[2],stop=times[2]+30,length=30),
       LCSParameters(boxradius=4, indexradius=0.25, pmax=1.4,
                     merge_heuristics=[combine_20_aggressive]),
       p=p, on_torus=false);

fig = plot_vortices(vortices, singularities, (-5, -40), (7.5, -28);
    bg=bg, title="DBS field and transport barriers", showlabel=false)

Computations on larger domains are of course possible. Here is a paper that computes similar structures as above using OceanTools.jl and CoherentStructures.jl on a global scale.