{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Leeway backtracking\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import os\nfrom datetime import timedelta\nimport cmocean\nimport pyproj\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport xarray as xr\nimport opendrift\nfrom opendrift.readers import reader_netCDF_CF_generic\nfrom opendrift.models.leeway import Leeway" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We try to find the likelihood of the origin of a found object by two different methods:\n1. backwards simulation from position where object is found ('Observation')\n2. forwards simulation from a uniform grid of possible initial locations, selecting the origins of particles actually hitting the observed target\n\nWe use 24 hours from the NorKyst ocean model (800m pixel size) and Arome atmospheric model (2.5km pixel size)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "o = Leeway(loglevel=50)\nreader_arome = reader_netCDF_CF_generic.Reader(o.test_data_folder() +\n '16Nov2015_NorKyst_z_surface/arome_subset_16Nov2015.nc')\nreader_norkyst = reader_netCDF_CF_generic.Reader(o.test_data_folder() +\n '16Nov2015_NorKyst_z_surface/norkyst800_subset_16Nov2015.nc')\no.add_reader([reader_norkyst, reader_arome])\n\nduration = timedelta(hours=24)\nstart_time = reader_norkyst.start_time\nend_time = start_time + duration\n\nobject_type = 26 # 26 = Life-raft, no ballast\noutfile = 'leeway.nc'\nilon = 4.3 # Incident position\nilat = 60.6\ntext = [{'s': 'Observation', 'x': ilon, 'y': ilat, 'fontsize': 20, 'color': 'g', 'zorder': 1000}]\n# Define domain of possible origin\nlons = np.arange(3.4, 5, .1/20)\nlats = np.arange(59.7, 60.8, .05/20)\n#lons = lons[0::2] # using every second, due to memory limitation on CircleCI\n#lats = lats[0::2]\ncorners = [lons[0], lons[-1], lats[0], lats[-1]]\nlons, lats = np.meshgrid(lons, lats)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simulating first backwards for 24 hours:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "o.seed_elements(lon=ilon, lat=ilat, radius=5000, radius_type='uniform', number=30000,\n time=end_time, object_type=object_type)\no.run(duration=duration, time_step=-900, time_step_output=3600, outfile=outfile)\nod = opendrift.open_xarray(outfile)\ndensity_backwards = od.get_histogram(pixelsize_m=5000).isel(time=-1).isel(origin_marker=0)\ndensity_backwards = density_backwards.where(density_backwards>0)\ndensity_backwards = density_backwards/density_backwards.sum()*100\nvmax = density_backwards.max()\no.plot(background=density_backwards, clabel='Probability of origin [%]', text=text, corners=corners,\n fast=True, markersize=.5, lalpha=.02, vmin=0, vmax=vmax)\nos.remove(outfile)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simulating then forwards, starting at a uniform grid 24 hours earlier (440 x 320 = 140800 elements at ~500m separation)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "o = Leeway(loglevel=50)\no.add_reader([reader_norkyst, reader_arome])\no.seed_elements(lon=lons, lat=lats, radius=0,\n time=start_time, object_type=object_type)\no.run(duration=duration, time_step=900, time_step_output=3600, outfile=outfile)\nprint(o)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finding the elements actually hitting the target (within 5 km) after 24 hours:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "lon, lat = o.get_lonlats()\nlonend = lon[:, -1]\nlatend = lat[:, -1]\ngeod = pyproj.Geod(ellps='WGS84')\non = np.ones(lonend.shape)\ndummy1, dummy2, dist2incident = geod.inv(lonend, latend, ilon*on, ilat*on)\nhits = np.where(dist2incident<5000)[0]\nhit_start_lons = lon[hits, 0]\nhit_start_lats = lat[hits, 0]\no_hit = opendrift.open(outfile, elements=hits)\n\no.animation(compare=o_hit, legend=['Elements not hitting target', 'Elements hitting target'],\n fast=True, corners=corners, text=text)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "o.plot(compare=o_hit, legend=['Elements not hitting target', 'Elements hitting target'],\n show_elements=False, fast=True, corners=corners, text=text)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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)\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "of = opendrift.open_xarray(outfile, elements=hits)\ndensity_forwards = of.get_histogram(pixelsize_m=5000).isel(time=0).isel(origin_marker=0)\ndensity_forwards = density_forwards.where(density_forwards>0)\nratio = density_forwards/density_forwards.sum()*100\no_hit.plot(background=ratio, clabel='Probability of origin [%]', text=text, corners=corners,\n fast=True, markersize=.5, lalpha=.02, vmin=0, vmax=vmax)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.6" } }, "nbformat": 4, "nbformat_minor": 0 }