GTFS - Missing bus stops in the output (geopackage)

[gigioLR]

Hi, while running the software for Porto Alegre, we noticed that only 338 of the 5524 bus stops in the GTFS file were included in the output (geopackage). What could it be?

GTFS - Porto Alegre, Brazil:
arquivo-gtfs.zip

[carlhiggs]

Hi @gigioLR ,

Thanks for sharing; would you also mind sharing your study region configuration file? - the dates or other aspects like study region boundary configured there also play a role in this, so good to check out both the data and how it's interpreted. If you share that too, I'll aim to check this out in coming couple of days. Cheers, Carl

[gigioLR]

No problem @carlhiggs, here it is:

• Study region configuration file (.yml): https://drive.google.com/file/d/1qU7-ZM8sf8ca_8GHm70Ung4JBeTzZmk4/view?usp=drive_link
• Study region boundary file (.gpkg): https://drive.google.com/file/d/1t9NPR1vYjerWNRu8SVuK_XeINUfkOYN4/view?usp=sharing

For the start date, we considered when the GTFS file was first published ("20190930"). We tested a few options for end dates ("20230314", "20230723", "20231225", "20231231"), but the quantity of bus stops changed little (334 to 338).

Thank you, Carl!

[carlhiggs]

Thanks @gigioLR ,

I've just had a go at running this, and once I got it going (see note below**), I got the same result re count of bus stops.

Looking at the GTFS file, I think I notice something that might explain it --- in stop_times.txt, some stops have arrival_time and departure_time records, but most do not.

For example, stop_id 2642 (not in the 388 included stops) does not,

216-1@1#540,,,2462,41

but stop 6127 (included in the processed 388) does

2543-2@1#1820,19:00:00,19:00:00,6127,65

At a very brief glance, that seems to explain where things may be programatically going wrong --- there are 390 unique stop ids in stop_times after dropping records with NA, but all up there are 5351 unique stop ids. (there are also 173 stop ids in stops.txt that aren't in stop_times.txt)

According to the GTFS reference specification, these are conditionally required --- so in principle, not having arrival and departure times shouldn't be a reason for exclusion. But it seems at some point in the logic effectively, that's what's happening.

So, there are a few considerations here:

• the function get_hlc_stop_frequency assumes departure times are present
• they could be interpolated (that's a possible external solution; and possibly it could be implented - but it would involve assumptions)
• regardless of interpolation, the function also aims to provide location of 'any pt' --- but currently, it excludes those stops without explicit departure and arrival, which as you've found is a big problem
  • so perhaps, before returning the 'stop_frequent' data frame, the stops with service frequency estimates could be left joined onto all stops within the analysed region, leaving some with nulls
  • the user could be alerted that X number of stop_ids were returned with null frequency estimates, and that interpolation of departure time may be required to accurately account for service frequency

I haven't made any change to GTFS code yet - but at a brief look, that is what appears to be the issue, and

• short term the solution would be to create a modified GTFS zip file with an edited stop_times.txt containing interpolated departure time. It not ideal, but it should work, and you would have control and oversight of the interpolation algorithm used to generate departure times.
• medium term, I think we should warn users when this appears to be an issue, and ensure that the total number of stops reflects the total number of stops or the total number of stops with stop_ids that are present in stop_times, regardless of departure time being present (that's a methodological choice that should be thought through). I wouldn't want to rush doing this, as while it might be simple to do, if rushed something could go wrong.
• longer term, we could look to implement an optional function that could be set in configuration to interpolate departure times
• or we find a supporting library that could otherwise calculate headway frequencies better than our in-house function, robust to absence of arrival and/or departure times.

Ours isn't the only GTFS library with this issue it seems -- gtfs-lite the library we currently use for parsing GTFS mentions the issue of interpolation for its trips_at_stops function, and that it won't work if arrival and/or departure aren't provided. There are other solutions around which we might be able to get inspiration or implementation from -- e.g. this feed parser (hasn't been developed in a while) had a github issue related to the same topic, and a fix was implemented - itinero/GTFS#20. A quick search suggested there are a few tools around that have done it, but they aren't necessarily current. If there is a general open algorithm that is straightforward we might be able to add it as an optional setting in configuration, and metadata could record that this has been done.

Hope this helps --- I need to focus on writing for remainder of week, so may not get to making changes related to this for a little while, but I'll let my colleagues know (@shiqin-liu @gboeing ) in case they have any additional thoughts on how to deal with this,

Thanks again for sharing @gigioLR; this is an important issue that needs to be thought through ,

Carl

---

**re the config I had an issue with the configured projection of EPSG 4674 -- on my end this seemed to be in interpreted in degrees not meters, which cause problems when buffering the study region (made a very large area!) which then resulted in a 'segmentation error' when trying to retrieve the OSM streetnetwork for this via osmnx ---- that was the error I was presented with, and I traced it back to the projection and over-large buffer.

Anyway, just wanted to let you know I used this projection instead - which seeemed appropriate for Porto Alegre -- but you would know better than me

## Metadados de Sistema de Referência de Coordenadas 
crs:
    ## nome do 	Sistema de Referência de Coordenadas, SIRGAS2000
    name: SIRGAS 2000 / UTM zone 22S
    ## acrônimo do catálogo padrão que define este SRC, e.g. EPSG
    standard: EPSG
    ## identificador de referência espacial que identifique este SRC de acordo com o padrão específico, e.g. 5635 (veja https://spatialreference.org/)
    srid: 31982

And for population, I noticed you're using the 2022 release of GHS-POP; just wanted to let you know they've released an updated version of this, which if you want you can download from URL in config below (i haven't yet implemented the suggestion of attempting to auto-download using configuration - its a good idea, and will try to get to it sometime):

## Metadados da população (raster ou vetor)
population:
    ## nome dos dados populacionais
    name: "Dados populacionais dos Assentamentos Humanos Globais: 2020, Mollweide (EU JRC, 2022)"
    ## caminho relativo ao diretório de dados do projeto para a pasta contining tifs ou para o arquivo vetorial
    data_dir: population_grids/GHS_POP_E2020_GLOBE_R2023A_54009_100_V1_0_R13_C14
    ## tipo de dado (por exemplo, "raster:Int64" ou "vector")
    data_type: raster:Float64
    ###########
    ## Campos específicos de dados vetoriais; descomente se estiver usando dados vetoriais (por exemplo, shp, geojson, gpkg)
    ## O campo da coluna com estimativas populacionais para o grupo populacional de interesse (pode ser o total ou para um subgrupo de interesse que seja o foco de seus indicadores)
    # vector_population_data_field:
    ## O campo com os dados de população total dentro dessa área (se você estiver interessado no total, deve ser o mesmo valor que vector_population_data_field)
    # population_denominator:
    ###########
    ###########
    ## Campos específicos de dados raster; comente se estiver usando dados raster (ou seja, tif)
    ## resolução da imagem, por exemplo, 100 m
    resolution: 100m
    ## a banda de imagem contendo os dados relevantes, por exemplo, para GHSL-POP, 1
    raster_band: 1
    ## Um valor na imagem que representa 'nenhum dado', por exemplo, para GHSL-POP, -200
    raster_nodata: -200
    ###########
    ## Pontos de amostra que intersectam com células de grade com população estimada menor que este valor serão excluídos da análise. Dependendo da resolução dos dados de população, você pode usar isso para excluir áreas com população muito baixa devido à incerteza de onde alguém pode morar nessa área, ou se moram de fato
    pop_min_threshold: 1
    ## Metadados do sistema de coordenadas de referência para os dados populacionais (por exemplo, Mollweide, ESRI, 54009)
    crs_name: Mollweide
    crs_standard: ESRI
    crs_srid: 54009
    ## URL de onde esses dados foram obtidos
    source_url: https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_POP_GLOBE_R2023A/GHS_POP_E2020_GLOBE_R2023A_54009_100/V1-0/tiles/GHS_POP_E2020_GLOBE_R2023A_54009_100_V1_0_R14_C13.zip
    ## metadados para citação
    ## quando foi publicado (aaaa), por exemplo, 2023
    year_published: 2022
    ## o ano que se pretende representar (aaaa), por exemplo, 2020
    year_target: 2020
    ## quando você o adquiriu (aaaammdd). Isso pode ser útil para registrar, pois os dados podem ser objeto de revisão. por exemplo, 20230627
    date_acquired: 20230222
    ## licença, por exemplo, "CC BY 4.0"
    licence: CC BY 4.0
    ## citação, por exemplo, "Schiavina, M; Freire, S; Carioli, A., MacManus, K (2023): GHS-POP R2023A - GHS population grid multitemporal (1975-2030). European Commission, Joint Research Centre (JRC) [Dataset] doi: 10.2905/D6D86A90-4351-4508-99C1-CB074B022C4A"
    citation: "Schiavina, Marcello; Freire, Sergio; MacManus, Kytt (2022): GHS-POP R2022A - GHS population grid multitemporal (1975-2030). European Commission, Joint Research Centre (JRC) [Dataset] doi: 10.2905/D6D86A90-4351-4508-99C1-CB074B022C4A"

[gigioLR]

Hi @carlhiggs, thanks for the feedback! We imagined that something was missing in the GTFS file...

About the CRS, yes, we actually tested both of them a few days ago, so your suggestion is correct. Thanks for pointing out about using a projected CRS (EPSG 31982) instead of a geographic one (EPSG 4674) in the config file!

On the GHL-POP 2023 update, thanks for reminding this, we used the updated version for other cities but not for Porto Alegre.

Giovani

[shiqin-liu]

Blank stop times have always been an issue for generating the transit network. The only tool I am currently award of is ArcPy form ArcGIS https://pro.arcgis.com/en/pro-app/latest/tool-reference/public-transit/gtfs-to-public-transit-data-model.htm. They have an option for interpolating the empty stop times.

I have not used the tool before, what we usually approach this issue is to manually change the GTFS stop_times with estimating departure and arrival time (might require some local knowledge), as Carl you also suggest.

I agree with adding this as a warning message to alert the users.

[gigioLR]

Thanks for the insight @shiqin-liu!

[carlhiggs]

I'll convert this to an issue now

[carlhiggs]

Thanks again for raising this issue @gigioLR and for sharing your experiences with this @shiqin-liu,

As per above discussion on #359 I've now ticked boxes for:

• Add a note to example and template region configuration files to advise that GTFS feeds missing departure times will require interpolation for valid results
• Add a warning and/or raise exception if departure times are missing
• Look into using stops.txt for stop locations, rather than stop_times (left join before returning), to ensure the list of 'any' stops derived from GTFS contians all, regardless of information on service frequency

Re the last point I'm curious to hear your thoughts and from @gboeing on whether you think this is a good idea (commit link to view code difference: global-healthy-liveable-cities@e6a6797)

We already did do a left join with stops, but previously we left-joined the unique 'stop_times' stops with departure IDs and frequency estimates with 'stops' (to associate with geometries). So, I just shifted the order around so the index are 'all stops' regardless of being listed in stop_times. So if a few stops are listed in stops.txt but not stop_times.txt (as is apparently the case with Las Palmas), the denominator will be correspondingly larger and thus shift the frequency estimates In the case of Las Palmas, only very slightly --- but in the case of Porto Alegre, for which most stops lack departure times, the denominator goes up from 338 to 5187, with the result that the GTFS 'any pt' is accurate, but percentage estimates for frequency of service are incorrect and can't be interpreted (previously the interpretation, tacit as it wasn't really intentional, was "the percentage of stops with departure times that have headway <= threshold mins" --- now its "the percentage of stops (regardless of presence of service information) that have a headway <= threshold mins".

Not sure if that explanation is clear but easy to see implification if you look at the output:

Las Palmas, before modification to retain stops without stop_times

Las Palmas de Gran Canaria analysis:
  - /home/ghsci/process/data/transit_feeds/Example/gtfs_es_las_palmas_de_gran_canaria_guaguas_20230222.zip)
  - 20230405 to 20230605)
  - ['07:00:00', '19:00:00']


     Tram                  0 stops identified (0.05 seconds)
     Metro                 0 stops identified (0.07 seconds)
     Rail                  0 stops identified (0.05 seconds)
     Bus                 788 stops identified (0.71 seconds)
     Ferry                 0 stops identified (0.05 seconds)
     Cable tram            0 stops identified (0.05 seconds)
     Aerial lift           0 stops identified (0.06 seconds)
     Funicular             0 stops identified (0.05 seconds)
     Trolleybus            0 stops identified (0.05 seconds)
     Monorail              0 stops identified (0.05 seconds)

Las Palmas de Gran Canaria summary (all feeds):
       tot_stops  headway<=30  headway<=20  pct_headway<=30  pct_headway<=20
mode
Bus          788          675          448            85.66            56.85
total        788          675          448            85.66            56.85

and after (minor differnce in 'tot_stops' denominator, hence small shift in percentage estiimate)

Las Palmas de Gran Canaria analysis:
  - /home/ghsci/process/data/transit_feeds/Example/gtfs_es_las_palmas_de_gran_canaria_guaguas_20230222.zip)
  - 20230405 to 20230605)
  - ['07:00:00', '19:00:00']


     Tram                  0/793 (0.0%) stops identified with departure times (0.07 seconds)
     Metro                 0/793 (0.0%) stops identified with departure times (0.05 seconds)
     Rail                  0/793 (0.0%) stops identified with departure times (0.05 seconds)
     Bus                 788/793 (99.4%) stops identified with departure times (0.71 seconds)
     Ferry                 0/793 (0.0%) stops identified with departure times (0.05 seconds)
     Cable tram            0/793 (0.0%) stops identified with departure times (0.07 seconds)
     Aerial lift           0/793 (0.0%) stops identified with departure times (0.08 seconds)
     Funicular             0/793 (0.0%) stops identified with departure times (0.05 seconds)
     Trolleybus            0/793 (0.0%) stops identified with departure times (0.05 seconds)
     Monorail              0/793 (0.0%) stops identified with departure times (0.08 seconds)

Las Palmas de Gran Canaria summary (all feeds):
       tot_stops  headway<=30  headway<=20  pct_headway<=30  pct_headway<=20
mode
Bus          793          675          448            85.12            56.49
total        793          675          448            85.12            56.49

For Porto Alegre, before switching the left join -- but including the warning:

**WARNING**: 5237 stops with null departure times found in stop_times.txt.  Service frequencies will likely be inaccurate for these locations.

Porto Alegre analysis:
  - /home/ghsci/process/data/transit_feeds/poa/arquivo-gtfs.zip)
  - 20190930 to 20230723)
  - ['07:00:00', '19:00:00']


     Tram                  0 stops identified (0.10 seconds)
     Metro                 0 stops identified (0.09 seconds)
     Rail                  0 stops identified (0.09 seconds)
     Bus                 338 stops identified (1.00 seconds)
     Ferry                 0 stops identified (0.10 seconds)
     Cable tram            0 stops identified (0.10 seconds)
     Aerial lift           0 stops identified (0.10 seconds)
     Funicular             0 stops identified (0.10 seconds)
     Trolleybus            0 stops identified (0.11 seconds)
     Monorail              0 stops identified (0.09 seconds)

Porto Alegre summary (all feeds):
       tot_stops  headway<=30  headway<=20  pct_headway<=30  pct_headway<=20
mode
Bus          338          191          147            56.51            43.49
total        338          191          147            56.51            43.49

Porto Alegre after retaining all stops ---

**WARNING**: 5237 stops with null departure times found in stop_times.txt.  Service frequencies will likely be inaccurate for these locations.

Porto Alegre analysis:
  - /home/ghsci/process/data/transit_feeds/poa/arquivo-gtfs.zip)
  - 20190930 to 20230723)
  - ['07:00:00', '19:00:00']


     Tram                  0/5187 (0.0%) stops identified with departure times (0.10 seconds)
     Metro                 0/5187 (0.0%) stops identified with departure times (0.11 seconds)
     Rail                  0/5187 (0.0%) stops identified with departure times (0.10 seconds)
     Bus                 338/5187 (6.5%) stops identified with departure times (0.98 seconds)
     Ferry                 0/5187 (0.0%) stops identified with departure times (0.11 seconds)
     Cable tram            0/5187 (0.0%) stops identified with departure times (0.09 seconds)
     Aerial lift           0/5187 (0.0%) stops identified with departure times (0.09 seconds)
     Funicular             0/5187 (0.0%) stops identified with departure times (0.09 seconds)
     Trolleybus            0/5187 (0.0%) stops identified with departure times (0.12 seconds)
     Monorail              0/5187 (0.0%) stops identified with departure times (0.09 seconds)

Porto Alegre summary (all feeds):
       tot_stops  headway<=30  headway<=20  pct_headway<=30  pct_headway<=20
mode
Bus         5187          191          147             3.68             2.83
total       5187          191          147             3.68             2.83

So, I guess my questions are

• rather than just provide a warning that many would miss (ie. running analysis script, you'd only see warning if you look at log file), would it be better to raise an error or simply sys.exit with a clear warning in a case like Porto Alegre's?
  • I think maybe it would be better to sys.exit as it would present a nice error warning that users would be forced to deal with if they want to proceed (ie. manually clean GTFS, or omit it from analysis)
• in the case of Las Palmas, do you think its valid to include 'all stops' in the frequency analysis -- or was it right before to exclude stops that aren't in stop_times at all (if its not in there, is it serviced at all or just historical? - may be of no utility...)?

Keen to hear your thoughts @gigioLR , @shiqin-liu and @gboeing - apologies for long post, hope it makes sense,

Carl

[carlhiggs]

@marcdmallafre and @xavidelclos you might be interested in this thread on some issues with GTFS analysis where feeds don't include explicit departure times --- love to hear your thoughts on the questions asked above, on how we should handle these cases, as I know you've been working closely with this kind of data!

[carlhiggs]

Hi @gigioLR @marcdmallafre @xavidelclos @shiqin-liu @gboeing ,

I've just now addressed this issue (#359) by,

1. in the case of missing departure times in the stop_times.txt file, the analysis will exit and warn the user that this feed used as is will provide inaccurate results, and suggest one option of doing this being to interpolate departure times using stop sequence start and end times within each trip_id
2. providing an optional function to provide a linear interpolation of stop_times as described above.

This addresses the Porto Alegre situation described above, and may be relevant for GTFS feeds in other cities - perhaps Spanish ones.

In the case of Porto Alegre, wtih interpolation, here are the results:

# python _10_gtfs_analysis.py BR_Porto_Alegre_2023_R01

/home/ghsci/process/data/transit_feeds/poa/arquivo-gtfs.zip

**Note**: 5237 stops with null departure times found in stop_times.txt.
This GTFS feed has been configured to to fill null arrival and departure values.based on a linear interpolation according to the provided stop sequence start and end times within each trip_id.  This is an approximation based on the available information, but results may still differ from the actual service frequencies at these stops.

Porto Alegre analysis:
  - /home/ghsci/process/data/transit_feeds/poa/arquivo-gtfs.zip)
  - 20190930 to 20230723)
  - ['07:00:00', '19:00:00']


     Tram                  0/5187 (0.0%) stops identified with departure times (0.10 seconds)
     Metro                 0/5187 (0.0%) stops identified with departure times (0.09 seconds)
     Rail                  0/5187 (0.0%) stops identified with departure times (0.09 seconds)
     Bus                4973/5187 (95.9%) stops identified with departure times (9.78 seconds)
     Ferry                 0/5187 (0.0%) stops identified with departure times (0.15 seconds)
     Cable tram            0/5187 (0.0%) stops identified with departure times (0.10 seconds)
     Aerial lift           0/5187 (0.0%) stops identified with departure times (0.10 seconds)
     Funicular             0/5187 (0.0%) stops identified with departure times (0.09 seconds)
     Trolleybus            0/5187 (0.0%) stops identified with departure times (0.10 seconds)
     Monorail              0/5187 (0.0%) stops identified with departure times (0.10 seconds)

Porto Alegre summary (all feeds):
       tot_stops  headway<=30  headway<=20  pct_headway<=30  pct_headway<=20
mode
Bus         5187         3930         3372            75.77            65.01
total       5187         3930         3372            75.77            65.01

So, with interpolation we see that 3,372/5,187 (65%) of bus stops within Porto Alegre are estimated to have a day time week day service frequency of 20 minutes or better.

This is now addressed in https://github.com/global-healthy-liveable-cities/global-indicators/releases/tag/4.4.10

[gigioLR]

Hi @carlhiggs, I'll tell Júlio and the rest of our team about it and we'll test it here too. Thanks!

Giovani

[marcdmallafre]

Hi @carlhiggs and @gigioLR,

I will take a look at this. Probably i can reproduce the analysis with our cities. When I've calculated it with the new version I'll tell you if I've had important changes.

Thanks,
Marc