sf: Session aborts when reading large geopackage into R

Question

I want to read a geopackage into R that is roughly ~ 600 mb in size. Every time I want to read it my R-Studio session aborts when doing the

data = read_sf(path_geo)

Is there any "trick" that I could use to get the file somehow into memory?

The file can be downloaded from here: https://drive.google.com/file/d/1HPR7XpYJdddACqsGTVOZ2rwAI2NyMrD0/view?usp=sharing

Answer 1

Here's how to get a subset of a large geopackage layer if you can do a simple filter by attribute:

This layer has 149 features:

> p = st_read("flu.gpkg","localauth",quiet=TRUE)
> dim(p)
[1] 149   4

Suppose that layer is too big to read without crashing my computer, but I only want the features where the population is over 1 million. I can use this SQL query expression, and only get back 7 features:

> p = st_read("flu.gpkg",query="select * from localauth where pop2016 > 1000000", quiet=TRUE)
> dim(p)
[1] 7 4

which is a spatial object as expected:

> p
Simple feature collection with 7 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 330453.6 ymin: 89700.1 xmax: 640168.7 ymax: 482748.7
Projected CRS: unnamed
       CODE          NAME pop2016                           geom
1 E08000025    Birmingham 1124569 MULTIPOLYGON (((399306.1 27...
2 E10000012         Essex 1455340 MULTIPOLYGON (((575268.1 18...
3 E10000014     Hampshire 1360426 MULTIPOLYGON (((431521.1 90...
4 E10000015 Hertfordshire 1176720 MULTIPOLYGON (((503035.8 19...
5 E10000016          Kent 1541893 MULTIPOLYGON (((633197.5 15...
6 E10000017    Lancashire 1198798 MULTIPOLYGON (((338034.1 43...
7 E10000030        Surrey 1176549 MULTIPOLYGON (((500486.3 17...

Note that the SQL is run at the GDAL driver level, so the query is done by efficient code that doesn't have to load all the data in. If the column you are querying on has a database index, it could also be really quick.

If you want to get a quick look at a large geopackage layer, just read the first few rows:

> p = st_read("flu.gpkg",query="select * from localauth limit 5", quiet=TRUE)
> p
Simple feature collection with 5 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 350784.7 ymin: 152770 xmax: 551943.8 ymax: 413058.6
Projected CRS: unnamed
       CODE                         NAME pop2016                           geom
1 E09000002         Barking and Dagenham  206460 MULTIPOLYGON (((543732.7 18...
2 E09000003                       Barnet  386083 MULTIPOLYGON (((528648.4 19...
3 E08000016                     Barnsley  241218 MULTIPOLYGON (((416179.2 39...
4 E06000022 Bath and North East Somerset  187751 MULTIPOLYGON (((379810.8 16...
5 E06000055                      Bedford  168751 MULTIPOLYGON (((497192.9 24...

You can also use various spatial SQL functions in the query. Maybe I only want large polygons:

> p = st_read("flu.gpkg",query="select * from localauth where st_area(geom) > 5000000000", quiet=TRUE)
> p
Simple feature collection with 6 features and 3 fields
Geometry type: MULTIPOLYGON

Note st_area here is not the R function of the same name!

I think you can possibly do complex spatial intersection tests in the query here, for example if I had another layer of rivers in the same geopackage, I think I could select areas that have a given river running through, using an SQL query, all running in the driver code and without the overhead of reading all the package into R.