[GIS] Simplifying and plotting polygons in Leaflet package in R

leafletpolygonrrgdalrgeos

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/abs@.nsf/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.

Best Answer

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.

Related Question