3

Question:

I'm attempting to run some simple plots on some shapefiles with a high number of complex polygons. These aren't massive files - for the example below, the shapefile is 70mb (15872 obs), but plots, buffers, etc. are running extremely slowly, e.g. over an hour in some cases for a single relatively simple plot. Can anyone recommend some ways to optimise these operations? I cannot simplify the geometry of the polygons, as I need to ensure accuracy, but this seems to render R unworkable as a platform for reproducible geospatial research.

Sample data:

Download data as ESRI Shapefile from page at: https://gateway.snh.gov.uk/natural-spaces/dataset.jsp?dsid=SSSI

Reproducible example:

(install libraries and data ahead of execution)

require(RCurl) # used for fetching reproducible datasets
require(sf) # new simplefeature data class, supercedes sp in many ways
# using GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
require(sp) # needed for proj4string, deprecated by sf()
require(rgdal) # version version: 1.3-6
require(rgeos) # used for buffering below
require(devtools)
require(ggplot)

# Download data as ESRI Shapefile from page at: https://gateway.snh.gov.uk/natural-spaces/dataset.jsp?dsid=SSSI

unzip("SSSI_SCOTLAND_ESRI.zip")
sssi_sf <- st_read("SSSI_SCOTLAND.shp")
sssi_sp <- readOGR("./", "SSSI_SCOTLAND") 

# First test out plots using spatialfeatures and spdf with core R

system.time(
plot(sssi_sf)
  )

system.time(
plot(sssi_sp)
  )

# Then test out plots using spatialfeatures and spdf with ggplot2

system.time(
  ggplot() +
    geom_sf(data = sssi_sf)
)

system.time(
  ggplot() +
    geom_polygon(data = sssi_sp)
)

A few notes:

I've updated R and packages as of this morning, so am working with the latest versions (of pre-compiled binaries) of relevant libraries.

Also worth noting is that I've seen indications elsewhere that slow plots might be related to XQuartz limitations, but I've tested and confirmed that plots are terribly slow (> 1hr for plot(sssi_sf) below for example) on both MacOS and Windows (running in virtualbox).

  • 1
    I've added the R tag (I think a lot of us search on that tag) and removed the ggplot tag. – Spacedman Feb 8 at 17:08
1

Well, On my laptop with an i7-4core and 24GB RAM, I was able to use a primitive plot on the sp class in 7 seconds and all 7 columns in a panel plot with an sf class in 62 seconds.

However, ggplot2 seems to be having some speed issues with the data regardless of what I do. As to buffering, you do have 15,872 features but, additionally it is taking longer than expected (ran for 15 mins before I killed it).

You may want to try simplifying your geometry using: rgeos::gSimplify or sf::st_simplify. Some think that this reduces precision but, often what it accomplishes is reducing the number of unnecessary nodes in a feature geometry. This is dictated by the tolerance factor used in the generalization algorythm. Also, keep in mind that if the end result is the product of a buffering operation there is quite a bit of introduced uncertainty in the precision of the resulting buffers. Besides, maintaining high precision while plotting 15,872 features seems a bit unnecessary.

Here is a worked example on your data. It speeded up buffering operations to a reasonable expectation and rendered a ggplot2 plot in ~1 min.

library(sf)
library(sp)
library(rgeos)
library(ggplot2)
library(rgdal)

sssi_sf <- st_read("SSSI_SCOTLAND.shp")
sssi_sp <- readOGR(getwd(), "SSSI_SCOTLAND")

Check validity of feature(s) geometry

rgeos::gIsValid(sssi_sp) 
valid <- st_is_valid(sssi_sf)
  valid[valid == FALSE]

Here we apply the Douglas-Peuker simplify algorithm to sp and sf objects, then buffer the results to 20 meters with a quadseg of 30 (checking for run times). For the sp object we plot a single polygon to check if the shape has changed. After applying a simplify operation, the buffer operation took ~62 seconds on the sf class object and ~30 seconds on the sp.

sssi_sp2 <- rgeos::gSimplify(sssi_sp, tol=3)
  par(mfrow=c(1,2))
    plot(sssi_sp[1,])
    plot(sssi_sp2[1,])

  system.time(
    sssi_b1000 <- rgeos::gBuffer(sssi_sp2, width = 20, quadsegs = 30)
  )  

sssi_sf2 <- sf::st_simplify(sssi_sf)
  system.time(
    sssi_b1000 <- sf::st_buffer(sssi_sf2, dist = 20)
  ) 

Render ggplot2 plot on simplified geomerty

ggplot(sssi_sf2) + geom_sf(aes(fill = PA_CODE))

Just to illustrate the uncertainty that is introduce to the precision in buffering operations (across 4 scales), here is a quick example. Note that I am using an unnecessarily high number of line segments (quadsegs) to use to approximate the buffer.

par(mfrow=c(2,2))
  plot(gBuffer(sssi_sp[1,], width = 10,  
       quadsegs = 30), main="Buffer = 10")
    plot(sssi_sp[1,], add=TRUE) 
  plot(gBuffer(sssi_sp[1,], width = 20,  
       quadsegs = 30), main="Buffer = 20")
    plot(sssi_sp[1,], add=TRUE)
  plot(gBuffer(sssi_sp[1,], width = 40,  
       quadsegs = 30), main="Buffer = 40")
    plot(sssi_sp[1,], add=TRUE)
  plot(gBuffer(sssi_sp[1,], width = 80,  
       quadsegs = 30), main="Buffer = 80")
    plot(sssi_sp[1,], add=TRUE)
  • Your point regarding gSimplify/st_simplify is really helpful, as that sets some of my anxieties at ease re: accuracy. In the background, I'm also running calculations using st_within() with some relevant point data, so I'll do a quick check later on whether there are affects on the output of these calculations. But for plots, especially given the scale (as you suggest), this seems like a very good approach. Do you mind sharing code you used to run "primitive plot on the sp class in 7 seconds and all 7 columns in a panel plot with an sf class in 62 seconds"? – Jeremy Kidwell Feb 8 at 18:32
  • PS, how long did that ggplot on sssi_sf2 take you? I'm still working on it after 15 mins, even though it suggests 7 seconds when wrapped in system.time() – Jeremy Kidwell Feb 8 at 19:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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