Bird migration analysis example¶
This tutorial uses data published on Movebank, specifically: Navigation experiments in lesser black-backed gulls (data from Wikelski et al. 2015)-gps.csv
This tutorial covers:
- Trajectory data preprocessing
- Loading movement data from common geospatial file formats
- Exploring spatial & non-spatial data distributions
- Converting GeoDataFrames into Trajectories describing continuous tracks of moving objects
- Trajectory data analysis
- Investigating individual trajectories
- Comparing different years
- Investigating trajectories of multiple individuals
In [1]:
import numpy as np
import pandas as pd
import geopandas as gpd
import movingpandas as mpd
import shapely as shp
import hvplot.pandas
import matplotlib.pyplot as plt
from geopandas import GeoDataFrame, read_file
from shapely.geometry import Point, LineString, Polygon
from datetime import datetime, timedelta
from holoviews import opts, dim
from os.path import exists
from urllib.request import urlretrieve
import warnings
warnings.filterwarnings('ignore')
plot_defaults = {'linewidth':5, 'capstyle':'round', 'figsize':(9,3), 'legend':True}
opts.defaults(opts.Overlay(active_tools=['wheel_zoom'], frame_width=300, frame_height=500))
hvplot_defaults = {'tiles':None, 'cmap':'Viridis', 'colorbar':True}
mpd.show_versions()
MovingPandas 0.17.0 SYSTEM INFO ----------- python : 3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:34:57) [MSC v.1936 64 bit (AMD64)] executable : H:\miniconda3\envs\mpd-ex\python.exe machine : Windows-10-10.0.19045-SP0 GEOS, GDAL, PROJ INFO --------------------- GEOS : None GEOS lib : None GDAL : 3.7.0 GDAL data dir: None PROJ : 9.2.1 PROJ data dir: H:\miniconda3\pkgs\proj-9.0.0-h1cfcee9_1\Library\share\proj PYTHON DEPENDENCIES ------------------- geopandas : 0.13.2 pandas : 2.0.3 fiona : 1.9.4 numpy : 1.24.4 shapely : 2.0.1 rtree : 1.0.1 pyproj : 3.6.0 matplotlib : 3.7.2 mapclassify: 2.5.0 geopy : 2.3.0 holoviews : 1.17.0 hvplot : 0.8.3 geoviews : 1.9.6 stonesoup : 1.0
Loading the bird movement data¶
In [2]:
%%time
df = read_file('../data/gulls.gpkg')
print(f"Finished reading {len(df)}")
Finished reading 89867 CPU times: total: 18.2 s Wall time: 18.5 s
This is what the data looks like:
In [3]:
df.head()
Out[3]:
| event-id | visible | timestamp | location-long | location-lat | sensor-type | individual-taxon-canonical-name | tag-local-identifier | individual-local-identifier | study-name | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1082620685 | true | 2009-05-27 14:00:00.000 | 24.58617 | 61.24783 | gps | Larus fuscus | 91732 | 91732A | Navigation experiments in lesser black-backed ... | POINT (24.58617 61.24783) |
| 1 | 1082620686 | true | 2009-05-27 20:00:00.000 | 24.58217 | 61.23267 | gps | Larus fuscus | 91732 | 91732A | Navigation experiments in lesser black-backed ... | POINT (24.58217 61.23267) |
| 2 | 1082620687 | true | 2009-05-28 05:00:00.000 | 24.53133 | 61.18833 | gps | Larus fuscus | 91732 | 91732A | Navigation experiments in lesser black-backed ... | POINT (24.53133 61.18833) |
| 3 | 1082620688 | true | 2009-05-28 08:00:00.000 | 24.58200 | 61.23283 | gps | Larus fuscus | 91732 | 91732A | Navigation experiments in lesser black-backed ... | POINT (24.58200 61.23283) |
| 4 | 1082620689 | true | 2009-05-28 14:00:00.000 | 24.58250 | 61.23267 | gps | Larus fuscus | 91732 | 91732A | Navigation experiments in lesser black-backed ... | POINT (24.58250 61.23267) |
In [4]:
df.plot()
Out[4]:
<Axes: >
Let's see how many individuals we have in the dataset:
In [5]:
df['individual-local-identifier'].unique()
Out[5]:
array(['91732A', '91733A', '91734A', '91735A', '91737A', '91738A',
'91739A', '91740A', '91741A', '91742A', '91743A', '91744A',
'91745A', '91746A', '91747A', '91748A', '91749A', '91750A',
'91751A', '91752A', '91754A', '91755A', '91756A', '91758A',
'91759A', '91761A', '91762A', '91763A', '91764A', '91765A',
'91766A', '91767A', '91769A', '91771A', '91774A', '91775A',
'91776A', '91777A', '91778A', '91779A', '91780A', '91781A',
'91782A', '91783A', '91785A', '91786A', '91787A', '91788A',
'91789A', '91794A', '91795A', '91797A', '91798A', '91799A',
'91800A', '91802A', '91803A', '91807A', '91809A', '91810A',
'91811A', '91812A', '91813A', '91814A', '91815A', '91816A',
'91819A', '91821A', '91823A', '91824A', '91825A', '91826A',
'91827A', '91828A', '91829A', '91830A', '91831A', '91832A',
'91835A', '91836A', '91837A', '91838A', '91839A', '91843A',
'91845A', '91846A', '91848A', '91849A', '91852A', '91854A',
'91858A', '91861A', '91862A', '91864A', '91865A', '91866A',
'91870A', '91871A', '91872A', '91875A', '91876A', '91877A',
'91878A', '91880A', '91881A', '91884A', '91885A', '91893A',
'91894A', '91897A', '91900A', '91901A', '91903A', '91907A',
'91908A', '91910A', '91911A', '91913A', '91916A', '91918A',
'91919A', '91920A', '91921A', '91924A', '91929A', '91930A'],
dtype=object)
The records per individual are not evenly distributed:
In [6]:
df['individual-local-identifier'].value_counts().plot(kind='bar', figsize=(17,3))
Out[6]:
<Axes: xlabel='individual-local-identifier'>
Finally, let's create trajectories:
In [7]:
tc = mpd.TrajectoryCollection(df, 'individual-local-identifier', t='timestamp', min_length=100)
tc
Out[7]:
TrajectoryCollection with 125 trajectories
And let's generalize them to speed up the following analyses:
In [8]:
tc = mpd.MinTimeDeltaGeneralizer(tc).generalize(tolerance=timedelta(days=1))
Investigating individual trajectories¶
Let's pick out a specific individual. For example, '91916A' is the individual with most records in our dataset:
In [9]:
filtered = tc.filter('individual-local-identifier', '91916A')
my_traj = filtered.trajectories[0].copy()
my_traj.df.head()
Out[9]:
| event-id | visible | location-long | location-lat | sensor-type | individual-taxon-canonical-name | tag-local-identifier | individual-local-identifier | study-name | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|
| timestamp | ||||||||||
| 2009-08-15 15:00:00 | 1082625177 | true | 7.91500 | 54.18533 | gps | Larus fuscus | 91916 | 91916A | Navigation experiments in lesser black-backed ... | POINT (7.91500 54.18533) |
| 2009-08-16 15:00:00 | 1082625181 | true | 9.44183 | 54.87233 | gps | Larus fuscus | 91916 | 91916A | Navigation experiments in lesser black-backed ... | POINT (9.44183 54.87233) |
| 2009-08-17 15:00:00 | 1082625185 | true | 9.44250 | 54.87283 | gps | Larus fuscus | 91916 | 91916A | Navigation experiments in lesser black-backed ... | POINT (9.44250 54.87283) |
| 2009-08-18 15:00:00 | 1082625189 | true | 9.41167 | 54.85550 | gps | Larus fuscus | 91916 | 91916A | Navigation experiments in lesser black-backed ... | POINT (9.41167 54.85550) |
| 2009-08-19 15:00:00 | 1082625193 | true | 9.39583 | 54.88000 | gps | Larus fuscus | 91916 | 91916A | Navigation experiments in lesser black-backed ... | POINT (9.39583 54.88000) |
In [10]:
my_traj.hvplot(title=f'Movement of {my_traj.id}', line_width=2)
Out[10]:
This individual has been travelling back and forth for quite a few years!
One way to take a closer look at this individual's travels is to split the overall track into yearly trips:
In [11]:
trips_by_year = mpd.TemporalSplitter(filtered).split(mode='year')
trips_by_year.to_traj_gdf()
Out[11]:
| individual-local-identifier | start_t | end_t | geometry | length | direction | |
|---|---|---|---|---|---|---|
| 0 | 91916A_2009-12-31 00:00:00 | 2009-08-15 15:00:00 | 2009-12-31 19:00:00 | LINESTRING (7.91500 54.18533, 9.44183 54.87233... | 5.352375e+06 | 131.939002 |
| 1 | 91916A_2010-12-31 00:00:00 | 2010-01-01 19:00:00 | 2010-12-31 07:00:00 | LINESTRING (39.18417 21.17583, 39.18833 21.170... | 1.232428e+07 | 281.047091 |
| 2 | 91916A_2011-12-31 00:00:00 | 2011-01-01 07:00:00 | 2011-12-31 04:00:00 | LINESTRING (39.17000 21.18017, 39.16883 21.156... | 1.441793e+07 | 189.238229 |
| 3 | 91916A_2012-12-31 00:00:00 | 2012-01-01 04:00:00 | 2012-12-31 19:00:00 | LINESTRING (39.16933 21.16667, 39.17567 21.142... | 1.324811e+07 | 62.132640 |
| 4 | 91916A_2013-12-31 00:00:00 | 2013-01-01 19:00:00 | 2013-12-31 13:00:00 | LINESTRING (39.18167 21.17333, 39.18100 21.173... | 1.293261e+07 | 273.165851 |
| 5 | 91916A_2014-12-31 00:00:00 | 2014-01-01 13:00:00 | 2014-12-31 19:00:00 | LINESTRING (39.17083 21.15400, 39.17100 21.157... | 1.300973e+07 | 33.742967 |
| 6 | 91916A_2015-12-31 00:00:00 | 2015-01-01 19:00:00 | 2015-08-27 09:00:00 | LINESTRING (39.18167 21.16967, 39.18233 21.169... | 6.551740e+06 | 343.887905 |
Now we can explore individual years:
In [12]:
one_year = trips_by_year.get_trajectory('91916A_2010-12-31 00:00:00')
one_year
Out[12]:
Trajectory 91916A_2010-12-31 00:00:00 (2010-01-01 19:00:00 to 2010-12-31 07:00:00) | Size: 354 | Length: 12324280.3m Bounds: (18.53417, 21.047, 39.203, 61.543) LINESTRING (39.18417 21.17583, 39.18833 21.17083, 39.18767 21.172, 39.1865 21.17117, 39.18817 21.170
In [13]:
one_year.add_speed(units=("km", "h"))
In [14]:
one_year.hvplot(title=f'Movement speed of {one_year.id}',
line_width=5.0, c='speed', cmap='RdYlGn', colorbar=True, clim=(0,20))
Out[14]:
Let's see where this individual was on a specific day:
In [15]:
def plot_location_at_timestamp(traj, t, fig_size=250):
loc = GeoDataFrame([traj.get_row_at(t)])
return (loc.hvplot(title=str(t), geo=True, tiles='OSM', size=200, color='red') *
traj.hvplot(line_width=1.0, color='black', tiles=False, width=fig_size, height=fig_size))
In [16]:
( plot_location_at_timestamp(one_year, datetime(2010,9,1)) +
plot_location_at_timestamp(one_year, datetime(2010,10,1)) +
plot_location_at_timestamp(one_year, datetime(2010,11,1)) )
Out[16]: