2

I am interested in learning about how to work with "Road Network" files - for example, I would like to find out the driving distance between two sets of geographical coordinates (i.e. longitude and latitude).

I found this shapefile for the Canadian Road Networks : https://www12.statcan.gc.ca/census-recensement/2011/geo/RNF-FRR/files-fichiers/lrnf000r22a_e.zip - now, I am trying to import this file into R.

Below, this is the code I am using to first download the shapefile to a temporary folder and then try to read the shapefile:

library(sf)
library(rgdal)
# Set the URL for the shapefile
url <- "https://www12.statcan.gc.ca/census-recensement/2011/geo/RNF-FRR/files-fichiers/lrnf000r22a_e.zip"

Create a temporary folder to download and extract the shapefile


temp_dir <- tempdir()
temp_file <- file.path(temp_dir, "lrnf000r22a_e.zip")


Download the shapefile to the temporary folder


download.file(url, temp_file)


Extract the shapefile from the downloaded zip file


unzip(temp_file, exdir = temp_dir)


Read the shapefile using the rgdal package


library(rgdal)
shapefile <- readOGR(dsn = temp_dir, layer = "lrnf000r22a_e")

But when trying to run the last line of code (readOGR), I get the following error:

OGR data source with driver: ESRI Shapefile 
Source: "C:\Users\me\AppData\Local\Temp\RtmpwDKofs", layer: "lrnf000r22a_e"
with 2246324 features
It has 21 fields
Integer64 fields read as strings:  OBJECTID 
Error: memory exhausted (limit reached?)
In addition: Warning messages:
1: OGR support is provided by the sf and terra packages among others 
2: OGR support is provided by the sf and terra packages among others 
3: OGR support is provided by the sf and terra packages among others 
4: OGR support is provided by the sf and terra packages among others 
5: OGR support is provided by the sf and terra packages among others 
6: OGR support is provided by the sf and terra packages among others

This seems to be a very large shapefile and my computer does not have enough memory to work with this file.

Has anyone ever encountered this problem before? Are there any strategies I can use to resolving this error?

Note: Perhaps its possible to read this in "chunks"? (e.g. Read N number of rows from shapefile using GeoPandas)

Vince
  • 20,017
  • 15
  • 45
  • 64
stats_noob
  • 145
  • 12
  • You loaded the sf package. It has at least one function for reading shapefiles, st_read. I’m not sure if it will be more efficient, but why not try that? – John Polo Apr 02 '23 at 16:30
  • Another method to consider is using PostGIS via the R API. https://cran.r-project.org/web/packages/rpostgis/index.html – John Polo Apr 02 '23 at 16:33
  • 2
    No shapefile can legally exceed 2Gb in the .shp or 2Gb in the .dbf, so most modern computers can actually handle all that could be loaded, though it's doubtful that they ever would. Keep in mind that shapefiles were created when 32Mb was a lot of RAM, and filesystems could not handle files in excess of 2Gb. "Very large" in modern parlance is 20-100 million features, which would blow out the maximum file size of a shapefile with complex geometry early on. See https://gis.stackexchange.com/questions/348557/how-does-shapefile-2gb-limit-equate-to-70-million-points – Vince Apr 02 '23 at 17:00

1 Answers1

3

Using the sf package you can apply an SQL SELECT query to a shapefile read using st_read with the query argument. For example if you have a huge road network to load you could maybe restrict it to main roads, or roads within a smaller bounding box that you are interested in.

For example to select only the highways:

hwy = st_read("./lrnf000r22a_e.shp",
       query="select * from lrnf000r22a_e where TYPE='HWY'")

That gets you 115619 features down from the 2.2 million in the full data set.

You can do a lot more filtering with st_read including using query= to select "chunks" by specifying ID numbers (I think this data set uses OBJECTID), and you can also pass a wkt_filter parameter to get only data within some polygonal bounds - see the examples in the help for st_read for full details.

Spacedman
  • 63,755
  • 5
  • 81
  • 115
  • @ Spacedman: Thank you so much for your answer! I guess I could also filter on : PRNAME_L (province)? – stats_noob Apr 03 '23 at 14:11
  • I think this is how you would use the code you wrote alongside mine: hwy = st_read(file.path(temp_dir, "lrnf000r22a_e.shp"), query="select * from lrnf000r22a_e where TYPE='HWY'") – stats_noob Apr 03 '23 at 14:14
  • Well yes, I had the shapefile downloaded to the working folder. Doesn't change anything with the select query though, does it? – Spacedman Apr 03 '23 at 15:51
  • @ Spaceman: Thank you so much for your answer! I posted a follow up question over here: https://gis.stackexchange.com/questions/456838/calculating-driving-distances-from-shapefiles ... can you please take a look at it if you have time? – stats_noob Apr 03 '23 at 19:38