Simplifying and plotting polygons in Leaflet package in R?

Question

I was finding the plotting of shapefiles very slow in R. After reading this (How to speed up the plotting of polygons in R?) I found all the tips were helpful for plotting in base R. The second one (creating a reduced shapefile by removing small polygons using a custom function) was particularly useful.

library(rgdal)
library(leaflet)
library(rgeos)

full.shapefile <- readOGR(dsn="170411_maps", layer="SA2_2011_AUST")

simplified.shapefile <- gSimplify(full.shapefile, tol=0.01, topologyPreserve=TRUE)
simplified.shapefile = SpatialPolygonsDataFrame(simplified.shapefile, data=full.shapefile@data)

# Remove the islands
getSmallPolys <- function(poly, minarea=0.01) {
  # Get the areas
  areas <- lapply(poly@polygons,
              function(x) sapply(x@Polygons, function(y) y@area))

  # Quick summary of the areas
  print(quantile(unlist(areas)))

  # Which are the big polygons?
  bigpolys <- lapply(areas, function(x) which(x > minarea))
  length(unlist(bigpolys))

  # Get only the big polygons
  for(i in 1:length(bigpolys)){
    if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]][1] >= 1){
      poly@polygons[[i]]@Polygons  <- poly@polygons[[i]]@Polygons[bigpolys[[i]]]
      poly@polygons[[i]]@plotOrder <- 1:(length(poly@polygons[[i]]@Polygons))
    }
  }
  return(poly)
}

reduced.shapefile <- getSmallPolys( simplified.shapefile , 0.01 )

This worked perfectly as long as a stayed with the basic R plot() function.

However, shifting to the Leaflet package caused problems. While the other SpatialPolygonsDataframes would work fine (albeit very slowly) in Leaflet, the reduced.shapefile would not.

leaflet(data = reduced.shapefile) %>%
  addPolygons(fillColor = "blue",
              fillOpacity = 1,
              color = "transparent",
              weight = 1)

Running that code would give the error: Error in pgons@Polygons[[index]] : subscript out of bounds.

Can anyone see the specific cause of the problem, and a way that I can restructure or adjust the function (getSmallPolys) such that the new SpatialPolygonsDataframe will work in Leaflet and not just plot()?

I've included some additional notes below:

• If anyone's interested, the SA2 ESRI shapefile I'm using can be downloaded here: http://www.abs.gov.au/AUSSTATS/[email protected]/DetailsPage/1270.0.55.001July%202011?OpenDocument
• Examples of the specific subgroups of polygons that cause the problems include reduced.shapefile[156,] and reduced.shapefile[981,]
• I suspect that the problem relates to a part of Leaflet trying to reference an element of i in reduced.shapefile[981,]@polygons[[1]]@Polygons[[i]] that is larger than the number of polygons in the reduced shapefile.

Answer 1

Recommend you try rmapshaper::ms_simplify to prepare simplified polygon data for plotting, it should handle everything correctly, and you'll only need one line of code.

Answer 2

i'm not fully sure as to your error message but I think it is to do with the presence of holes in the shapefile features and that the function does not take them into account (therefore invalidating some geometries or 'messing' with the internal structure of the features after processing). Sorry for the lack of detail here as it's just a hunch.

Personally, i'd disaggregate all multipart polygons and process individually:

packs <- c("sp","rgdal","rgeos","data.table","leaflet")
lapply(packs, require, character.only = TRUE)

# get data
full.shapefile <- readOGR(dsn="170411_maps", layer="SA2_2011_AUST")

simplified.shapefile <- gSimplify(full.shapefile, tol=0.01, topologyPreserve=TRUE)
simplified.shapefile = SpatialPolygonsDataFrame(simplified.shapefile, data=full.shapefile@data)

something like this:

getSmallPolys <- function(poly, minarea=0.01) {

  # disaggregate multipart polygons
  d <- disaggregate(simplified.shapefile)

  # add true area (holes taken into account) & provide a unique ID to each part (needed later)
  d@data$area <- gArea(d, byid = TRUE)
  d@data$uniq <- 1:nrow(d@data)

  # convert to data table for subsetting
  dt <- as.data.table(d@data)

  # add a count of polygons per original feature ID
  dt[,COUNT:=.N,by=.(SA2_MAIN11)]

  # add a flag that states either over minarea or max of its ID
  dt[,flag:=ifelse(area>minarea | area==max(area),1,0),by=.(SA2_MAIN11)]

  # subset to keep only records with flag
  dt2 <- dt[flag==1]

  # subset the disaggregated shapefile, keeping only original column attributes
  d2 <- d[d@data$uniq %in% dt2[,uniq],names(sa)]

  # reaggregate to multipart polygons
  a <- aggregate(d2,by=names(d2))

  return(a)
}

reduced.shapefile <- getSmallPolys( simplified.shapefile , 0.01 )

I've made the assumption that you want to keep at least one polygon per SA2_MAIN11 ID, even if it's smaller than 0.01, otherwise you lose most of the data.

This seems to now plot happily in leaflet, down from 55 to 60 secs on the original to about 6 secs:

leaflet(data = reduced.shapefile) %>%
  addPolygons(fillColor = "blue",
              fillOpacity = 1,
              color = "transparent",
              weight = 1)

Finally, by transforming the result to EPSG:4326 (the native CRS of leaflet), we can shave another 5 to 10% off the plot time:

rsll <- spTransform(reduced.shapefile, CRS("+init=epsg:4326"))

hopefully something useful in here, sorry i couldn't solve exactly why your function caused errors.