Computing intersections with polygons¶
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.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
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: >
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