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.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: >
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_2009-02-04 10:43:04.222205 2 POINT (116.36799 39.90468) abc
2009-02-04 10:43:05.000000 3226 774 3_2009-02-04 10:43:04.222205 2 POINT (116.36798 39.90471) abc
2009-02-04 10:43:07.000000 3227 775 3_2009-02-04 10:43:04.222205 2 POINT (116.36795 39.9048) abc
2009-02-04 10:43:09.000000 3228 776 3_2009-02-04 10:43:04.222205 2 POINT (116.36793 39.90487) abc
2009-02-04 10:43:11.000000 3229 777 3_2009-02-04 10:43:04.222205 2 POINT (116.36794 39.90495) abc
... ... ... ... ... ... ...
2009-02-04 10:48:14.000000 3390 938 3_2009-02-04 10:43:04.222205 2 POINT (116.36535 39.90589) abc
2009-02-04 10:48:15.000000 3391 939 3_2009-02-04 10:43:04.222205 2 POINT (116.36525 39.90589) abc
2009-02-04 10:48:16.000000 3392 940 3_2009-02-04 10:43:04.222205 2 POINT (116.36516 39.90589) abc
2009-02-04 10:48:17.000000 3393 941 3_2009-02-04 10:43:04.222205 2 POINT (116.36507 39.90589) abc
2009-02-04 10:48:17.399995 3393 941 3_2009-02-04 10:43:04.222205 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_2009-02-04 10:43:04.222205 2 POINT (116.36799 39.90468) abc
2009-02-04 10:43:05.000000 3226 774 3_2009-02-04 10:43:04.222205 2 POINT (116.36798 39.90471) abc
2009-02-04 10:43:07.000000 3227 775 3_2009-02-04 10:43:04.222205 2 POINT (116.36795 39.9048) abc
2009-02-04 10:43:09.000000 3228 776 3_2009-02-04 10:43:04.222205 2 POINT (116.36793 39.90487) abc
2009-02-04 10:43:11.000000 3229 777 3_2009-02-04 10:43:04.222205 2 POINT (116.36794 39.90495) abc
... ... ... ... ... ... ...
2009-03-10 11:08:05.000000 4934 672 4_2009-03-10 11:01:41.826068 2 POINT (116.36539 39.90588) abc
2009-03-10 11:08:06.000000 4935 673 4_2009-03-10 11:01:41.826068 2 POINT (116.36529 39.90588) abc
2009-03-10 11:08:07.000000 4936 674 4_2009-03-10 11:01:41.826068 2 POINT (116.36519 39.90588) abc
2009-03-10 11:08:08.000000 4937 675 4_2009-03-10 11:01:41.826068 2 POINT (116.36509 39.90588) abc
2009-03-10 11:08:08.575754 4937 675 4_2009-03-10 11:01:41.826068 2 POINT (116.36504 39.90588) abc

375 rows × 6 columns

In [ ]: