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 

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

SYSTEM INFO
-----------
python     : 3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:34:57) [MSC v.1936 64 bit (AMD64)]
executable : H:\miniconda3\envs\mpd-ex\python.exe
machine    : Windows-10-10.0.19045-SP0

GEOS, GDAL, PROJ INFO
---------------------
GEOS       : None
GEOS lib   : None
GDAL       : 3.7.0
GDAL data dir: None
PROJ       : 9.2.1
PROJ data dir: H:\miniconda3\pkgs\proj-9.0.0-h1cfcee9_1\Library\share\proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 0.13.2
pandas     : 2.0.3
fiona      : 1.9.4
numpy      : 1.24.4
shapely    : 2.0.1
rtree      : 1.0.1
pyproj     : 3.6.0
matplotlib : 3.7.2
mapclassify: 2.5.0
geopy      : 2.3.0
holoviews  : 1.17.0
hvplot     : 0.8.3
geoviews   : 1.9.6
stonesoup  : 1.0
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

Clipping a TrajectoryCollection¶

Alternatively, using TrajectoryCollection:

In [5]:
clipped = tc.clip(polygon)
clipped
Out[5]:
TrajectoryCollection with 2 trajectories
In [6]:
clipped.plot()
Out[6]:
<Axes: >
No description has been provided for this image

Computing intersections¶

In [7]:
polygon_feature = {
    "geometry": polygon,
    "properties": {'field1': 'abc'}
}
In [8]:
my_traj = tc.trajectories[2]
intersections = my_traj.intersection(polygon_feature)
intersections
Out[8]:
TrajectoryCollection with 1 trajectories
In [9]:
intersections.to_point_gdf()
Out[9]:
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.90467) 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.90480) 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

In [ ]: