2

Where can I find a list of all roads in Ontario?

I'm looking for a publicly available list, such as a website or downloadable PDF/spreadsheet. The intended audience aren't GIS professionals, so a GIS layer wouldn't be suitable.

Example use case:

Find all the roads in Ontario that have the word marble in their names.

1 Answer 1

1

If there is no publicly available list, let's prepare it :). The best source will be OSM Data (http://download.geofabrik.de/north-america/canada.html) from where we will take a bit. In below example I'm using R and osmextract package. Full list of available attributes for roads in OSM you can find under link https://wiki.openstreetmap.org/wiki/Key:highway.

Let's load the library and create subdirectory for files storage:

library(osmextract)
if (!dir.exists("data")) {
  dir.create("data")
}

Let's check, if there is a subset of data or Ontario

on <- oe_match(place = "Ontario")
#> The input place was matched with: Ontario

Let's download the file (here I'm using wget for download).

download.file(
  url = on$url, 
  destfile = "data/ontario-latest.osm.pbf",
  method = "wget",
  extra = "-c"
  )

In below steps we are extracting the necessary "layers" of geospatial data for further analysis:

oe_vectortranslate(
  file_path = "data/ontario-latest.osm.pbf",
  layer = "lines",
  extra_tags = c("highway", "name")
  )

oe_vectortranslate(
  file_path = "data/ontario-latest.osm.pbf",
  layer = "multipolygons",
  extra_tags = c("boundary", "admin_level")
)

Let's read the highways (roads) from the data file:

highways <- oe_read(
  file_path = "data/ontario-latest.osm.pbf",
  layer = "lines", 
  query = "SELECT osm_id, highway, name, geometry FROM lines WHERE highway IS NOT NULL"
  )
#> Simple feature collection with 1135393 features and 3 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -95.36875 ymin: 41.72912 xmax: -74.32076 ymax: 56.13301
#> Geodetic CRS:  WGS 84

This step is just for fun. Analyzing spatial data without drawing it is like decaffeinated coffee.

boundary <- oe_read(
  file_path = "data/ontario-latest.osm.pbf",
  layer = "multipolygons", 
  query = "SELECT osm_id, name, geometry FROM multipolygons \
             WHERE boundary = 'administrative' AND admin_level = '4'"
)

#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -95.15602 ymin: 41.67656 xmax: -74.32011 ymax: 56.86135
#> Geodetic CRS:  WGS 84


boundary |>
  sf::st_geometry() |>
  plot()

highways |>
  subset(highway %in% c("motorway", "primary")) |>
  sf::st_geometry() |>
  plot(col = "blue", add = TRUE)

Now, only the roads, which you are interested in:

highways |>
  subset(grepl("marble", name, ignore.case = TRUE))

#> Simple feature collection with 49 features and 3 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -79.82654 ymin: 43.40498 xmax: -75.28377 ymax: 45.54217
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>         osm_id      highway                 name                       geometry
#> 3950   5064974 unclassified     Marble Rock Road LINESTRING (-76.1998 44.373...
#> 4975   5379107  residential Marble Arch Crescent LINESTRING (-79.29217 43.73...
#> 4977   5379116  residential Marble Arch Crescent LINESTRING (-79.29426 43.73...
#> 7042   9397164  residential      Marblehead Road LINESTRING (-79.55451 43.69...
#> 31503 28615578  residential   Marblethorne Court LINESTRING (-79.59533 43.64...
#> 39875 31382638  residential         Marble Place LINESTRING (-79.48716 44.06...
#> 56675 33852914  residential  Marblehead Crescent LINESTRING (-79.74659 43.73...
#> 61350 33873128  residential  Marbleseed Crescent LINESTRING (-79.75081 43.75...
#> 88686 33974318  residential         Marble Court LINESTRING (-77.66915 44.46...
#> 91993 33979087     tertiary    Marble Point Road LINESTRING (-77.69978 44.48...

To write it down (without geometries) you can use write.csv() function, like:

highways |>
  sf::st_drop_geometry() |>
  subset(grepl("marble", name, ignore.case = TRUE)) |>
  write.csv(file = "data/marble_roads.csv")

To add other attributes, like number of lanes, surface just modify extra_tags = c("highway", "name") in oe_vectortranslate() and then modify the query in highways <- oe_read(...) chunk.

Regards, Grzegorz

Created on 2022-10-04 with reprex v2.0.2

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.