Computing intersections with polygons¶

No description has been provided for this image

Binder IPYNB HTML

Clipping and intersection functions can be used to extract trajectory segments that are located within an area of interest polygon.

In [1]:
import pandas as pd
import geopandas as gpd
import movingpandas as mpd
import shapely as shp
import hvplot.pandas
import folium

from geopandas import GeoDataFrame, read_file
from shapely.geometry import Point, LineString, Polygon
from datetime import datetime, timedelta
from holoviews import opts

import warnings

warnings.filterwarnings("ignore")

opts.defaults(
    opts.Overlay(active_tools=["wheel_zoom"], frame_width=500, frame_height=400)
)

mpd.show_versions()

MovingPandas 0.22.4

SYSTEM INFO
-----------
python     : 3.10.19 | packaged by conda-forge | (main, Jan 26 2026, 23:45:08) [GCC 14.3.0]
executable : /home/anita/miniforge3/envs/mpd-ex/bin/python
machine    : Linux-6.8.0-107-generic-x86_64-with-glibc2.39

PROJ INFO
-----------
PROJ       : 9.6.2
PROJ data dir: /home/anita/miniforge3/envs/mpd-ex/share/proj

PYTHON DEPENDENCIES
-------------------
numpy      : 1.23.1
geopandas  : 1.0.1
geopy      : 2.4.1
geoviews   : 1.15.1
holoviews  : 1.22.1
hvplot     : 0.12.2
mapclassify: 2.8.1
matplotlib : 3.10.8
pandas     : 2.3.3
pyproj     : 3.7.1
shapely    : 2.1.2
stonesoup  : 1.8
In [2]:
gdf = read_file("../data/geolife_small.gpkg")
tc = mpd.TrajectoryCollection(gdf, "trajectory_id", t="t")

Clipping a Trajectory¶

In [3]:
xmin, xmax, ymin, ymax = 116.365035, 116.3702945, 39.904675, 39.907728
polygon = Polygon(
    [(xmin, ymin), (xmin, ymax), (xmax, ymax), (xmax, ymin), (xmin, ymin)]
)

my_traj = tc.trajectories[2]
intersections = my_traj.clip(polygon)

print("Found {} intersections".format(len(intersections)))

Found 1 intersections
In [4]:
ax = my_traj.plot()
gpd.GeoSeries(polygon).plot(ax=ax, color="lightgray")
intersections.plot(ax=ax, color="red", linewidth=5)
Out[4]:
<Axes: >
No description has been provided for this image
In [5]:
m = my_traj.explore(color="blue", style_kwds={"weight": 4}, name="Trajectory")

intersections.explore(m=m, color="red", style_kwds={"weight": 4}, name="Intersection")

folium.TileLayer("OpenStreetMap").add_to(m)
folium.LayerControl().add_to(m)

m
Out[5]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Clipping a TrajectoryCollection¶

Alternatively, using TrajectoryCollection:

In [6]:
clipped = tc.clip(polygon)
clipped
Out[6]:
TrajectoryCollection with 2 trajectories
In [7]:
clipped.plot()
Out[7]:
<Axes: >
No description has been provided for this image
In [8]:
clipped.explore(
    column="trajectory_id",
    cmap="cool",
    tiles="CartoDB positron",
    style_kwds={"weight": 4},
)
Out[8]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Computing intersections for a Trajectory¶

In [9]:
polygon_feature = {"geometry": polygon, "properties": {"field1": "abc"}}
In [10]:
my_traj = tc.trajectories[2]
intersections = my_traj.intersection(polygon_feature)
intersections
Out[10]:
TrajectoryCollection with 1 trajectories
In [11]:
intersections.plot()
Out[11]:
<Axes: >
No description has been provided for this image
In [12]:
intersections.explore(color="blue", style_kwds={"weight": 4})
Out[12]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [13]:
intersections.to_point_gdf()
Out[13]:
id sequence trajectory_id tracker geometry intersecting_field1
t
2009-02-04 10:43:04.222205 3225 773 3_0 2 POINT (116.36799 39.90468) abc
2009-02-04 10:43:05.000000 3226 774 3_0 2 POINT (116.36798 39.90471) abc
2009-02-04 10:43:07.000000 3227 775 3_0 2 POINT (116.36795 39.9048) abc
2009-02-04 10:43:09.000000 3228 776 3_0 2 POINT (116.36793 39.90487) abc
2009-02-04 10:43:11.000000 3229 777 3_0 2 POINT (116.36794 39.90495) abc
... ... ... ... ... ... ...
2009-02-04 10:48:14.000000 3390 938 3_0 2 POINT (116.36535 39.90589) abc
2009-02-04 10:48:15.000000 3391 939 3_0 2 POINT (116.36525 39.90589) abc
2009-02-04 10:48:16.000000 3392 940 3_0 2 POINT (116.36516 39.90589) abc
2009-02-04 10:48:17.000000 3393 941 3_0 2 POINT (116.36507 39.90589) abc
2009-02-04 10:48:17.399995 3393 941 3_0 2 POINT (116.36504 39.90589) abc

170 rows × 6 columns

Computing intersections for a TrajectoryCollection¶

In [14]:
intersections = tc.intersection(polygon_feature)
intersections
Out[14]:
TrajectoryCollection with 2 trajectories
In [15]:
intersections.plot()
Out[15]:
<Axes: >
No description has been provided for this image
In [16]:
intersections.explore(
    column="trajectory_id",
    cmap="autumn",
    tiles="CartoDB positron",
    style_kwds={"weight": 4},
)
Out[16]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [17]:
intersections.to_point_gdf()
Out[17]:
id sequence trajectory_id tracker geometry intersecting_field1
t
2009-02-04 10:43:04.222205 3225 773 3_0 2 POINT (116.36799 39.90468) abc
2009-02-04 10:43:05.000000 3226 774 3_0 2 POINT (116.36798 39.90471) abc
2009-02-04 10:43:07.000000 3227 775 3_0 2 POINT (116.36795 39.9048) abc
2009-02-04 10:43:09.000000 3228 776 3_0 2 POINT (116.36793 39.90487) abc
2009-02-04 10:43:11.000000 3229 777 3_0 2 POINT (116.36794 39.90495) abc
... ... ... ... ... ... ...
2009-03-10 11:08:05.000000 4934 672 4_0 2 POINT (116.36539 39.90588) abc
2009-03-10 11:08:06.000000 4935 673 4_0 2 POINT (116.36529 39.90588) abc
2009-03-10 11:08:07.000000 4936 674 4_0 2 POINT (116.36519 39.90588) abc
2009-03-10 11:08:08.000000 4937 675 4_0 2 POINT (116.36509 39.90588) abc
2009-03-10 11:08:08.575754 4937 675 4_0 2 POINT (116.36504 39.90588) abc

375 rows × 6 columns

In [ ]: