{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Horse collar data exploration\n", "\n", "\n", "\n", "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/movingpandas/movingpandas-examples/main?filepath=2-analysis-examples/horse-collar.ipynb)\n", "[![IPYNB](https://img.shields.io/badge/view-ipynb-hotpink)](https://github.com/movingpandas/movingpandas-examples/blob/main/2-analysis-examples/horse-collar.ipynb)\n", "[![HTML](https://img.shields.io/badge/view-html-green)](https://movingpandas.github.io/movingpandas-website/2-analysis-examples/horse-collar.html)\n", "\n", "This notebook presents a systematic movement data exploration workflow. \n", "The workflow consists of five main steps:\n", "\n", "1. **Establishing an overview** by visualizing raw input data records\n", "2. **Putting records in context** by exploring information from consecutive movement data records (such as: time between records, speed, and direction)\n", "3. **Extracting trajectories, locations & events** by dividing the raw continuous tracks into individual trajectories, locations, and events\n", "4. **Exploring patterns** in trajectory and event data by looking at groups of the trajectories or events\n", "5. **Analyzing outliers** by looking at potential outliers and how they may challenge preconceived assumptions about the dataset characteristics\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import geopandas as gpd\n", "import movingpandas as mpd\n", "import shapely as shp\n", "import hvplot.pandas \n", "import matplotlib.pyplot as plt\n", "\n", "from pyproj import CRS\n", "from geopandas import GeoDataFrame, read_file\n", "from shapely.geometry import Point, LineString, Polygon\n", "from datetime import datetime, timedelta\n", "from holoviews import opts, dim, Layout\n", "from os.path import exists\n", "from urllib.request import urlretrieve\n", "\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "\n", "plot_defaults = {'linewidth':5, 'capstyle':'round', 'figsize':(9,3), 'legend':True}\n", "opts.defaults(opts.Overlay(active_tools=['wheel_zoom'], frame_width=300, frame_height=400))\n", "hvplot_defaults = {'tiles':None, 'cmap':'Viridis', 'colorbar':True}\n", "\n", "mpd.show_versions()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "mpd.__version__" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Raw data import" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df = read_file('../data/horse_collar.gpkg')\n", "df['t'] = pd.to_datetime(df['timestamp'])\n", "df = df.set_index('t').tz_localize(None)\n", "print(\"This dataset contains {} records.\\nThe first lines are:\".format(len(df)))\n", "df.head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df.columns" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df = df.drop(columns=['LMT_Date', 'LMT_Time',\n", " 'Origin', 'SCTS_Date', 'SCTS_Time', 'Latitude [?]', 'Longitude [?]',\n", " 'FixType', 'Main [V]', 'Beacon [V]', 'Sats', 'Sat',\n", " 'C/N', 'Sat_1', 'C/N_1', 'Sat_2', 'C/N_2', 'Sat_3', 'C/N_3', 'Sat_4',\n", " 'C/N_4', 'Sat_5', 'C/N_5', 'Sat_6', 'C/N_6', 'Sat_7', 'C/N_7', 'Sat_8',\n", " 'C/N_8', 'Sat_9', 'C/N_9', 'Sat_10', 'C/N_10', 'Sat_11', 'C/N_11',\n", " 'Easting', 'Northing',], axis=1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df.head()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "collar_id = df['CollarID'].unique()[0]\n", "print(\"There is only one collar with ID {}.\".format(collar_id))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df['Activity'].unique()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "original_crs = df.crs\n", "original_crs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Establishing an overview\n", "\n", "The first step in our proposed EDA workflow can be performed directly on raw\n", "input data since it does not require temporally ordered data. It is therefore suitable\n", "as a first exploratory step when dealing with new data.\n", "\n", "### Q1.1 Geographic extent: Is the geographical extent as expected and are there holes in the spatial coverage?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df.to_crs({'init': 'epsg:4326'}).hvplot(title='Geographic extent of the dataset', geo=True, tiles='OSM')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The main area (the horste's pasture?) is located south of Nykobing Strandhuse.\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pd.DataFrame(df).sort_values('lat').tail(2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "If we remove these first two records from the dataset, the remainder of the records are located in a small area:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df = df[2:].to_crs({'init': 'epsg:4326'})\n", "( df.hvplot(title='OSM showing paths and fences', size=2, geo=True, tiles='OSM') +\n", " df.hvplot(title='Imagery showing land cover details', size=2, color='red', geo=True, tiles='EsriImagery') )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It looks like the horse generally avoids areas without green vegetation since point patterns in these area appear more sparse than in other areas." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "temp = df.to_crs(CRS(25832))\n", "temp['geometry'] = temp['geometry'].buffer(5)\n", "total_area = temp.dissolve(by='CollarID').area \n", "total_area = total_area[collar_id]/10000\n", "print('The total area covered by the data is: {:,.2f} ha'.format(total_area))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q1.2 Temporal extent: Is the temporal extent as expected and are there holes in the temporal coverage?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"The dataset covers the time between {} and {}.\".format(df.index.min(), df.index.max()))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(\"That's {}\".format(df.index.max() - df.index.min()))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df['No'].resample('1d').count().hvplot(title='Number of records per day')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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). \n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q1.3 Spatio-temporal gaps: Does the geographic extent vary over time or do holes appear during certain times? \n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "df['Y-M'] = df.index.to_period('M') \n", "Layout([df[df['Y-M']==i].hvplot(title=str(i), size=2, geo=True, tiles='OSM') for i in df['Y-M'].unique()])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The largest change between months seems to be that the southernmost part of the pasture wasn't used in August and September 2019. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Putting records in context\n", "\n", "The second exploration step puts movement records in their temporal and geographic\n", "context. The exploration includes information based on consecutive movement data\n", "records, such as time between records (sampling intervals), speed, and direction.\n", "Therefore, this step requires temporally ordered data. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q2.1 Sampling intervals: Is the data sampled at regular or irregular intervals?\n", "\n", "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 ." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "t = df.reset_index().t\n", "df = df.assign(delta_t=t.diff().values)\n", "df['delta_t'] = df['delta_t'].dt.total_seconds()/60\n", "pd.DataFrame(df).hvplot.hist('delta_t', title='Histogram of intervals between consecutive records (in minutes)', bins=60, bin_range=(0, 60))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The time delta between consecutive records is usually around 30 minutes. \n", "\n", "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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q2.2 Speed values: Are there any unrealistic movements? \n", "\n", "For example: Does the data contain unattainable speeds?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tc = mpd.TrajectoryCollection(df, 'CollarID')\n", "traj = tc.trajectories[0]\n", "traj.add_speed()\n", "max_speed = traj.df.speed.max()\n", "print(\"The highest computed speed is {:,.2f} m/s ({:,.2f} km/h)\".format(max_speed, max_speed*3600/1000))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q2.3 Movement patterns: Are there any patterns in movement direction or speed?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pd.DataFrame(traj.df).hvplot.hist('speed', title='Histogram of speeds (in meters per second)', bins=90)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The speed distribution shows no surprising patterns." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "traj.add_direction(overwrite=True)\n", "pd.DataFrame(traj.df).hvplot.hist('direction', title='Histogram of directions', bins=90)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is some variation in movement directions but no directions stand out in the histogram.\n", "\n", "Let's look at spatial patterns of direction and speed!\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q2.4 Temporal context: Does the movement make sense in its temporal context? \n", "\n", "For example: Do nocturnal animal tracks show movement at night?\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pd.DataFrame(traj.df).hvplot.heatmap(title='Mean speed by hour of day and month of year', \n", " x='t.hour', y='t.month', C='speed', reduce_function=np.mean)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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. \n", "\n", "#### Temperature context\n", "\n", "In addition to time, the dataset also contains temperature information for each record:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "traj.df['n'] = 1\n", "pd.DataFrame(traj.df).hvplot.heatmap(title='Record count by temperature and month of year', \n", " x='Temp [?C]', y='t.month', C='n', reduce_function=np.sum)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pd.DataFrame(traj.df).hvplot.heatmap(title='Mean speed by temperature and month of year', \n", " x='Temp [?C]', y='t.month', C='speed', reduce_function=np.mean)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Q2.5 Geographic context: Does the movement make sense in its geographic context? \n", "\n", "For example: Do vessels follow traffic separation schemes defined in maritime maps? Are there any ship trajectories crossing land?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "traj.df['dir_class'] = ((traj.df['direction']-22.5)/45).round(0)\n", "temp = traj.df\n", "Layout([temp[temp['dir_class']==i].hvplot(geo=True, tiles='OSM', size=2, width=300, height=300, title=str(int(i*45))+\"°\") for i in sorted(temp['dir_class'].unique())])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are no obvious spatial movement direction patterns." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "traj.df['speed_class'] = (traj.df['speed']*2).round(1)\n", "temp = traj.df\n", "plots = []\n", "for i in sorted(temp['speed_class'].unique()):\n", " filtered = temp[temp['speed_class']==i]\n", " if len(filtered) < 10: \n", " continue\n", " plots.append(filtered.hvplot(geo=True, tiles='EsriImagery', color='red', size=2, width=300, height=300, title=str(i/2)))\n", "Layout(plots)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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. \n", "\n", "Medium speed records (classes 0.1 and 0.15 m/s) seem to be more common along paths and channels. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Extracting trajectories & locations / events\n", "\n", "The third exploration step looks at individual trajectories. It therefore requires that\n", "the continuous tracks are split into individual trajectories. Analysis results depend on\n", "how the continuous streams are divided into trajectories, locations, and events. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Trajectory lines: Do the trajectory lines look plausible or are there indications of out of sequence positions or other unrealistic location jumps?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tc.hvplot() " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Due to the 30 minute reporting interval, the trajectories are rather sparse. \n", "\n", "The trajectories mostly stay within the (fenced?) area. However, there are a few cases of positions outside the area." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Movement during week #1" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "daily = mpd.TemporalSplitter(tc).split(mode='day')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Layout([daily.trajectories[i].hvplot(title=daily.trajectories[i].id, c='speed', line_width=2, cmap='RdYlBu') for i in range(0,7)])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 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?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "daily_starts = daily.get_start_locations()\n", "daily_starts.set_index(pd.to_datetime(daily_starts['timestamp']), inplace=True) \n", "daily_starts['month'] = daily_starts.index.month\n", "daily_starts.hvplot(c='month', geo=True, tiles='EsriImagery', cmap='autumn')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is no clear preference for a certain home location where the horse would tend to spend the night. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "moving = mpd.TrajectoryCollection(traj.df[traj.df['speed'] > 0.05], 'CollarID')\n", "moving = mpd.ObservationGapSplitter(moving).split(gap=timedelta(minutes=70))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "moving.get_start_locations().hvplot(c='month', geo=True, tiles='EsriImagery', color='red')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.3 Trajectory length" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "daily_lengths = [traj.get_length() for traj in daily]\n", "daily_t = [traj.get_start_time() for traj in daily]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "daily_lengths = pd.DataFrame(daily_lengths, index=daily_t, columns=['length'])\n", "daily_lengths.hvplot(title='Daily trajectory length')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.)\n", "\n", "The seasonal trends agree well with the previously discovered seasonal movement speed patterns: winter trajectories tend to be shorter than summer trajectories. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.4 Covered area" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Method 1: Convex hulls around trajectory" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "daily_areas = [(traj.id, traj.to_crs(CRS(25832)).to_linestring().convex_hull.area/10000) for traj in daily]\n", "daily_areas = pd.DataFrame(daily_areas, index=daily_t, columns=['id', 'area'])\n", "daily_areas.hvplot(title='Daily covered area [ha]', y='area')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Method 2: Buffered trajectory" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "daily_areas = [(traj.id, traj.to_crs(CRS(25832)).to_linestring().buffer(15).area/10000) for traj in daily]\n", "daily_areas = pd.DataFrame(daily_areas, index=daily_t, columns=['id', 'area'])\n", "daily_areas.hvplot(title='Daily covered area [ha]', y='area')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The ten smallest areas are:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "daily_areas.sort_values(by='area')[:10]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "daily_areas = daily_areas.drop(datetime(2018,11,14,12,30,8))\n", "daily_areas = daily_areas.drop(datetime(2019,11,7,0,0,9))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The smallest area for a complete day was observed on 2018-11-19 with only 1.2 ha:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Layout([daily.get_trajectory(i).hvplot(title=i, c='speed', line_width=2, cmap='RdYlBu', width=300, height=300) for i in daily_areas.sort_values(by='area')[:3].id])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.5 Stop detection\n", "\n", "Instead of splitting the continuous track into daily trajectories, an alternative approach is to split it at stops. Stops can be defined as parts of the track where the moving object stays within a small area for a certain duration. \n", "\n", "Let's have a look at movement of one day and how stop detection parameter settings affect the results:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "MAX_DIAMETER = 100\n", "MIN_DURATION = timedelta(hours=3)\n", "\n", "one_day = daily.get_trajectory('30788_2018-11-17 00:00:00')\n", "one_day_stops = mpd.TrajectoryStopDetector(one_day).get_stop_segments(\n", " min_duration=MIN_DURATION, max_diameter=MAX_DIAMETER)\n", "\n", "( one_day.hvplot(title='Stops in Trajectory {}'.format(one_day.id), line_width=3.0, color='slategray') * \n", " one_day_stops.hvplot(line_width=5, tiles=None, color='deeppink') *\n", " one_day_stops.get_start_locations().hvplot(geo=True, size=200, color='deeppink') )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's apply stop detection to the whole dataset:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "stops = mpd.TrajectoryStopDetector(tc).get_stop_points(min_duration=MIN_DURATION, max_diameter=MAX_DIAMETER)\n", "len(stops)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The spatial distribution reveals preferred stop locations:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "stops.hvplot(geo=True, tiles='OSM', color='deeppink', size=MAX_DIAMETER, alpha=0.2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "stops['duration_h'] = (stops['end_time']-stops['start_time']).dt.total_seconds() / 3600\n", "pd.DataFrame(stops)['duration_h'].hvplot.hist(title='Stop duration histogram', xlabel='Duration [hours]', ylabel='n', bins=12)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Continue exploring MovingPandas\n", "\n", "1. [Bird migration analysis](bird-migration.ipynb)\n", "1. [Ship data analysis](ship-data.ipynb)\n", "1. [Horse collar data exploration](horse-collar.ipynb)\n", "1. [OSM traces](osm-traces.ipynb)\n", "1. [Soccer game](soccer-game.ipynb)\n", "1. [Mars rover & heli](mars-rover.ipynb)" ] } ], "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.10.8" } }, "nbformat": 4, "nbformat_minor": 4 }