GTFS data segmentation using R

[subha-nair]

Hi, is there a way to cut the GTFS shapes file into segments where they overlap, and visualise the overlapped lines as a single line with frequency of bus routes on them (shown as line thickness). I am new to R so a step by step process would be helpful.

[polettif]

There is a way but unfortunately there is no easy step-by-step workflow to do this. The main issue is identifying overlaps. If your feed has precise shapes.txt coordinates you can use sf::st_segmentize() to create segments. You then compare all segments. I'd guess you'd like to summarise shapes that are "close-by but not equal", then you'd have to use st_distance(..., which="Hausdorff") or something like that.

So it's mainly a not-so-trivial geospatial issue. As you're new to R I'd suggest looking at Geocomputation with R, I don't think I can help you with that further.

Summarising the number of trips (not routes, as a route might have more than one trip with a corresponding shape) on a shape can be done like this.

library(tidytransit)
library(dplyr)

g = gtfs_duke |>
    gtfs_as_sf()

trip_counts = g$trips |>
    distinct(trip_id, shape_id) |>
    count(shape_id, name = "n_trips")

shapes = left_join(g$shapes, trip_counts, "shape_id")

library(ggplot2)
ggplot(shapes) +
    geom_sf(aes(linewidth = n_trips))

^{Created on 2024-12-17 with reprex v2.1.1}

Note that in this plot there are several overlapping shapes so the linewidth doesn't really show the right number of trips. So this is just a starting point and doesn't solve the issue. To fix overlaps you'd have to solve the geospatial issue.

[idshklein]

if you could give an example workflow to the hausdorf method mentioned it will be appriciated

[polettif]

As I said, it's not a trivial example to do. If you want to do matching you need to link segments with shapes (which then can be linked to trips and routes). You can create segments like this:

library(dplyr)
library(tidytransit)

g = gtfs_as_sf(gtfs_duke, crs = 32119) # EPSG code for the feed area

shape_segments = g$shapes |>
	sf::st_segmentize(50) |> # create segments with at most 50 m distance 
	nngeo::st_segments() |> # split shape linestrings into segmented linestrings
	mutate(segment_id = paste0(shape_id, "_", row_number()), .by = shape_id)

After that it essentially becomes a map matching problem and thus heavily depends on how your data is structured. Do you want to match each direction independently? What's the maximum distance you consider two shape segments as "equal"? Do you want to match routes to existing path segments (if so, which segment should be used?), to newly created segments or another data source?

I'm sorry if this is not really helpful but in my experience gtfs feeds rarely provide clean enough data to use without quite some experimenting and fiddling. I did try it with the duke dataset for a bit but didn't get it to work. Also, my initial guess with Hausdorff distance probably doesn't work since two consecutive segments have a distance of 0 but you probably don't want to match those.

[polettif]

In the spirit of advent of code I couldn't help but tinker with this a bit more 🎄🎅 🎁

The following code matches nearby shape points and finds overlapping shapes along segments. It looks ok but I didn't inspect the plausibility of the example dataset further and there might be pitfalls with other dataset. So it's not a catch-all workflow but it should provide a starting point.

Snap nearby shape coordinates to allow matching

library(dplyr)
library(sf)
library(tidytransit)

# Add interpolated points along shape lines
add_shape_points = function(gtfs_shapes, crs, max_segment_dist) {
    gtfs_shapes |>
        shapes_as_sf(crs) |>
        st_segmentize(max_segment_dist) |>
        tidytransit:::sf_lines_to_df() |> # direct access fn in tidytransit::sf_as_tbl()
        as_tibble()
}

# Build hierarchical clusters with all gtfs-shape points and combine them
# within max_dist. The returned `shapes_snapped` has the same number of shape points
# but their lonlat coordinates have been change to their cluster's center.
snap_simplify_shapes = function(gtfs_shapes, crs, max_cluster_dist) {
    stopifnot(inherits(gtfs_shapes, "data.frame"))
    stopifnot(!inherits(gtfs_shapes, "sf"))

    # We don't want to convert shapes to actual lines (with gtfs_sf_sf() or shapes_as_sf()),
    # the point coordinates in shapes.txt are actually already fine for our use case.
    # However, we need to convert them from lon/lat to a projected coordinate system,
    # for the duke example feed it's EPSG:32119
    shapes = gtfs_shapes |>
        st_as_sf(coords = c("shape_pt_lon", "shape_pt_lat"), crs = 4326, remove = FALSE) |>
        st_transform(crs)

    # add point coordinates as columns, round coordinates
    shapes <- shapes |>
        bind_cols(round(st_coordinates(shapes), 1)) |>
        st_drop_geometry() |>
        mutate(shape_point_id = paste0(X, "-", Y))

    # use only unique point coordinates, reduce clustering runtime
    coords = shapes |>
        distinct(shape_point_id, X, Y)

    # hierarchical clustering
    htree = cluster::agnes(coords[,c("X", "Y")])

    cluster_coords = function(htree, coords, k) {
        coords$cluster <- cutree(htree, k = k)

        # Use mean coordinates as cluster center
        # Maybe this data is already available somewhere with the cluster package
        coords_clustered = coords |>
            group_by(cluster) |>
            mutate(cluster_X = mean(X), cluster_Y = mean(Y)) |>
            ungroup()

        # calculate euclidian distance to the cluster center
        cluster_stats = coords_clustered |>
            mutate(dist_to_mean = sqrt((X-cluster_X)^2+(Y-cluster_Y)^2)) |>
            summarise(min_dist = min(dist_to_mean), max_dist = max(dist_to_mean)) |>
            ungroup()

        attributes(coords_clustered)$cluster_stats <- cluster_stats
        attributes(coords_clustered)$max_dist <- max(cluster_stats$max_dist)
        coords_clustered
    }

    # find k value where distance to cluster center is less than max_dist (for all clusters)
    # first increase k in bigger intervals, afterwards finetune k back in 1-step-decrements
    for(k in seq(1, nrow(coords), 100)) {
        coords_clustered = cluster_coords(htree, coords, k)
        if(attributes(coords_clustered)$max_dist < max_cluster_dist) break
    }
    for(k in seq(k, 1, -1)) {
        coords_clustered = cluster_coords(htree, coords, k)
        if(attributes(coords_clustered)$max_dist > max_cluster_dist) break
    }

    # join shapes with new clustered coordinates and transform back to WGS 84
    shapes_snapped = shapes |>
        inner_join(coords_clustered, c("shape_point_id", "X", "Y")) |>
        select(-shape_pt_lat, -shape_pt_lon, -X, -Y) |>
        st_as_sf(coords = c("cluster_X", "cluster_Y"), crs = crs) |>
        st_transform(4326)

    # use new cluster lon/lat coordinates and trim columns to original columns
    shapes_snapped <- shapes_snapped |>
        bind_cols(round(st_coordinates(shapes_snapped), 6)) |> # up to 11 cm accuracy
        st_drop_geometry() |>
        rename(shape_pt_lon = X, shape_pt_lat = Y)

    shapes_snapped <- shapes_snapped[,colnames(gtfs_shapes)]

    return(shapes_snapped)
}

With the functions out of the way, let's simplify the shapes for the example gtfs_duke dataset:

# Option 1: snap shapes with existing shape point coordinates
# new_shapes = snap_simplify_shapes(gtfs_duke$shapes, 32119, 15)

# Option 2: add more points along shape line to better match lines
# between original shape points
new_shapes = gtfs_duke$shapes |>
    add_shape_points(32119, 20) |>
    snap_simplify_shapes(32119, 30)

Compare original shapes to snapped shapes:

mapview::mapview(shapes_as_sf(gtfs_duke$shapes), color = "blue", legend = F)@map |>
    leaflet::setView(-78.934, 36.005, 17)

mapview::mapview(shapes_as_sf(new_shapes), color = "red", legend = F)@map |>
    leaflet::setView(-78.934, 36.005, 17)

Find overlapping segments of shapes

# Create "segments" between individual shape points and create a table linking
# shape_ids and newly created segment_ids. Each shape has 1:n segments (probably
# way more than 1) and each segment has 1:m shapes. Segments are one-directional.
# You can build segments with the original shapes table (i.e. without snapping the
# shape coordinates in the previous step if data is already good enough
get_shape_segments = function(gtfs_shapes) {
    gtfs_shapes_offset = gtfs_shapes |>
        mutate(shape_pt_sequence = shape_pt_sequence-1) |>
        filter(shape_pt_sequence > 0) |>
        select(shape_id, shape_pt_lon, shape_pt_lat, shape_pt_sequence)

    # use linestring wkt as id
    # Note that there are overlapping segments with different directions (A->B and B->A) 
    # which are handled separately in the following workflow. To change this we'd need to
    # create the same id for both directions (for example by sorting the coordinates before hashing)
    create_segment_wkt = function(x1,y1,x2,y2) {
        sprintf("LINESTRING (%.6f %.6f, %.6f %.6f)", x1, y1, x2, y2)
    }

    gtfs_shape_segments = inner_join(gtfs_shapes, gtfs_shapes_offset, c("shape_id", "shape_pt_sequence")) |>
        mutate(segment_id = create_segment_wkt(shape_pt_lon.x, shape_pt_lat.x, shape_pt_lon.y, shape_pt_lat.y)) |>
        distinct(shape_id, segment_id)

    return(gtfs_shape_segments)
}

(shape_segments = get_shape_segments(new_shapes))
#> # A tibble: 8,355 × 2
#>    shape_id segment_id                                             
#>    <chr>    <chr>                                                  
#>  1 p_111287 LINESTRING (-78.914605 36.006000, -78.914605 36.006000)
#>  2 p_111287 LINESTRING (-78.914605 36.006000, -78.914368 36.005660)
#>  3 p_111287 LINESTRING (-78.914368 36.005660, -78.914368 36.005660)
#>  4 p_111287 LINESTRING (-78.914368 36.005660, -78.914435 36.005321)
#>  5 p_111287 LINESTRING (-78.914435 36.005321, -78.914435 36.005321)
#>  6 p_111287 LINESTRING (-78.914435 36.005321, -78.914761 36.005071)
#>  7 p_111287 LINESTRING (-78.914761 36.005071, -78.914761 36.005071)
#>  8 p_111287 LINESTRING (-78.914761 36.005071, -78.915102 36.005272)
#>  9 p_111287 LINESTRING (-78.915102 36.005272, -78.915102 36.005272)
#> 10 p_111287 LINESTRING (-78.915102 36.005272, -78.915571 36.005285)
#> # ℹ 8,345 more rows

segments = shape_segments |>
    distinct(segment_id) |>
    mutate(segment_wkt = segment_id) |>
    st_as_sf(wkt = "segment_wkt", crs = 4326)

Analysis workflow, find number of trips per segment

# trips per day
trips_per_day = gtfs_duke$trips |>
    count(service_id, name = "n_trips_on_service") |>
    inner_join(gtfs_duke$.$dates_services, "service_id") |>
    summarise(n_trips = sum(n_trips_on_service), .by = date)

# Example with how many trips run on a given date
n_trips_per_shape = gtfs_duke$.$dates_services |>
    filter(date == "2019-09-05") |>
    inner_join(gtfs_duke$trips, "service_id") |>
    count(shape_id, name = "n_trips")

n_trips_per_segment = shape_segments |>
    inner_join(n_trips_per_shape, "shape_id") |>
    summarise(n_trips = sum(n_trips), .by = segment_id) |>
    inner_join(segments, "segment_id") |>
    st_as_sf() |>
    arrange(n_trips)

# still a bit misleading because overlapping directional segments
library(ggplot2)
ggplot(n_trips_per_segment) +
    theme_bw() +
    geom_sf(aes(linewidth = n_trips, color = n_trips), lineend = "round")

^{Created on 2024-12-19 with reprex v2.1.1}

[mpadge]

@polettif Thanks so much for your detailed consideration here. It's a problem I've had on my to-do list for a while. Almost all shapes contained in GTFS files exist in Open Street Map, so it should be entirely possible to write a function to do this kind of extraction and analysis. The result would then be geographically precise, and not involve any of the inherit inaccuracies that arise from snap-to-grid kind of approaches. The idea would be to input a GTFS feed, with or without a "shapes" file, and get a network representation back with nodes as stops, and full geographic shapes in between, I'd suggest by default in sfnetworks format.

[polettif]

Thanks for your reply @mpadge it is indeed an interesting topic to tackle.

Almost all shapes contained in GTFS files exist in Open Street Map

I'm not quite sure I understand what you mean here. Do you mean the route taken by a vehicle is generally available as roads or rails in OSM? Then I guess this also relates to the discussions in #156 about pfaedle. I looked into map-matching public transit routes for MATSim back at uni but it's been a while and I'm not up to date on research or applications. Or do you mean the routes taken exist as relationships for transit lines in OSM?

The result would then be geographically precise, and not involve any of the inherit inaccuracies that arise from snap-to-grid kind of approaches.

Yes, this would be very desirable and the second part of the described workflow ("Find overlapping segments of shapes") can then be applied to the precise shapes.

I'd just note that no abstraction might not be a desirable outcome for the issue that started this discussion, especially when it comes to visualisation. I think there will always be some level of abstraction/generalisation necessary, starting with bus lanes mapped separately or trains using different platforms. Also, the paths generated by an algorithm might not match the actual route taken by a vehicle and the level of detail of OSM and GTFS (stop locations) can vary quite a bit depending on region. But that's more of a general caveat.

That's just my 2 cents though, I might be overthinking this.

[mpadge]

responses to your questions @polettif:

I had forgotten about pfaedle, which should be considered a reference work in this context. That said, I would also advocate for direct usage of OSM relation data wherever they exist. In those cases, their usage is far simpler and more direct than anything else.

do you mean the routes taken exist as relationships for transit lines in OSM?

Yes, that's exactly what I mean. Examples are https://www.openstreetmap.org/relation/6808203 or https://www.openstreetmap.org/relation/1153295. Those are bus lines, which are generally least likely to exist. Lines for other modes like train, tram, or ferry, may be generally presumed to exist in OSM. Your example of Duke Uni is from Durham NC, which does not seem to have any bus route data (see https://www.openstreetmap.org/node/6509747385, for example). I wouldn't consider your example of trains using different platforms to be especially important - train data are very generally accurately mapped in OSM, including in NC. Mis-mapped platforms are very easily corrected in OSM, and so are unlikely to be that common. Missing bus routes are a far bigger issue, I'd suggest.

I agree with your points about abstraction, but note as always that the best solution for any mismatches between GTFS and OSM data is to update the latter. In a perfect world, both would always agree, and anybody updating OSM to improve matches with actual routes from GTFS publishers is helping to achieve that. But yeah, that kind of general preaching doesn't actually help practical applications...

If we were to proceed and develop some kind of generic algorithm for extracting "shapes" files from OSM data, then I'd suggest an initial routine should by default fail any time that route segment data are missing, and only succeed where all data exist. We could then refine in further iterations to permit partial construction of "shapes", for example by omitting bus routes where these don't exist.

[polettif]

I'm not sure how to proceed here @mpadge as I'm quite overwhelmed by all the possible data quality issues with GTFS and OSM we might run into. For your examples I don't know how to match GTFS routes with the OSM routes with some level of certainty.

I agree with your points about abstraction, but note as always that the best solution for any mismatches between GTFS and OSM data is to update the latter. In a perfect world, both would always agree, and anybody updating OSM to improve matches with actual routes from GTFS publishers is helping to achieve that.

Just on a conceptual level, OSM route relations are related to GTFS routes, not trips. Shapes on the other hand belong to trips and a GTFS route can contain more than one trip/shape. So I'm not sure this synchronised updating works even in a perfect world.

If we were to proceed and develop some kind of generic algorithm for extracting "shapes" files from OSM data, then I'd suggest an initial routine should by default fail any time that route segment data are missing, and only succeed where all data exist. We could then refine in further iterations to permit partial construction of "shapes", for example by omitting bus routes where these don't exist.

A step-by-step algorithm from "best quality" to "we can't even match this route to a network" with possible iterations sounds indeed like the way to go. Let me know if I can help with any development.