Note
Go to the end to download the full example code.
Leeway backtracking
import os
from datetime import timedelta
import cmocean
import pyproj
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import opendrift
from opendrift.readers import reader_netCDF_CF_generic
from opendrift.models.leeway import Leeway
We try to find the likelihood of the origin of a found object by two different methods: 1. backwards simulation from position where object is found (‘Observation’) 2. forwards simulation from a uniform grid of possible initial locations, selecting the origins of particles actually hitting the observed target
We use 24 hours from the NorKyst ocean model (800m pixel size) and Arome atmospheric model (2.5km pixel size)
o = Leeway(loglevel=50)
reader_arome = reader_netCDF_CF_generic.Reader(o.test_data_folder() +
'16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc')
reader_norkyst = reader_netCDF_CF_generic.Reader(o.test_data_folder() +
'16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc')
o.add_reader([reader_norkyst, reader_arome])
duration = timedelta(hours=24)
start_time = reader_norkyst.start_time
end_time = start_time + duration
object_type = 26 # 26 = Life-raft, no ballast
outfile = 'leeway.nc'
ilon = 4.3 # Incident position
ilat = 60.6
text = [{'s': 'Observation', 'x': ilon, 'y': ilat, 'fontsize': 20, 'color': 'g', 'zorder': 1000}]
# Define domain of possible origin
lons = np.arange(3.4, 5, .1/20)
lats = np.arange(59.7, 60.8, .05/20)
corners = [lons[0], lons[-1], lats[0], lats[-1]]
lons, lats = np.meshgrid(lons, lats)
/opt/conda/envs/opendrift/lib/python3.11/site-packages/pyproj/crs/crs.py:1286: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
proj = self._crs.to_proj4(version=version)
/opt/conda/envs/opendrift/lib/python3.11/site-packages/pyproj/crs/crs.py:1286: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
proj = self._crs.to_proj4(version=version)
Simulating first backwards for 24 hours:
o.seed_elements(lon=ilon, lat=ilat, radius=5000, radius_type='uniform', number=30000,
time=end_time, object_type=object_type)
o.run(duration=duration, time_step=-900, time_step_output=3600, outfile=outfile)
od = opendrift.open_xarray(outfile)
density_backwards = od.get_histogram(pixelsize_m=5000).isel(time=-1).isel(origin_marker=0)
density_backwards = density_backwards.where(density_backwards>0)
density_backwards = density_backwards/density_backwards.sum()*100
vmax = density_backwards.max()
o.plot(background=density_backwards, clabel='Probability of origin [%]', text=text, corners=corners, fast=True, markersize=.5, lalpha=.02, vmin=0, vmax=vmax)
os.remove(outfile)
/root/project/opendrift/readers/interpolation/interpolators.py:17: RuntimeWarning: overflow encountered in cast
data[mask] = np.finfo(np.float64).min
/opt/conda/envs/opendrift/lib/python3.11/site-packages/numpy/ma/core.py:467: RuntimeWarning: invalid value encountered in cast
fill_value = np.array(fill_value, copy=False, dtype=ndtype)
14:23:34 DEBUG opendrift.config:161: Adding 16 config items from __init__
14:23:34 DEBUG opendrift.config:171: Overwriting config item readers:max_number_of_fails
14:23:34 DEBUG opendrift.config:161: Adding 5 config items from __init__
14:23:34 INFO opendrift.models.basemodel:533: OpenDriftSimulation initialised (version 1.11.2 / v1.11.2-43-gce9f6bb)
14:23:34 DEBUG opendrift.config:161: Adding 6 config items from leeway
14:23:34 DEBUG opendrift.config:171: Overwriting config item seed:jibe_probability
14:23:34 DEBUG opendrift.export.io_netcdf:226: Importing with Xarray from leeway.nc
14:23:35 INFO opendrift:118: Returning <class 'opendrift.models.leeway.Leeway'> object
14:23:35 INFO opendrift.models.basemodel:3854: calculating for origin_marker 0...
14:23:35 DEBUG opendrift.models.basemodel:2333: Setting up map: corners=[3.4, 4.9949999999999655, 59.7, 60.797499999999005], fast=True, lscale=None
14:23:35 WARNING opendrift.models.basemodel:2379: Plotting fast. This will make your plots less accurate.
Simulating then forwards, starting at a uniform grid 24 hours earlier (440 x 320 = 140800 elements at ~500m separation)
o = Leeway(loglevel=50)
o.add_reader([reader_norkyst, reader_arome])
o.seed_elements(lon=lons, lat=lats, radius=0,
time=start_time, object_type=object_type)
o.run(duration=duration, time_step=900, time_step_output=3600, outfile=outfile)
print(o)
14:25:08 DEBUG opendrift.config:161: Adding 10 config items from environment
14:25:08 DEBUG opendrift.config:161: Adding 5 config items from environment
/root/project/opendrift/readers/interpolation/interpolators.py:17: RuntimeWarning: overflow encountered in cast
data[mask] = np.finfo(np.float64).min
/opt/conda/envs/opendrift/lib/python3.11/site-packages/numpy/ma/core.py:467: RuntimeWarning: invalid value encountered in cast
fill_value = np.array(fill_value, copy=False, dtype=ndtype)
===========================
--------------------
Reader performance:
--------------------
/root/project/tests/test_data/16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc
0:00:22.0 total
0:00:00.1 preparing
0:00:03.8 reading
0:00:03.3 interpolation
0:00:00.0 interpolation_time
0:00:17.7 rotating vectors
0:00:00.1 masking
--------------------
/root/project/tests/test_data/16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc
0:00:21.0 total
0:00:00.1 preparing
0:00:03.4 reading
0:00:02.9 interpolation
0:00:00.0 interpolation_time
0:00:17.2 rotating vectors
0:00:00.1 masking
--------------------
global_landmask
0:00:01.2 total
0:00:00.0 preparing
0:00:01.0 reading
0:00:00.0 masking
--------------------
Performance:
1:27.1 total time
0.6 configuration
10.2 preparing main loop
0.0 moving elements to ocean
1:15.6 main loop
12.7 updating elements
0.5 cleaning up
--------------------
===========================
Model: Leeway (OpenDrift version 1.11.2)
97169 active LeewayObj particles (43631 deactivated, 0 scheduled)
-------------------
Environment variables:
-----
x_sea_water_velocity
y_sea_water_velocity
1) /root/project/tests/test_data/16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc
-----
x_wind
y_wind
1) /root/project/tests/test_data/16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc
-----
land_binary_mask
1) global_landmask
Discarded readers:
Time:
Start: 2015-11-16 00:00:00 UTC
Present: 2015-11-17 00:00:00 UTC
Calculation steps: 96 * 0:15:00 - total time: 1 day, 0:00:00
Output steps: 25 * 1:00:00
-------------------
"Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind']", "Missing variables: ['x_wind', 'y_wind']", "Missing variables: ['x_wind', 'y_wind']", "Missing variables: ['x_wind', 'y_wind']", "Missing variables: ['x_wind', 'y_wind']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_wind', 'y_wind', 'x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']", "Missing variables: ['x_sea_water_velocity', 'y_sea_water_velocity']"
===========================
Finding the elements actually hitting the target (within 5 km) after 24 hours:
lon, lat = o.get_lonlats()
lonend = lon[:, -1]
latend = lat[:, -1]
geod = pyproj.Geod(ellps='WGS84')
on = np.ones(lonend.shape)
dummy1, dummy2, dist2incident = geod.inv(lonend, latend, ilon*on, ilat*on)
hits = np.where(dist2incident<5000)[0]
hit_start_lons = lon[hits, 0]
hit_start_lats = lat[hits, 0]
o_hit = opendrift.open(outfile, elements=hits)
o.animation(compare=o_hit, legend=['Elements not hitting target', 'Elements hitting target'], fast=True, corners=corners, text=text)
14:26:48 DEBUG opendrift.config:161: Adding 16 config items from __init__
14:26:48 DEBUG opendrift.config:171: Overwriting config item readers:max_number_of_fails
14:26:48 DEBUG opendrift.config:161: Adding 5 config items from __init__
14:26:48 INFO opendrift.models.basemodel:533: OpenDriftSimulation initialised (version 1.11.2 / v1.11.2-43-gce9f6bb)
14:26:48 DEBUG opendrift.config:161: Adding 6 config items from leeway
14:26:48 DEBUG opendrift.config:171: Overwriting config item seed:jibe_probability
14:26:48 DEBUG opendrift.export.io_netcdf:281: Importing from leeway.nc
/opt/conda/envs/opendrift/lib/python3.11/site-packages/numpy/ma/core.py:467: RuntimeWarning: invalid value encountered in cast
fill_value = np.array(fill_value, copy=False, dtype=ndtype)
14:26:49 DEBUG opendrift.models.basemodel:1683: No elements to deactivate
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: environment:constant:x_wind -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: environment:fallback:x_wind -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: environment:constant:y_wind -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: environment:fallback:y_wind -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: environment:constant:x_sea_water_velocity -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: environment:fallback:x_sea_water_velocity -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: environment:constant:y_sea_water_velocity -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: environment:fallback:y_sea_water_velocity -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: environment:constant:land_binary_mask -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: environment:fallback:land_binary_mask -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: general:use_auto_landmask -> True
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:current_uncertainty -> 0
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:current_uncertainty_uniform -> 0
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:max_speed -> 5.0
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: readers:max_number_of_fails -> 1
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: general:coastline_action -> stranding
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: general:time_step_minutes -> 10.0
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: general:time_step_output_minutes -> 60.0
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: seed:ocean_only -> True
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: seed:number -> 1
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:max_age_seconds -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:advection_scheme -> euler
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:horizontal_diffusivity -> 0
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:profiles_depth -> 50
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:wind_uncertainty -> 0
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:relative_wind -> False
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:deactivate_north_of -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:deactivate_south_of -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:deactivate_east_of -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: drift:deactivate_west_of -> None
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: seed:origin_marker -> 0
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: seed:z -> 0
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: seed:jibe_probability -> 0.04
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: seed:capsized -> 0
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: seed:current_drift_factor -> 1
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: seed:object_type -> Person-in-water (PIW), unknown state (mean values)
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: capsizing -> False
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: capsizing:leeway_fraction -> 0.4
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: capsizing:wind_threshold -> 30
14:26:49 DEBUG opendrift.export.io_netcdf:380: Setting imported config: capsizing:wind_threshold_sigma -> 5
14:26:49 INFO opendrift:78: Returning <class 'opendrift.models.leeway.Leeway'> object
14:26:49 DEBUG opendrift.models.basemodel:2333: Setting up map: corners=[3.4, 4.9949999999999655, 59.7, 60.797499999999005], fast=True, lscale=None
14:26:49 WARNING opendrift.models.basemodel:2379: Plotting fast. This will make your plots less accurate.
/opt/conda/envs/opendrift/lib/python3.11/site-packages/cartopy/mpl/geoaxes.py:1683: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
result = super().scatter(*args, **kwargs)
14:26:51 DEBUG opendrift.models.basemodel:3010: Saving animation..
14:26:52 INFO opendrift.models.basemodel:4570: Saving animation to /root/project/docs/source/gallery/animations/example_leeway_backtrack_0.gif...
14:28:16 DEBUG opendrift.models.basemodel:4608: MPLBACKEND = agg
14:28:16 DEBUG opendrift.models.basemodel:4609: DISPLAY = None
14:28:16 DEBUG opendrift.models.basemodel:4610: Time to save animation: 0:01:24.900426
14:28:16 INFO opendrift.models.basemodel:3003: Time to make animation: 0:01:27.260269
o.plot(compare=o_hit, legend=['Elements not hitting target', 'Elements hitting target'], show_elements=False, fast=True, corners=corners, text=text)
14:28:16 DEBUG opendrift.models.basemodel:2333: Setting up map: corners=[3.4, 4.9949999999999655, 59.7, 60.797499999999005], fast=True, lscale=None
14:28:16 WARNING opendrift.models.basemodel:2379: Plotting fast. This will make your plots less accurate.
(<GeoAxes: title={'center': 'OpenDrift - Leeway (LIFE-RAFT-NB-1)\n2015-11-16 00:00 to 2015-11-17 00:00 UTC (25 steps)'}>, <Figure size 793.298x1100 with 1 Axes>)
Plot the initial density of elements that actually hit the target after 24 hours. To be compared with the density figure from backwards simulation (see top)
of = opendrift.open_xarray(outfile, elements=hits)
density_forwards = of.get_histogram(pixelsize_m=5000).isel(time=0).isel(origin_marker=0)
density_forwards = density_forwards.where(density_forwards>0)
o_hit.plot(background=density_forwards/density_forwards.sum()*100, clabel='Probability of origin [%]', text=text, corners=corners, fast=True, markersize=.5, lalpha=.02, vmin=0, vmax=vmax)
14:37:08 DEBUG opendrift.config:161: Adding 10 config items from environment
14:37:08 DEBUG opendrift.config:161: Adding 5 config items from environment
14:37:09 DEBUG opendrift.config:161: Adding 16 config items from __init__
14:37:09 DEBUG opendrift.config:171: Overwriting config item readers:max_number_of_fails
14:37:09 DEBUG opendrift.config:161: Adding 5 config items from __init__
14:37:09 INFO opendrift.models.basemodel:533: OpenDriftSimulation initialised (version 1.11.2 / v1.11.2-43-gce9f6bb)
14:37:09 DEBUG opendrift.config:161: Adding 6 config items from leeway
14:37:09 DEBUG opendrift.config:171: Overwriting config item seed:jibe_probability
14:37:09 DEBUG opendrift.export.io_netcdf:226: Importing with Xarray from leeway.nc
14:37:09 INFO opendrift:118: Returning <class 'opendrift.models.leeway.Leeway'> object
14:37:09 INFO opendrift.models.basemodel:3854: calculating for origin_marker 0...
14:37:09 DEBUG opendrift.models.basemodel:2333: Setting up map: corners=[3.4, 4.9949999999999655, 59.7, 60.797499999999005], fast=True, lscale=None
14:37:09 WARNING opendrift.models.basemodel:2379: Plotting fast. This will make your plots less accurate.
(<GeoAxes: title={'center': 'OpenDrift - Leeway (LIFE-RAFT-NB-1)\n2015-11-16 00:00 to 2015-11-17 00:00 UTC (25 steps)'}>, <Figure size 793.298x1100 with 2 Axes>)
Total running time of the script: (14 minutes 22.349 seconds)