Fill gaps, smooth, compare and convert trajectories

Question. How do I fill gaps in, smooth, compare, convert and store trajectories, and how do I put them in the right time zone?

You will get a trajectory with its recording gaps filled, a position at a chosen time, a smoothed track, a distance between two tracks, local-time timestamps, H3 grid cells, a grid built with a tiler and a spatial join to it, groups of trajectories with similar shapes, a GeoDataFrame, and a file round trip.

Prerequisites. Start with a trajectory. The trajectory table in this notebook is a TrajDataFrame, the library's typed trajectory wrapper; every function also accepts a plain pandas or Polars dataframe.

Data. GeoLife, the first 5 trips of users 000 and 001 (about 20,500 fixes, about 2 MB of the 313 MB archive), read with fastmob.data.load_geolife. Timestamps are UTC. The grid used for the spatial join is a synthetic square grid built by a tiler. The notebook stops if the download fails.

%pip install -q 'fastmob>=0.3.4' pandas matplotlib geopandas shapely polars

import tempfile
from pathlib import Path

Load and inspect

load_geolife returns a TrajDataFrame. .to_pandas() gives an ordinary dataframe when you need one.

import fastmob
import numpy as np
import pandas as pd
from fastmob.data import load_geolife
from fastmob.trajectory import interpolate, interpolate_at, smooth, trajectory_distance

tdf = load_geolife(users=["000", "001"], max_trips_per_user=5)

fixes = tdf.to_pandas()
print(len(fixes), "fixes;", fixes["uid"].nunique(), "users")
fixes.head(3)

20500 fixes; 2 users

Put timestamps in local time

GeoLife is recorded in UTC but collected in Beijing (UTC+8). Anything that depends on the hour of the day, such as home and work or daily patterns, must use local time. timezone_conversion converts the datetime column in place from one time zone to another; the result has no time-zone information attached (it is tz-naive).

local = fastmob.TrajDataFrame(fixes.copy())
local.timezone_conversion(from_timezone="UTC", to_timezone="Asia/Shanghai")

pd.concat(
    [
        fixes["datetime"].head(2).rename("utc"),
        local.to_pandas()["datetime"].head(2).rename("beijing"),
    ],
    axis=1,
)

The Beijing column is eight hours ahead of UTC: a fix at 02:53 UTC happened at 10:53 local time. The conversion changes the table in place and drops the time-zone label, so do not run it twice on the same table.

Keep the data in time order

Most measures sort internally, but sort_by_uid_and_datetime returns a table ordered by user and then time, which is what you want before writing or exporting a table.

ordered = tdf.sort_by_uid_and_datetime()
type(ordered).__name__, len(ordered)

('TrajDataFrame', 20500)

The table has the same 20,500 rows, now ordered by user and time.

Fill recording gaps

GPS recording is irregular. interpolate looks for gaps between consecutive fixes that are longer than sampling_rate_s seconds (here 60) and inserts one new position in each, by linear interpolation between the two neighbouring fixes. The original fixes are all kept, so the result is never thinner than the input and is not a regular time grid. The same operation is available as the function interpolate and as the method TrajDataFrame.interpolate. To thin a trajectory use compress or simplify (see the previous notebook).

filled = interpolate(fixes, method="linear", sampling_rate_s=60)
filled_m = tdf.interpolate(method="linear", sampling_rate_s=60)

gaps = int((fixes.groupby("uid")["datetime"].diff().dt.total_seconds() > 60).sum())
print(
    f"{len(fixes)} fixes -> {len(filled)} positions: {len(filled) - len(fixes)} inserted, {gaps} gaps longer than 60 s"
)
print("the method gives the same number of rows:", len(filled_m) == len(filled))
filled.head(3)

20500 fixes -> 20589 positions: 89 inserted, 89 gaps longer than 60 s


the method gives the same number of rows: True

89 positions were inserted (20,500 to 20,589), exactly one for each of the 89 gaps longer than 60 seconds. GeoLife records every 1 to 5 seconds, so most fixes are already denser than the threshold and nothing is added there. The first rows of the result are the original fixes. The method and the function give the same number of rows.

Ask for a position at a given time

interpolate_at returns each user's estimated position at the timestamps you list, with a valid flag that is False when the user was not recorded at that time. Pass a list of timestamps (strings or datetime64 values).

at_time = interpolate_at(fixes, at=["2008-10-23 03:00:00"])
at_time_m = tdf.interpolate_at(["2008-10-23 03:00:00"])

print(
    "the method gives the same positions:",
    at_time_m.to_pandas()["lat"].fillna(0).round(6).tolist()
    == at_time["lat"].fillna(0).round(6).tolist(),
)
at_time

the method gives the same positions: True

User 000 has an estimated position at 03:00 UTC on 23 October 2008 (valid is True). User 001 has no estimate (NaN, valid is False), because that time lies outside what was recorded for this user. Always filter on valid before using the coordinates.

Smooth GPS noise

smooth replaces each position with a Kalman-filter estimate at the same timestamp, so the number of rows does not change. process_noise_std_km sets how much the track may deviate from straight-line motion: larger values smooth less.

smoothed = smooth(fixes, method="kalman_cv")
smoothed_m = tdf.smooth()

print(
    "the method gives the same smoothed latitudes:",
    bool(np.allclose(smoothed_m.to_pandas()["lat"], smoothed["lat"])),
)
shift_m = (
    np.hypot(
        smoothed["lat"] - fixes["lat"],
        (smoothed["lng"] - fixes["lng"]) * np.cos(np.radians(fixes["lat"])),
    )
    * 111_320
)
print(
    f"median shift {shift_m.median():.1f} m, 95th percentile {shift_m.quantile(0.95):.1f} m"
)

the method gives the same smoothed latitudes: True
median shift 0.3 m, 95th percentile 1.5 m

The filter moves the median point by about 0.3 m and 95% of points by less than 1.5 m. That is light denoising: consumer GPS error is several metres, so the result stays close to the measurements. Increase the noise setting to smooth less, decrease it to smooth more, and check the effect on the measure you care about.

Compare two trajectories

trajectory_distance measures how different two tracks are. The default dtw (dynamic time warping) lets the two tracks run at different speeds and returns a distance in kilometres summed over the best alignment. It compares exactly two single tracks, so slice one user each. We use the first 300 fixes of each.

first, second = (
    fixes[fixes["uid"] == "000"].head(300),
    fixes[fixes["uid"] == "001"].head(300),
)
distance = trajectory_distance(first, second, method="dtw")
distance_m = fastmob.TrajDataFrame(first).trajectory_distance(
    fastmob.TrajDataFrame(second), method="dtw"
)

print(
    f"DTW distance between user 000 and user 001: {distance:.0f} km (method: {distance_m:.0f} km)"
)

DTW distance between user 000 and user 001: 496 km (method: 496 km)

The distance is about 496 km summed over the best alignment of 300 fixes against 300 fixes, an average of about 1.7 km between aligned points: the two users were not moving along the same route. The value grows with the number of points, so compare tracks of similar length.

Summaries from the same object

Several measures are also methods of the TrajDataFrame, so a quick summary needs no imports: jump_lengths, radius_of_gyration, compress, stay_locations and work_location. They are the same functions used in the earlier notebooks.

radius = tdf.radius_of_gyration()
jumps = tdf.jump_lengths()
compressed = tdf.compress(spatial_radius_km=0.2)
stops = tdf.stay_locations(minutes_for_a_stop=20, spatial_radius_km=0.2)
work = tdf.work_location()

print(len(compressed), "fixes after compress;", len(stops), "stops")
print("work location of each user:", work.to_pandas().round(4).values.tolist())
print(
    "mean jump length per user (km):",
    [round(sum(j) / len(j), 3) for j in jumps.to_pandas()["jump_lengths"]],
)
radius.to_pandas()

629 fixes after compress; 32 stops
work location of each user: [['000', 39.9946, 116.3268], ['001', 40.0138, 116.3065]]
mean jump length per user (km): [np.float64(0.02), np.float64(0.009)]

The methods return the same results as the functions: compression leaves 629 of 20,500 fixes, 32 stops are found, and the radius of gyration is about 5.5 km for user 000 and 6.3 km for user 001 (in kilometres). Measures that return one row per user come back as an Arrow table; .to_pandas() converts it.

Move between tables

settings_from copies the column settings of one TrajDataFrame to another (useful after building a new table from an old one). to_geodataframe converts to a GeoPandas table with point geometry, which needs geopandas; it is installed by the first cell of this notebook.

clone = fastmob.TrajDataFrame(fixes.copy())
clone.settings_from(tdf)
gdf = tdf.to_geodataframe()

print(type(gdf).__name__, gdf.crs, len(gdf))
gdf.head(2)

GeoDataFrame EPSG:4326 20500

The GeoDataFrame has one point per fix in longitude and latitude (EPSG:4326), so it can be used with any GeoPandas tool.

Build a grid and snap positions to it

A tiler covers an area with cells and returns a GeoPandas table with a tile_id column and polygon geometry. They are created through a registry: tilers.tiler.create(name, ...) (and get, which behaves the same) builds a tessellation with the tiler registered under name (squared, h3_tessellation or voronoi), and register_tiler adds your own under a new name. We cover the extent of the data with squares of 2 km. mapping then assigns each fix to the tile that contains it.

import geopandas as gpd
from fastmob.tessellation import tilers
from shapely.geometry import box

extent = gpd.GeoDataFrame(
    geometry=[
        box(
            fixes["lng"].min(),
            fixes["lat"].min(),
            fixes["lng"].max(),
            fixes["lat"].max(),
        )
    ],
    crs="EPSG:4326",
)
tessellation = tilers.tiler.create("squared", base_shape=extent, meters=2000)
same_grid = tilers.tiler.get("squared", base_shape=extent, meters=2000)
tilers.tiler.register_tiler("squared_copy", tilers.SquaredTessellationTiler())
mapped = tdf.mapping(tessellation, remove_na=True)

print(
    len(tessellation),
    "square tiles;",
    len(same_grid),
    "from get;",
    sorted(tilers.tiler._tilers),
)
mapped.to_pandas().head(3)

196 square tiles; 196 from get; ['h3_tessellation', 'squared', 'squared_copy', 'voronoi']

The 2 km tiler covers the extent of the data with 196 square tiles, and get builds the same grid. After register_tiler the registry lists a fourth tiler, squared_copy, next to h3_tessellation, squared and voronoi. Each fix now has the tile_id of the tile that contains it (the first fixes fall in tile 133). The grid is synthetic: its ids mean nothing outside this example.

Group trajectories by shape

cluster_trajectory_shapes groups trajectories by the shape of their path, regardless of where they are, how large they are or how they are oriented (a signature of distances between points of each path, clustered with DBSCAN). It takes a list of trajectories, each a table of one trip, and returns one cluster number per trajectory (-1 means no cluster). cluster_trajectory_shapes_from_segments does the same for a table with a segment id column, such as the output of segment. A segment needs at least two points (we keep the ones with at least 100). epsilon is the neighbourhood size in the normalised signature space: larger values merge more.

from fastmob.trajectory import (
    cluster_trajectory_shapes,
    cluster_trajectory_shapes_from_segments,
)

segmented = fastmob.segment(fixes, method="observation_gap", gap_s=30 * 60)
sizes = segmented.groupby(["uid", "segment_id"]).size()
long_segments = (
    segmented.set_index(["uid", "segment_id"])
    .loc[sizes[sizes >= 100].index]
    .reset_index()
)
by_segment = cluster_trajectory_shapes_from_segments(
    long_segments, segment_col="segment_id", uid_col="uid", epsilon=0.2
)
trips = [trip for _, trip in long_segments.groupby(["uid", "segment_id"])]
labels = cluster_trajectory_shapes(trips, epsilon=0.2)

print(
    len(trips),
    "segments of at least 100 fixes; cluster sizes:",
    by_segment["shape_cluster"].value_counts().to_dict(),
)
print(
    "same result as the list version:",
    by_segment["shape_cluster"].tolist() == labels.tolist(),
)
by_segment.head()

22 segments of at least 100 fixes; cluster sizes: {0: 11, -1: 4, 1: 3, 2: 2, 3: 2}
same result as the list version: True

Of 22 segments with at least 100 fixes, 11 share one shape cluster (0), four are labelled -1 (no cluster) and the rest form three small clusters, at epsilon = 0.2. The list version and the segment version give the same labels. Shape clustering says which trips look alike (for example straight against looping), not where they are.

Convert and compare typed tables

Every typed table (TrajDataFrame, Staypoints, ...) derives from BaseDataFrame, which provides to_pandas, to_polars and to_native (the library's own storage, here an Arrow table) and compare_to, which applies an evaluation metric to a column of two tables (the default metric is the Wasserstein distance). The same methods also build the interactive charts of the optional fastmob[vis] extra (plot_bar, plot_scatter, ...), described in the visualization reference.

as_polars = tdf.to_polars()
native = tdf.to_native()
comparison = tdf.compare_to(fastmob.TrajDataFrame(smoothed), "lat")

print(
    type(as_polars).__module__,
    "|",
    type(native).__module__,
    "| distance between the latitudes and their smoothed version:",
    round(comparison.overall, 6),
)

polars.dataframe.frame | pyarrow.lib | distance between the latitudes and their smoothed version: 2e-06

to_polars returns a Polars dataframe and to_native the Arrow table that stores the data. The Wasserstein distance between the latitudes and their smoothed version is about 0.000002 degrees (a few tenths of a metre), the same light smoothing seen above.

H3 cells

H3 is a global hexagonal grid. latlng_to_h3 adds a cell id to every row at a chosen resolution (9 is a cell roughly 200 m across) and h3_to_latlng converts cell ids back to the cell centres. Use it to count visits per place without defining zones.

cells_df = fastmob.latlng_to_h3(fixes, resolution=9)
centres = fastmob.h3_to_latlng(cells_df.drop_duplicates("h3_cell")[["h3_cell"]])

print(cells_df["h3_cell"].nunique(), "distinct cells at resolution 9")
centres.head(3)

340 distinct cells at resolution 9

The 20,500 fixes fall in 340 distinct cells at resolution 9. h3_to_latlng gives each cell's centre, which is a convenient label for a place; the centre is not a visited position.

Store and reload

io.write saves a table to CSV (or another supported format, chosen by the file extension) and io.read loads it back. io.load_geolife_trajectories reads GeoLife .plt files from a folder, which is what load_geolife calls after downloading.

from fastmob import io

path = Path(tempfile.mkdtemp()) / "geolife_subset.csv"
io.write(tdf, path)
reloaded = io.read(path)
from_plt = io.load_geolife_trajectories(
    Path(tempfile.gettempdir()) / "fastmob_data" / "geolife", user_ids=["000"]
)

print(
    len(reloaded),
    "rows reloaded;",
    len(from_plt),
    "rows read from the cached .plt files of user 000",
)

20500 rows reloaded; 3424 rows read from the cached .plt files of user 000

The round trip returns all 20,500 rows. Reading the cached .plt files of user 000 directly returns 3,424 rows, the five trips that load_geolife extracted for that user.

Which datasets does the library list?

list_datasets and get_dataset_info describe the datasets of the scikit-mobility collection (name, type, source, size).

from fastmob.data import get_dataset_info, list_datasets

names = list_datasets()
print(names)
{
    k: get_dataset_info("foursquare_nyc")[k]
    for k in ("name", "data_type", "download_format")
}

['flow_foursquare_nyc', 'foursquare_nyc', 'nyc_boundaries', 'parking_san_francisco', 'taxi_san_francisco']





{'name': 'Foursquare_NYC', 'data_type': 'trajectory', 'download_format': 'zip'}

Five datasets are listed. foursquare_nyc is a trajectory dataset distributed as a zip file. Only the metadata is read here; the data themselves are not downloaded.

Interpretation and limits

• Gap filling creates positions that were never observed, and only one per gap. Interpolated points are suitable for plots, not as evidence of where someone was, and the result is not a regular time series: resample explicitly (for example with interpolate_at on a time grid you choose) if you need one.
• smooth replaces measured positions with model estimates; the filter assumes roughly constant velocity and handles sharp turns and stops poorly.
• Dynamic time warping compares shapes; it ignores clock time and the scale is in kilometres summed along the alignment, so it grows with track length.
• mapping gives a tile id only to fixes inside a polygon; with remove_na=True fixes outside the grid are dropped. Tile ids and sizes depend on the tiler and meters.
• Shape clustering depends on epsilon: with a small value most trajectories are labelled -1 (noise), with a large one they all merge.
• load_dataset (the scikit-mobility download) is not demonstrated: its source host did not resolve when this notebook was prepared, so there is no executed example of it.

Next steps

• Build stops, trips and tours from these fixes: Build the trajectory hierarchy.
• Count flows between zones: Flows between places.

Related API

interpolate · interpolate_at · smooth · trajectory_distance · cluster_trajectory_shapes · cluster_trajectory_shapes_from_segments · TessellationTilers · SquaredTessellationTiler · H3TessellationTiler · VoronoiTessellationTiler · to_pandas · to_polars · to_native · compare_to · timezone_conversion · sort_by_uid_and_datetime · interpolate · interpolate_at · smooth · trajectory_distance · jump_lengths · radius_of_gyration · compress · stay_locations · work_location · settings_from · to_geodataframe · mapping · latlng_to_h3 · h3_to_latlng · read · write · load_geolife_trajectories · list_datasets · get_dataset_info