In my research I frequently work with large datasets. Sometimes that means datasets that cover the entire globe, and other times it means working with lots of micro-level event data. Usually, my computer is powerful enough to load and manipulate all of the data in R without issue. When my computer’s fallen short of the task at hand, my solution has often been to throw it at a high performance computing cluster. However, I finally ran into a situation where the data proved too large even for that approach.

As a result, I finally had to teach myself how to break large spatial datasets into more manageable chunks. In the process a learned a little SQL and a lot about the underlying software libraries that power the r-spatial ecosystem of R packages. In this post, I walk through the workflow I developed for this task and explain the logic behind each step.

On disk

The general idea is to work with data ‘on disk’ instead of ‘in memory’. Normally, when you load a dataset into R, your computer reads it from whatever storage media it uses (hard drive or solid state drive) into memory (RAM). Memory is considerably faster to read from and write to than storage, which is what lets you complete simple operations in R in the blink of an eye. Most consumer computers have much more storage than RAM (my 2015 MacBook Pro has 256 GB of storage and 8 GB of memory) so it’s very possible to end up with a dataset larger than your computer’s memory. In fact, it doesn’t have to be anywhere near the size of your computer’s memory to bump into this limit because every other application you have running uses up memory as well.

To deal with this issue, you can extract just the parts of a dataset you need to work with at any given time; this subset will be loaded into memory, and the rest remain on disk and invisible to R¹. There are a couple of R packages that exist for dealing with this issue, such as bigmemory for basic R data types like numerics or disk.frame for dplyr-compatible operations, but neither supports spatial data.

I’m going to use the cshapes to illustrate and explain this workflow². You can download and extract it from within R:

## download cshapes dataset
download.file('http://downloads.weidmann.ws/cshapes/Shapefiles/cshapes_0.6.zip',
              'cshapes.zip')

## extract cshapes dataset
unzip('cshapes.zip')

## check that dataset extracted correctly
list.files(path = '.', pattern = 'cshapes')

## [1] "cshapes_shapefile_documentation.txt" "cshapes.dbf"                        
## [3] "cshapes.prj"                         "cshapes.shp"                        
## [5] "cshapes.shx"                         "cshapes.zip"

Then use the sf package to load the data and check them out:

## load packages
library(tidyverse)
library(sf)

## read in cshapes
cshapes <- st_read('cshapes.shp')

The cshapes dataset is specifically designed to be easy to load and manipulate on a conventional laptop computer. To do this, it sacrifices a significant degree of detail in the polygons that represent each individual state. For many analyses, this is fine and won’t affect the results. However, sometimes you need to measure the length of borders between states, and the coastline paradox dictates that you use the most high resolution spatial data possible. In that case, the data might be too large for your computer to hold in memory. If that’s the case, then it’s time to start thinking about leaving the data on disk and only loading what you really need at any given point.

SQL

Luckily the sf package supports SQL queries to filter the data on disk and only read in a subset of the total data. SQL is a language for interacting with relational databases, and is incredibly fast compared to loading data into R and then filtering it. SQL has many variants, referred to as dialects, and the sf package uses one called OGR SQL dialect to interact with spatial datasets. The basic structure of a SQL call is SELECT col FROM "table" WHERE cond.

• SELECT tells the database what columns (fields in SQL parlance) we want
• FROM tells the database what table (databases can have many tables) to select those columns from
• WHERE tells that database we only want rows where some condition is true

If you use the tidyverse a lot, this may seem familiar to you because it’s pretty similar to dplyr syntax, except dplyr already knows which data frame you want to work with. If we want to only load one polygon at a time into R, then we need to know the field (or combination of fields) that uniquely identifies a polygon. To demonstrate, let’s load just the polygon for Morocco that begins in 1976 when it annexed the Northern part of Western Sahara. Let’s cheat by looking at the data I’ve loaded into R:

## filter to Morocco beginning in 1976
cshapes %>% filter(CNTRY_NAME == 'Morocco', GWSYEAR == 1976)

## Simple feature collection with 1 feature and 24 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -15.22687 ymin: 23.11465 xmax: -1.011809 ymax: 35.91916
## geographic CRS: WGS 84
##   CNTRY_NAME     AREA CAPNAME CAPLONG CAPLAT FEATUREID COWCODE COWSYEAR COWSMONTH COWSDAY COWEYEAR
## 1    Morocco 576351.8   Rabat   -6.83  34.02       220     600     1976         4       1     1979
##   COWEMONTH COWEDAY GWCODE GWSYEAR GWSMONTH GWSDAY GWEYEAR GWEMONTH GWEDAY ISONAME ISO1NUM ISO1AL2
## 1         8       4    600    1976        4      1    1979        8      4 Morocco     504      MA
##   ISO1AL3                       geometry
## 1     MAR MULTIPOLYGON (((-4.420418 3...

The cshapes dataset records when states change territorial boundaries or capital locations, so the combination of a state name or identifier and a start or end date uniquely identifies all rows in the data. Since, this polygon begins on April 1, 1976 and the Gleditsch and Ward code for Morocco is 600, plugging it all into the query argument to st_read() gets us:

## read in morocco polygon
morocco <- st_read('cshapes.shp',
                   query = 'SELECT * FROM "cshapes" WHERE GWCODE = 600 AND GWSYEAR = 1976 AND GWSMONTH = 4 AND GWSDAY = 1')

## verify country name
morocco$CNTRY_NAME

## [1] "Morocco"

Awesome! We were able to read in just one polygon from the cshapes dataset. Note that * means all columns. As I mentioned above, this is cheating, since we had to read the whole dataset into R with a standard st_read() call to learn the names and values of the variables we then filtered on.

Sneaking a peek

When this isn’t an option, we can sneak a peak at the data by loading just the first observation into R. This requires significantly less memory than loading an entire dataset, and can give us the information we need to filter the full dataset and read in one observation at a time. Most SQL implementations don’t have row numbers, so it’s hard to just grab the first row of the data for this purpose. However, the OGR SQL dialect documentation notes that it implements a special field called FID that is a feature ID, i.e., a row number. We can take advantage of FID to select the first polygon from the data using the query argument to st_read() again:

## read in first row of the data
cshapes_row <- st_read('cshapes.shp', query = 'SELECT * FROM "cshapes" WHERE FID = 1')

## inspect
cshapes_row

## Simple feature collection with 1 feature and 24 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -58.0714 ymin: 1.836245 xmax: -53.98612 ymax: 6.001809
## geographic CRS: WGS 84
##   CNTRY_NAME     AREA    CAPNAME CAPLONG   CAPLAT FEATUREID COWCODE COWSYEAR COWSMONTH COWSDAY
## 1   Suriname 145952.3 Paramaribo   -55.2 5.833333         1     115     1975        11      25
##   COWEYEAR COWEMONTH COWEDAY GWCODE GWSYEAR GWSMONTH GWSDAY GWEYEAR GWEMONTH GWEDAY  ISONAME
## 1     2016         6      30    115    1975       11     25    2016        6     30 Suriname
##   ISO1NUM ISO1AL2 ISO1AL3                 _ogr_geometry_
## 1     740      SR     SUR POLYGON ((-55.12796 5.82217...

Even if we knew that the data had an ID column and start and end dates, we wouldn’t know the precise formatting (capitalization, underscores or dashes) of column names, or whether start and end dates are stored as one column or sets of three like they are here.

Making a list

We still need more information if we want to iterate through the polygons in the data and load them one at a time. We know what columns uniquely identify the rows, but what don’t know all the values they take on. Without that, we we’re stuck. What (usually) makes spatial data big is not the tabular data themselves, but the spatial features they’re attached to. This is particularly the case with polygons, which can be incredibly large in size for complex features. So, the goal here is to get the data we care about (ID column and start date) and ditch everything else, loading only the bare minimum into memory.

To do this, we’ll use the ogr2ogr() function in the gdalUtils package³. ogr2ogr() converts between different spatial data formats. It also offers two features that we’re going to use to cut down the data to the bare minimum. The select argument is a SQL selection, so we’re going to create a comma separated list of our key columns. The nlt argument specifies what type of geometry to create in the output. Conveniently it accepts NONE as a value, which will yield a plain table of data with none of the memory-hogging geometries:

## load package
library(gdalUtils)

## convert to nonspatial geometry
ogr2ogr(src_datasource_name = 'cshapes.shp', dst_datasource_name = 'cshapes_no_geom',
        select = 'GWCODE,GWSYEAR,GWSMONTH,GWSDAY', nlt = 'NONE')

This will create a shapefile in the new directory cshapes_no_geom called cshapes. The usual .shp and .shx components of a shapefile are missing, but the .dbf part is there, and that’s the one we care about. Load it up with st_read() and we’ll have what we need:

## load non-geometry table
cshapes_id <- st_read('cshapes_no_geom/cshapes.dbf')

## inspect
head(cshapes_id)

##   GWCODE GWSYEAR GWSMONTH GWSDAY
## 1    110    1966        5     26
## 2    115    1975       11     25
## 3     52    1962        8     31
## 4    101    1946        1      1
## 5    990    1962        1      1
## 6    972    1970        6      4

Now you can load polygons one at a time and perform whatever geometric operations you need to. To illustrate, I’ll load the first four polygons in the dataset, calculate their area, and then plot them.

## set up four panel plot
par(mfrow = c(1, 4), mar = c(6.1, 4.1, 4.1, 4.1))

## read in each polygon and plot 
for (i in 1:4) {
  
  ## build SQL query
  query_str <- str_c('SELECT * FROM "cshapes" WHERE GWCODE = ', cshapes_id$GWCODE[i],
                     ' AND GWSYEAR = ', cshapes_id$GWSYEAR[i],
                     ' AND GWSMONTH = ', cshapes_id$GWSMONTH[i],
                     ' AND GWSDAY = ', cshapes_id$GWSDAY[i])
  
  ## read in data
  pol <- st_read('cshapes.shp', query = query_str)
  
  ## plot data
  pol %>%
    st_geometry() %>% 
    plot(main = pol$CNTRY_NAME,
         sub = str_c(round(units::set_units(st_area(pol), 'km^2'), digits = 0),
                      ' km^2'))
  
}

Won’t you be my neighbor?

Sometimes (oftentimes in spatial analysis) we need not just a polygon, but also its neighbors. That means loading just one polygon is insufficient. If your data are already in R, this is easy with the st_filter() function, but it’s much trickier if you’re trying to filter data before loading them into R⁴. Luckily, st_read() as you covered! The wkt_filter accepts a well-known text string that can be used to filter the data before loading them into R⁵. Well-known text is a standard string representation of geometry, and is actually how the sf package prints geometry in R:

st_point(c(1, 2))

We want to use the wkt_filter argument to only load polygons that intersect with our Morocco polygon into R. To do that, we need to convert our polygon to a well-known text string with the st_as_text() function, then pass it to st_read(). However, st_as_text() only accepts sfc and sfg objects, not sf objects:

## create well known text object to filter cshapes on disk
morocco_wkt <- st_as_text(morocco)

## Error in UseMethod("st_as_text"): no applicable method for 'st_as_text' applied to an object of class "c('sf', 'data.frame')"

To get around this, we need to drop the data on morocco and extract just the geometry of the polygon with st_geometry():

## create well known text object to filter cshapes on disk
morocco_wkt <- morocco %>% 
  st_geometry() %>% # convert to sfc
  st_as_text() # convert to well known text

## plot morocco and neighbors
st_read('cshapes.shp', wkt_filter = morocco_wkt) %>%
  st_geometry() %>%
  plot(main = morocco$CNTRY_NAME)

## add morocco polygon on top
morocco %>% 
  st_geometry() %>%
  plot(add = T, col = rgb(0, 1, 0, .5))

Notice that there are multiple polygon boundaries within the green area of our green Morocco polygon. That’s because there are 4 Morocco polygons in the data starting in 1956, 1958, 1976, and 1979. Be sure to filter the dataset, either as part of the SQL query or in a dplyr::filter() so that you only get polygons that existed contemporaneously with your polygon of interest.

Wrapping up

So far, we’ve covered:

• How to extract the first polygon for a spatial dataset and learn the names of identifier columns
• How to strip the geometry from a spatial dataset and extract just a table of these columns
• How to use these columns to iterate through the polygons in the dataset and import them one at a time, or along with their neighbors

You can technically skip the first two steps and just move the .shp and .shx files out of the directory before loading the .dbf file with st_read(), but that kind of feels like cheating to me⁶ and it only works with shapefiles. If you have another type of spatial dataset, read on.

This time for real

In my research, I often need to work with spatial data that’s measured at or aggregated up to different administrative divisions (ADMs). GADM helpfully provides a global dataset of ADMs. Although you can download ADMs for specific countries, I work with data in enough different countries that I finally decided to just download the entire dataset. While the cshapes example above just illustrated how to implement a pipeline for working with spatial data on disk, you may actually need to use one with these data depending on your machine’s hardware.

This master dataset comes as a GeoPackage. Most importantly for us, that means we can’t just delete a few component files to load the non-spatial table from the dataset; we have to convert it from a spatial dataset to a non-spatial one with ogr2ogr(). The GeoPackage contains ADMs from level 0 (countries) all the way down to level 5. Each level is stored as a separate layer in the .gpkg, and we can get a list of available layers with the st_layers() function:

## get layers
st_layers('~/Dropbox/Datasets/GADM/gadm34_levels_gpkg/gadm34_levels.gpkg')

## Driver: GPKG 
## Available layers:
##   layer_name geometry_type features fields
## 1     level0 Multi Polygon      256      2
## 2     level1 Multi Polygon     3610     10
## 3     level2 Multi Polygon    45962     13
## 4     level3 Multi Polygon   147427     16
## 5     level4 Multi Polygon   138053     14
## 6     level5 Multi Polygon    51427     15

We want to work with the third-order administrative divisions (cities, towns, and other municipalities in the US context), so we need the level3 layer. Where we just used the name of the dataset in our SQL call before, this time we’ll use level3. Now we just follow the same workflow as with the cshapes dataset above:

## get first observation
level3 <- st_read('~/Dropbox/Datasets/GADM/gadm34_levels_gpkg/gadm34_levels.gpkg',
                  query = 'SELECT * FROM "level3" WHERE FID = 1', layer = 'level3')

## inspect
level3

## Simple feature collection with 1 feature and 16 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 13.08792 ymin: -8.010127 xmax: 13.59943 ymax: -7.708598
## geographic CRS: WGS 84
##   GID_0 NAME_0   GID_1 NAME_1 NL_NAME_1     GID_2 NAME_2 NL_NAME_2       GID_3 NAME_3 VARNAME_3
## 1   AGO Angola AGO.1_1  Bengo      <NA> AGO.1.1_1 Ambriz      <NA> AGO.1.1.1_1 Ambriz      <NA>
##   NL_NAME_3  TYPE_3 ENGTYPE_3 CC_3 HASC_3                           geom
## 1      <NA> Commune   Commune <NA>   <NA> MULTIPOLYGON (((13.12764 -7...

This time we have a single column that uniquely identifies observations, GID_3, so we only have to extract one column from the dataset. We use the ogr2ogr() function as before, but we have to specify the layer = 'level3' argument since the GeoPackage has more than one layer and we want to work with a specific one. Since GID_3 is our identifier column, that’s what we select from the dataset:

## convert to nonspatial geometry
ogr2ogr(src_datasource_name = '/Users/Rob/Dropbox/Datasets/GADM/gadm34_levels_gpkg/gadm34_levels.gpkg',
        dst_datasource_name = 'gadm34_levels_no_geom',
        layer = 'level3',
        select = 'GID_3',
        nlt = 'NONE')

## load non-geometry table
gadm_ids <- st_read('gadm34_levels_no_geom/level3.dbf')

## inspect
head(gadm_ids)

##         GID_3
## 1 AGO.1.1.1_1
## 2 AGO.1.1.2_1
## 3 AGO.1.1.3_1
## 4 AGO.1.2.1_1
## 5 AGO.1.2.2_1
## 6 AGO.1.2.3_1

And we can again read the polygons into R one at a time and perform whatever spatial operations we need. Since our identifying column is a string this time, we need to enclose it quotes in our SQL call. SQL is very picky about quotation mark types, so while we needed to surround our layer name with double quotes, we need to surround our identifier variable with single quotes. I’m already using single quotes to define the character string for the SQL call, so I need to escape the single quotes around the identifier. You can do this with a single backslash (\). Thus, you can include single quotes in a single-quoted string like this: 'this is a string \'this is another part of a string\''. Other than that wrinkle, things are pretty much the same as with cshapes:

## for reproducibility
set.seed(27599)

## set up four panel plot
par(mfrow = c(1, 4), mar = c(2.1, 4.1, 4.1, 4.1))

## read in each polygon and plot
for (i in sample(1:nrow(gadm_ids), 4, replace = F)) { # mix it up
  
  ## build SQL query
  query_str <- str_c('SELECT * FROM "level3" WHERE GID_3 = \'',
                     gadm_ids$GID_3[i], '\'')
  
  ## read in polygon for ADM3 i
  adm3 <- st_read('~/Dropbox/Datasets/GADM/gadm34_levels_gpkg/gadm34_levels.gpkg',
                  query = query_str, layer = 'level3')
  
  ## plot polygon and label with full name
  print(plot(adm3$geom,
             main = adm3 %>%
               select(starts_with('NAME_')) %>% # get all name variables
               st_drop_geometry() %>% # drop geometry
               rev() %>% # reverse order of names to 3, 2, 1, 0
               str_c(collapse = ', '), # collapse w/ commas
             cex.main = .6))
  
}

Spatially filtering the GADM dataset is just as easy as with cshapes. To illustrate, I’m going to pull out a random polygon and use it to filter the data. However, these are third-order administrative divisions, and so it’s possible that even capturing all adjacent polygons won’t cover a very large area. To deal with this concern, we can buffer the polygon with the st_buffer() function before we convert it to well-known text:

## import single polygon
adm3 <- st_read('~/Dropbox/Datasets/GADM/gadm34_levels_gpkg/gadm34_levels.gpkg',
                  query = str_c('SELECT * FROM "level3" WHERE FID = 63130'))

## create well known text object to filter GADM on disk
adm3_wkt <- adm3 %>% 
  st_geometry() %>% # convert to sfc
  st_buffer(.025) %>% # buffer .05 decimal degrees
  st_as_text() # convert to well known text

## plot Dakkoun and neighbors w/in .05 decimal degrees
st_read('~/Dropbox/Datasets/GADM/gadm34_levels_gpkg/gadm34_levels.gpkg',
        layer = 'level3', wkt_filter = adm3_wkt) %>%
  st_geometry() %>%
  plot(main = adm3 %>%
               select(starts_with('NAME_')) %>%
               st_drop_geometry() %>%
               rev() %>%
               str_c(collapse = ', '))

## plot Dakkoun and highlight
adm3 %>%
  st_geometry() %>%
  plot(add = T, col = 'green')

## plot buffered polygon used to filter GADM on disk
adm3 %>% 
  st_geometry() %>% 
  st_buffer(.025) %>% 
  st_cast('LINESTRING') %>%
  plot(add = T, col = 'blue')

The green polygon above is Dakkoun, the 63,130th polygon in the the dataset. The blue line is the extent of the .025 decimal degree buffer applied to it to before filtering the dataset. This workflow can speed things up when working with these data, considering there are 884,562 third-order administrative division polygons in the dataset.

Making data manageable

The query and wkt_filter arguments to st_read() can help you work with large spatial datasets that are either too big to load into memory, or too slow to work with once loaded. While this is less of a concern with low resolution datasets created by social scientists, it can be incredibly useful if you ever have to work with super high resolution data created by remote sensing technologies or actual cartographers and geographers.