This tutorial uses AIS data published by the Danish Maritime Authority. The AIS record sample extracted for this tutorial covers vessel traffic on the 5th July 2017 near Gothenburg.
This tutorial covers:
%matplotlib inline
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:100% !important; }</style>"))
import urllib
import os
import numpy as np
import pandas as pd
from geopandas import GeoDataFrame, read_file
from shapely.geometry import Point, LineString, Polygon
from datetime import datetime, timedelta
import matplotlib.pyplot as plt
import sys
sys.path.append("..")
import movingpandas as mpd
import warnings
warnings.simplefilter("ignore")
%%time
df = read_file('data/demodata_ais.gpkg')
wgs84 = df.crs
print("Finished reading {}".format(len(df)))
Let's see what the data looks like:
df.head()
df.plot()
To convert the DataFrame to Trajectories we need to create a temporal index:
df['t'] = pd.to_datetime(df['Timestamp'], format='%d/%m/%Y %H:%M:%S')
df = df.set_index('t')
If we look at the data distributions, we can see that there are a lot of records with speed over ground (SOG) values of zero in this dataframe:
df['SOG'].hist(bins=100, figsize=(15,3))
Let's get rid of these rows with zero SOG:
print("Original size: {} rows".format(len(df)))
df = df[df.SOG>0]
print("Reduced to {} rows after removing 0 speed records".format(len(df)))
df['SOG'].hist(bins=100, figsize=(15,3))
Let's see what kind of ships we have in our dataset:
df['ShipType'].value_counts().plot(kind='bar', figsize=(15,3))
Finally, let's create trajectories:
%%time
MIN_LENGTH = 100 # meters
traj_collection = mpd.TrajectoryCollection(df, 'MMSI', min_length=MIN_LENGTH)
print("Finished creating {} trajectories".format(len(traj_collection)))
traj_collection = mpd.MinTimeDeltaGeneralizer(traj_collection).generalize(tolerance=timedelta(minutes=1))
Let's give the most common ship types distinct colors. The remaining ones will be just grey:
shiptype_to_color = {'Passenger': 'blue', 'HSC': 'green', 'Tanker': 'red', 'Cargo': 'orange'}
traj_collection.plot(column='ShipType', column_to_color=shiptype_to_color, linewidth=1, capstyle='round')
passenger = traj_collection.filter('ShipType', 'Passenger')
passenger.hvplot(title='Passenger ferries', line_width=2)
We can also plot individual trajectories to better visualize their properties, such as the changes in NavStatus:
my_traj = traj_collection.trajectories[0]
my_traj.df.head()
my_traj.df.tail()
my_traj.hvplot(title='Trajectory {}'.format(str(my_traj.id)), height=300, line_width=5.0, c='NavStatus', cmap='Dark2')
We can find ships passing under the bridge based on trajectory intersections with the bridge area.
area_of_interest = Polygon([(11.89935, 57.69270), (11.90161, 57.68902), (11.90334, 57.68967), (11.90104, 57.69354), (11.89935, 57.69270)])
intersecting = traj_collection.get_intersecting(area_of_interest)
print("Found {} intersections".format(len(intersecting)))
bridge_traj = intersecting.trajectories[0]
bridge_traj.hvplot(title='Trajectory {}'.format(str(bridge_traj.id)), height=300, line_width=5.0, c='NavStatus', cmap='Dark2')
bridge_traj.df.head()
Since AIS records with a speed over ground (SOG) value of zero have been removed from the dataset, we can use the split_by_observation_gap()
function to split the continuous observations into individual trips:
trips = mpd.ObservationGapSplitter(passenger).split(gap=timedelta(minutes=5))
print("Extracted {} individual trips from {} continuous vessel tracks".format(len(trips), len(passenger)))
Let's plot the resulting trips!
trips.hvplot(title='Passenger ferry trips', line_width=2)
Compared to plotting the original continuous observations, this visualization is much cleaner because there are no artifacts at the border of the area of interest.
Next, let's get the trip origins:
origins = trips.get_start_locations()
origins.hvplot(title='Trip origins by ship type', c='Name', geo=True, tiles='OSM')
In our data sample, trip origins can be:
origins.hvplot(title='Origins by speed', c='SOG', geo=True, tiles='OSM')
trips = mpd.ObservationGapSplitter(traj_collection).split(gap=timedelta(minutes=5))
area_of_interest = Polygon([(11.86815, 57.68273), (11.86992, 57.68047), (11.87419, 57.68140), (11.87288, 57.68348), (11.86815, 57.68273)])
We can identify vessels that start their trip within a given area of interest by intersecting trip starting locations with our area of interest:
departures = []
for traj in trips:
if traj.get_start_location().intersects(area_of_interest) and traj.get_length() > 100:
departures.append(traj)
print("Found {} departures".format(len(departures)))
departures[1].hvplot(title='Trajectory {}'.format(departures[1].id), line_width=5, c='Name', cmap='Dark2')
Let's see what kind of ships depart from here:
for traj in departures:
print("{} vessel '{}' departed at {}".format(traj.df['ShipType'].iloc[0], traj.df['Name'].iloc[0], traj.get_start_time()))
Of course, the same works for arrivals:
arrivals = []
for traj in trips:
if traj.get_end_location().intersects(area_of_interest) and traj.get_length() > 100:
arrivals.append(traj)
print("Found {} arrivals".format(len(arrivals)))
for traj in arrivals:
print("{} vessel '{}' arrived at {}".format(traj.df['ShipType'].iloc[0], traj.df['Name'].iloc[0], traj.get_end_time()))
To run this section, you need to have the scikit-learn package installed.
from sklearn.cluster import DBSCAN
from geopy.distance import great_circle
from shapely.geometry import MultiPoint
origins = trips.get_start_locations()
origins['lat'] = origins.geometry.y
origins['lon'] = origins.geometry.x
matrix = origins[['lat','lon']].values
kms_per_radian = 6371.0088
epsilon = 0.1 / kms_per_radian
db = DBSCAN(eps=epsilon, min_samples=1, algorithm='ball_tree', metric='haversine').fit(np.radians(matrix))
cluster_labels = db.labels_
num_clusters = len(set(cluster_labels))
clusters = pd.Series([matrix[cluster_labels == n] for n in range(num_clusters)])
print('Number of clusters: {}'.format(num_clusters))
origins['cluster'] = cluster_labels
def get_centermost_point(cluster):
centroid = (MultiPoint(cluster).centroid.x, MultiPoint(cluster).centroid.y)
centermost_point = min(cluster, key=lambda point: great_circle(point, centroid).m)
return Point(tuple(centermost_point)[1], tuple(centermost_point)[0])
centermost_points = clusters.map(get_centermost_point)
origins.hvplot(title='Clustered origins', c='cluster', geo=True, tiles='OSM', cmap='glasbey_dark')
origins_by_cluster = pd.DataFrame(origins).groupby(['cluster'])
summary = origins_by_cluster['ShipType'].unique().to_frame(name='types')
summary['n'] = origins_by_cluster.size()
summary['symbol_size'] = summary['n']*10 # for visualization purposes
summary['sog'] = origins_by_cluster['SOG'].mean()
summary['geometry'] = centermost_points
summary = summary[summary['n']>1].sort_values(by='n', ascending=False)
summary.head()
cluster_of_interest_id = 28
origins[origins['cluster']==cluster_of_interest_id].hvplot(
title='Cluster {}'.format(cluster_of_interest_id), c='ShipType', geo=True, tiles='OSM', height=500)
( trips.hvplot(title='Origin clusters by speed', color='gray', line_width=1) *
GeoDataFrame(summary, crs=wgs84).hvplot(c='sog', size='symbol_size', geo=True, cmap='RdYlGn')
)