Computing intersections with polygons

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

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.20.0

SYSTEM INFO
-----------
python     : 3.10.15 | packaged by conda-forge | (main, Oct 16 2024, 01:15:49) [MSC v.1941 64 bit (AMD64)]
executable : c:\Users\Agarkovam\AppData\Local\miniforge3\envs\mpd-ex\python.exe
machine    : Windows-10-10.0.19045-SP0

GEOS, GDAL, PROJ INFO
---------------------
GEOS       : None
GEOS lib   : None
GDAL       : None
GDAL data dir: None
PROJ       : 9.5.0
PROJ data dir: C:\Users\Agarkovam\AppData\Local\miniforge3\envs\mpd-ex\Library\share\proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 1.0.1
pandas     : 2.2.3
fiona      : None
numpy      : 1.23.1
shapely    : 2.0.6
pyproj     : 3.7.0
matplotlib : 3.9.2
mapclassify: 2.8.1
geopy      : 2.4.1
holoviews  : 1.20.0
hvplot     : 0.11.1
geoviews   : 1.13.0
stonesoup  : 1.4

gdf = read_file("../data/geolife_small.gpkg")
tc = mpd.TrajectoryCollection(gdf, "trajectory_id", t="t")

Clipping a Trajectory

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

ax = my_traj.plot()
gpd.GeoSeries(polygon).plot(ax=ax, color="lightgray")
intersections.plot(ax=ax, color="red", linewidth=5)

<Axes: >

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

Make this Notebook Trusted to load map: File -> Trust Notebook