Skip to content

Closest Point of Approach #40

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Closest Point of Approach #40

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
38c91ed
cleanup and lint
SiggyF Apr 13, 2025
7ab67bb
remove all output
SiggyF Apr 13, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
269 changes: 269 additions & 0 deletions 2-analysis-examples/ship-cpa.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,269 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"id": "30f0a27a-919d-466f-9625-30d486ede07b",
"metadata": {},
"outputs": [],
"source": [
"import datetime\n",
"from typing import Iterator\n",
"\n",
"import folium\n",
"import geopandas as gpd\n",
"import hvplot.pandas\n",
"import matplotlib.colors\n",
"import matplotlib.pyplot as plt\n",
"import movingpandas as mpd\n",
"import movingpandas.cpa\n",
"import numpy as np\n",
"import pandas as pd\n",
"import shapely\n",
"from tqdm.notebook import tqdm"
]
},
{
"cell_type": "markdown",
"id": "460a21ed-9bee-4e8a-a862-b8cb443e9cf3",
"metadata": {},
"source": [
"# Closest point of approach example\n",
"This notebook demonstrates the CPA calculation for two trajectories. In maritime navigation, CPA stands for Closest Point of Approach. It is a concept used to assess the risk of collision between vessels. CPA refers to the minimum distance that will occur between two ships if they continue on their current courses and speeds. It is often calculated using data from the Automatic Identification System (AIS), which tracks vessel movements. Here the public Danish AIS dataset is used. \n",
"\n",
"There is an example where two ships meet on the evening of july 5th of 2017. The ships that approach eachother are MMSI: 265550210 and MMSI: 265410000. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7753cbbe-08e2-4f91-87f3-141bfccebeeb",
"metadata": {},
"outputs": [],
"source": [
"mmsi_a = 265550210\n",
"mmsi_b = 265410000\n",
"\n",
"color_a = (0, 0.3, 0.3)\n",
"color_b = (0.3, 0, 0.3)\n",
"t_0 = datetime.datetime(2017, 7, 5, 19, 0, 0)\n",
"t_1 = datetime.datetime(2017, 7, 5, 19, 13, 55)\n",
"\n",
"\n",
"# we expect to find this as the closest point and interval of overlapping time window\n",
"expected = {\n",
" \"t\": datetime.datetime(2017, 7, 5, 19, 8, 31),\n",
" \"p\": shapely.Point(669670.780244901, 6397464.37325368),\n",
" \"q\": shapely.Point(669566.587642144, 6397098.724455765),\n",
" \"t0_min\": datetime.datetime(2017, 7, 5, 19, 0, 1),\n",
" \"t1_max\": datetime.datetime(2017, 7, 5, 19, 13, 53),\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "db13d8d3-9e00-4206-9568-f3f060274b9b",
"metadata": {},
"source": [
"# Read the data\n",
"Here we read the data from a geopackage file. We convert the timestamps to a real date object. \n",
"For visualisation purposes we add a variable `t_min_t_0` so that we can visualize the time since start of the trajectories. \n",
"We convert the data to UTM 32N so that we can do our computations in Euclidean space, in meters.\n",
"\n",
"We'll subset the data using the 2 ship identifiers and our time window of interest. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "baa1b901-c744-43fc-9dfe-4f98853636d3",
"metadata": {},
"outputs": [],
"source": [
"ais_gdf = gpd.read_file(\"../data/ais.gpkg\")\n",
"ais_gdf[\"t\"] = pd.to_datetime(ais_gdf[\"Timestamp\"], format=\"%d/%m/%Y %H:%M:%S\")\n",
"ais_gdf[\"t_min_t_0\"] = (ais_gdf[\"t\"] - t_0).dt.total_seconds()\n",
"\n",
"# convert to local UTM so we can do metric calculations\n",
"ais_gdf = ais_gdf.to_crs(\"EPSG:32632\")\n",
"ais_gdf = ais_gdf.sort_values([\"MMSI\", \"t\"])"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "f2c4b3cb-36c1-44f3-8770-2ba267354df3",
"metadata": {},
"outputs": [],
"source": [
"# Filter on MMSI of interest and time window of interest\n",
"ais_gdf = ais_gdf[\n",
" np.logical_and.reduce(\n",
" [\n",
" ais_gdf[\"MMSI\"].isin([mmsi_a, mmsi_b]),\n",
" ais_gdf[\"t\"] >= t_0,\n",
" ais_gdf[\"t\"] < t_1,\n",
" ]\n",
" )\n",
"]"
]
},
{
"cell_type": "markdown",
"id": "7fb044db-4631-4718-a046-c9d1ce6c8bf8",
"metadata": {},
"source": [
"# Convert to trajectories\n",
"Here we convert our data into trajectories. We'll extract the 2 trajectories as seperate trajectories. \n",
"The CPA calculation uses a pair of trajectories as input. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6be51834-5ee6-447f-ba79-804f174f501a",
"metadata": {},
"outputs": [],
"source": [
"ais_gdf[\"traj_id\"] = ais_gdf[\"MMSI\"]\n",
"ais_gdf[\"obj_id\"] = ais_gdf[\"MMSI\"]\n",
"traj_collection = mpd.TrajectoryCollection(\n",
" ais_gdf.set_index(\"t\"), traj_id_col=\"traj_id\", obj_id_col=\"obj_id\", min_length=100\n",
")\n",
"traj_collection.add_speed(overwrite=True, units=(\"m\", \"s\"))\n",
"traj_collection.add_direction(overwrite=True)\n",
"\n",
"\n",
"traj_a = traj_collection.get_trajectory(mmsi_a)\n",
"traj_b = traj_collection.get_trajectory(mmsi_b)"
]
},
{
"cell_type": "markdown",
"id": "f34c1260-b1df-4c4d-a08a-a3b3a4210e4e",
"metadata": {},
"source": [
"# Visualisation\n",
"We'll create a visualization using the `hvplot` code from holoviews. This shows that the trajectories overlapped. \n",
"But how close did the ships get to eachother? That is what we'll compute here. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e4a69b49-c978-49a3-a30b-4c07eada1d1f",
"metadata": {},
"outputs": [],
"source": [
"# create a variant in wgs84 to show a plot\n",
"traj_collection_wgs84 = mpd.TrajectoryCollection(\n",
" ais_gdf.to_crs(\"EPSG:4326\").set_index(\"t\"),\n",
" traj_id_col=\"traj_id\",\n",
" obj_id_col=\"obj_id\",\n",
" min_length=100,\n",
")\n",
"# show trajectories with number of seconds since t_0\n",
"traj_collection_wgs84.hvplot(c=\"t_min_t_0\", cmap=\"Blues\")"
]
},
{
"cell_type": "markdown",
"id": "76c90f0e-1451-47bf-b237-85dfd874f38e",
"metadata": {},
"source": [
"# Calculate the CPA\n",
"Now we have the functions defined we can compute a dataframe of all the cpa calculations for each sequential time step pair in the timesteps of ship a and b. The table has a lot of variables that we'll probably remove. We'll probably keep the position of the closest point (actualy a line, pq), the distance (dist) and the time (t) and maybe the status or a code. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5075ff5d-9222-4b49-814b-c436b9c7821c",
"metadata": {},
"outputs": [],
"source": [
"cpa_calc = movingpandas.cpa.CPACalculator(traj_a, traj_b)\n",
"cpa_df = cpa_calc.segments_gdf()\n",
"cpa = cpa_calc.min()\n",
"cpa_df.sort_values(\"dist\").head()"
]
},
{
"cell_type": "markdown",
"id": "bdaed393-b497-49b1-92a9-a76bc943e659",
"metadata": {},
"source": [
"# Visualisation\n",
"This visualization shows what has been computed. All the blue lines show the segments for which a CPA has been computed. \n",
"Between these sets a red line was found (using interpolation) where the ships were closest to eachother (assuming their current course, so not just in distance). This red line is what is referred to as the Closest Point of Approach. The distance was 380m and the time of this CPA was 2017-07-05 19:08:31. We'll also visualize the line of closest point of approach using a holoviews map."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8027c913-1bdd-4f3b-b2d0-9d4e97b9d24b",
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
"\n",
"ax.set_xlabel(\"x [m, UTM32N]\")\n",
"ax.set_ylabel(\"y [m, UTM32N]\")\n",
"ax.set_title(\n",
" f\"Closest point of approach\\ndistance: {cpa.dist:.0f}m\\nMMSI: {mmsi_a} & {mmsi_b}\\n{cpa.t_at}\"\n",
")\n",
"\n",
"# compute colors ranging from\n",
"N = matplotlib.colors.Normalize(0, cpa_df[\"dist\"].max())\n",
"cmap = matplotlib.cm.RdYlGn\n",
"colors = cmap(N(cpa_df[\"dist\"]))\n",
"thin = 15\n",
"# don't show all the computed lines\n",
"cpa_df[::thin].plot(ax=ax, color=colors[::thin])\n",
"# show the computed closest point of approach\n",
"gpd.GeoSeries(cpa[\"geometry\"]).plot(ax=ax, color=\"black\", linestyle=\":\")\n",
"traj_a.plot(ax=ax)\n",
"traj_b.plot(ax=ax);"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "79146a33-d758-4c5e-a9dc-b2df5d5f5c5b",
"metadata": {},
"outputs": [],
"source": [
"# We can alos add the Closest Point of Approach to our map.\n",
"cpa_wgs84 = gpd.GeoDataFrame(\n",
" geometry=gpd.GeoSeries(cpa[\"geometry\"]), crs=\"EPSG:32632\"\n",
").to_crs(\"EPSG:3857\")\n",
"# show trajectories with with CPA line\n",
"(\n",
" traj_collection_wgs84.hvplot(c=\"t_min_t_0\", cmap=\"Blues\")\n",
" * cpa_wgs84.hvplot(line_color=(0.8, 0.3, 0.3), line_width=3, line_alpha=0.8)\n",
")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.12.10"
}
},
"nbformat": 4,
"nbformat_minor": 5
}