Tutorial 3: Horse collar data exploration¶
This notebook presents a systematic movement data exploration workflow. The proposed workflow consists of five main steps:
- Establishing an overview by visualizing raw input data records
- Putting records in context by exploring information from consecutive movement data records (such as: time between records, speed, and direction)
- Extracting trajectories, locations & events by dividing the raw continuous tracks into individual trajectories, locations, and events
- Exploring patterns in trajectory and event data by looking at groups of the trajectories or events
- Analyzing outliers by looking at potential outliers and how they may challenge preconceived assumptions about the dataset characteristics
The workflow is demonstrated using horse collar tracking data provided by Prof. Lene Fischer (University of Copenhagen) and the Center for Technology & Environment of Guldborgsund Municiplaity in Denmark but should be generic enough to be applied to other tracking datasets.
The workflow is implemented in Python using Pandas, GeoPandas, and MovingPandas (http://movingpandas.org).
For an interactive version of this notebook visit https://mybinder.org/v2/gh/anitagraser/movingpandas/master.
Setup¶
%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
import geopandas as gpd
from geopandas import GeoDataFrame, read_file
from datetime import datetime, timedelta
from pyproj import CRS
import sys
sys.path.append("..")
import movingpandas as mpd
import warnings
warnings.simplefilter("ignore")
import hvplot.pandas # seems to be necessary for the following import to work
from holoviews import opts
opts.defaults(opts.Overlay(active_tools=['wheel_zoom']))
Raw data import¶
df = read_file('data/demodata_horse_collar.gpkg')
df['t'] = pd.to_datetime(df['timestamp'])
df = df.set_index('t').tz_localize(None)
print("This dataset contains {} records.\nThe first lines are:".format(len(df)))
df.head()
df.columns
df = df.drop(columns=['LMT_Date', 'LMT_Time',
'Origin', 'SCTS_Date', 'SCTS_Time', 'Latitude [?]', 'Longitude [?]',
'FixType', 'Main [V]', 'Beacon [V]', 'Sats', 'Sat',
'C/N', 'Sat_1', 'C/N_1', 'Sat_2', 'C/N_2', 'Sat_3', 'C/N_3', 'Sat_4',
'C/N_4', 'Sat_5', 'C/N_5', 'Sat_6', 'C/N_6', 'Sat_7', 'C/N_7', 'Sat_8',
'C/N_8', 'Sat_9', 'C/N_9', 'Sat_10', 'C/N_10', 'Sat_11', 'C/N_11',
'Easting', 'Northing',], axis=1)
df.head()
collar_id = df['CollarID'].unique()[0]
print("There is only one collar with ID {}.".format(collar_id))
df['Activity'].unique()
original_crs = df.crs
original_crs
1. Establishing an overview¶
The first step in our proposed EDA workflow can be performed directly on raw input data since it does not require temporally ordered data. It is therefore suitable as a first exploratory step when dealing with new data.
Q1.1 Geographic extent: Is the geographical extent as expected and are there holes in the spatial coverage?¶
df.to_crs({'init': 'epsg:4326'}).hvplot(title='Geographic extent of the dataset', geo=True, tiles='OSM', width=500, height=500)
The main area (the horste's pasture?) is located south of Nykobing Strandhuse.
However, we find also find two records on the road north west of the main area. Both points have been recorded on 2018-11-14 which is the first day of the dataset.
pd.DataFrame(df).sort_values('lat').tail(2)
A potential hypothesis for the origin of these two records is that the horse (or the collar) was transported on 2018-11-14, taking the road from Nykobing Falster south to the pasture.
If we remove these first two records from the dataset, the remainder of the records are located in a small area:
df = df[2:].to_crs({'init': 'epsg:4326'})
( df.hvplot(title='OSM showing paths and fences', size=2, geo=True, tiles='OSM', width=500, height=500) +
df.hvplot(title='Imagery showing land cover details', size=2, color='red', geo=True, tiles='EsriImagery', width=500, height=500) )
It looks like the horse generally avoids areas without green vegetation since point patterns in these area appear more sparse than in other areas.
temp = df.to_crs(CRS(25832))
temp['geometry'] = temp['geometry'].buffer(5)
total_area = temp.dissolve(by='CollarID').area
total_area = total_area[collar_id]/10000
print('The total area covered by the data is: {:,.2f} ha'.format(total_area))
Q1.2 Temporal extent: Is the temporal extent as expected and are there holes in the temporal coverage?¶
print("The dataset covers the time between {} and {}.".format(df.index.min(), df.index.max()))
print("That's {}".format(df.index.max() - df.index.min()))
df['No'].resample('1d').count().hvplot(title='Number of records per day')
On most days there are 48 (+/- 1) records per day. However, there are some days with more records (in Nov 2018 and later between Mai and August 2019).
There is one gap: On 2019-10-18 there are no records in the dataset and the previous day only contains 37 and the following day 27 records.
Q1.3 Spatio-temporal gaps: Does the geographic extent vary over time or do holes appear during certain times?¶
Considering that the dataset covers a whole year, it may be worthwhile to look at the individual months using small multiples map plots, for example:
df['Y-M'] = df.index.to_period('M')
a = None
for i in df['Y-M'].unique():
plot = df[df['Y-M']==i].hvplot(title=str(i), size=2, geo=True, tiles='OSM', width=300, height=300)
if a: a = a + plot
else: a = plot
a
The largest change between months seems to be that the southernmost part of the pasture wasn't used in August and September 2019.
2. Putting records in context¶
The second exploration step puts movement records in their temporal and geographic context. The exploration includes information based on consecutive movement data records, such as time between records (sampling intervals), speed, and direction. Therefore, this step requires temporally ordered data.
Q2.1 Sampling intervals: Is the data sampled at regular or irregular intervals?¶
For example, tracking data of migratory animals is expected to exhibit seasonal changes. Such changes in vehicle tracking systems however may indicate issues with data collection .
t = df.reset_index().t
df = df.assign(delta_t=t.diff().values)
df['delta_t'] = df['delta_t'].dt.total_seconds()/60
pd.DataFrame(df).hvplot.hist('delta_t', title='Histogram of intervals between consecutive records (in minutes)', bins=60, bin_range=(0, 60))
The time delta between consecutive records is usually around 30 minutes.
However, it seems that sometimes the intererval has been decreased to around 15 minutes. This would explain why some days have more than the usual 48 records.
Q2.2 Speed values: Are there any unrealistic movements?¶
For example: Does the data contain unattainable speeds?
tc = mpd.TrajectoryCollection(df, 'CollarID')
traj = tc.trajectories[0]
traj.add_speed()
max_speed = traj.df.speed.max()
print("The highest computed speed is {:,.2f} m/s ({:,.2f} km/h)".format(max_speed, max_speed*3600/1000))
Q2.3 Movement patterns: Are there any patterns in movement direction or speed?¶
pd.DataFrame(traj.df).hvplot.hist('speed', title='Histogram of speeds (in meters per second)', bins=90)
The speed distribution shows no surprising patterns.
traj.add_direction(overwrite=True)
pd.DataFrame(traj.df).hvplot.hist('direction', title='Histogram of directions', bins=90)
There is some variation in movement directions but no directions stand out in the histogram.
Let's look at spatial patterns of direction and speed!
Q2.4 Temporal context: Does the movement make sense in its temporal context?¶
For example: Do nocturnal animal tracks show movement at night?
pd.DataFrame(traj.df).hvplot.heatmap(title='Mean speed by hour of day and month of year',
x='t.hour', y='t.month', C='speed', reduce_function=np.mean)
The movement speed by hour of day shows a clear pattern throughout the year with earlier and longer fast movements during the summer months and later and slower movements during the winter months.
Temperature context¶
In addition to time, the dataset also contains temperature information for each record:
traj.df['n'] = 1
pd.DataFrame(traj.df).hvplot.heatmap(title='Record count by temperature and month of year',
x='Temp [?C]', y='t.month', C='n', reduce_function=np.sum)
pd.DataFrame(traj.df).hvplot.heatmap(title='Mean speed by temperature and month of year',
x='Temp [?C]', y='t.month', C='speed', reduce_function=np.mean)
Q2.5 Geographic context: Does the movement make sense in its geographic context?¶
For example: Do vessels follow traffic separation schemes defined in maritime maps? Are there any ship trajectories crossing land?
traj.df['dir_class'] = ((traj.df['direction']-22.5)/45).round(0)
a = None
temp = traj.df
for i in sorted(temp['dir_class'].unique()):
plot = temp[temp['dir_class']==i].hvplot(geo=True, tiles='OSM', size=2, width=300, height=300, title=str(int(i*45))+"°")
if a: a = a + plot
else: a = plot
a
There are no obvious spatial movement direction patterns.
traj.df['speed_class'] = (traj.df['speed']*2).round(1)
a = None
temp = traj.df
for i in sorted(temp['speed_class'].unique()):
filtered = temp[temp['speed_class']==i]
if len(filtered) < 10:
continue
plot = filtered.hvplot(geo=True, tiles='EsriImagery', color='red', size=2, width=300, height=300, title=str(i/2)) # alpha=max(0.05, 50/len(filtered)),
if a: a = a + plot
else: a = plot
a
Low speed records (classes 0.0 and 0.05 m/s) are distributed over the whole area with many points on the outline (fence?) of the area.
Medium speed records (classes 0.1 and 0.15 m/s) seem to be more common along paths and channels.
3. Extracting trajectories & locations / events¶
The third exploration step looks at individual trajectories. It therefore requires that the continuous tracks are split into individual trajectories. Analysis results depend on how the continuous streams are divided into trajectories, locations, and events.
3.1 Trajectory lines: Do the trajectory lines look plausible or are there indications of out of sequence positions or other unrealistic location jumps?¶
tc.hvplot()
Due to the 30 minute reporting interval, the trajectories are rather sparse.
The trajectories mostly stay within the (fenced?) area. However, there are a few cases of positions outside the area.
Movement during week #1¶
daily = tc.split_by_date(mode='day')
a = None
for i in range(0,7):
if a: a = a + daily.trajectories[i].hvplot(title=daily.trajectories[i].id, c='speed', line_width=2, cmap='RdYlBu', width=300, height=300)
else: a = daily.trajectories[i].hvplot(title=daily.trajectories[i].id, c='speed', line_width=2, cmap='RdYlBu', width=300, height=300)
a
3.2 Home/depot locations: Do day trajectories start and end at the same home (for human and animal movement) or depot (for logistics applications) location?¶
daily_starts = daily.get_start_locations()
daily_starts['month'] = daily_starts.t.dt.month
daily_starts.hvplot(c='month', geo=True, tiles='EsriImagery', cmap='autumn', width=500, height=500)
There is no clear preference for a certain home location where the horse would tend to spend the night.
Instead of spliting by date, we can also specify a minimum movement speed and then split the continuous observation when this minimum speed is not reached for a certain time:
moving = mpd.TrajectoryCollection(traj.df[traj.df['speed'] > 0.05], 'CollarID')
moving = moving.split_by_observation_gap(timedelta(minutes=70))
moving.get_start_locations().hvplot(c='month', geo=True, tiles='EsriImagery', color='red', width=500, height=500)
3.3 Trajectory length¶
daily_lengths = [traj.get_length() for traj in daily]
daily_t = [traj.get_start_time() for traj in daily]
daily_lengths = pd.DataFrame(daily_lengths, index=daily_t, columns=['length'])
daily_lengths.hvplot(title='Daily trajectory length')
The length of the daily trajectories varies between 1.6 and 6.2 km. (It is worth noting that this has to be considered a lower bound of the movement due to the sparseness of the tracking data.)
The seasonal trends agree well with the previously discovered seasonal movement speed patterns: winter trajectories tend to be shorter than summer trajectories.
3.4 Covered area¶
Method 1: Convex hulls around trajectory¶
daily_areas = [(traj.id, traj.to_crs(CRS(25832)).to_linestring().convex_hull.area/10000) for traj in daily]
daily_areas = pd.DataFrame(daily_areas, index=daily_t, columns=['id', 'area'])
daily_areas.hvplot(title='Daily covered area [ha]', y='area')
Method 2: Buffered trajectory¶
daily_areas = [(traj.id, traj.to_crs(CRS(25832)).to_linestring().buffer(15).area/10000) for traj in daily]
daily_areas = pd.DataFrame(daily_areas, index=daily_t, columns=['id', 'area'])
daily_areas.hvplot(title='Daily covered area [ha]', y='area')
The ten smallest areas are:
daily_areas.sort_values(by='area')[:10]
The days with the smallest areas cover areas include the first and the last observation day (since they are only partially recorded). We can remove those:
daily_areas = daily_areas.drop(datetime(2018,11,14,12,30,8))
daily_areas = daily_areas.drop(datetime(2019,11,7,0,0,9))
The smallest area for a complete day was observed on 2018-11-19 with only 1.2 ha:
a = None
for i in daily_areas.sort_values(by='area')[:8].id:
traj = daily.get_trajectory(i)
if a: a = a + traj.hvplot(title=i, c='speed', line_width=2, cmap='RdYlBu', width=300, height=300)
else: a = traj.hvplot(title=i, c='speed', line_width=2, cmap='RdYlBu', width=300, height=300)
a
4. Exploring patterns in trajectory and event data¶
The fourth exploration step looks at the set of all or subsets of the trajectories and events extracted from the movement data. This step commonly involves many of the computationally more expensive operations, such as trajectory similarity computations.
TODO
5. Analyzing outliers¶
The fifth and final exploration steps looks at potential outliers and how they may challenge preconceived assumptions about the dataset characteristics. This step in particular regularly requires domain knowledge that exceeds the information that can be discovered from the movement data itself.
TODO
Conclusion¶
TODO