# Simplifying and plotting polygons in Leaflet package in R?

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) %>%
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.
• Welcome to GIS SE! As a new user please take the tour to learn about our focused Q&A format.
– Midavalo
Apr 10, 2017 at 15:31
• do you want to retain at least one polygon per feature?
– Sam
Apr 11, 2017 at 8:51
• i'm think your problem has to do with the presence of holes in about 57 polygons, causing problems in the shapefile. You might be invalidating the geometries by removing the non-hole bit for some polygons. (bit of a guess).
– Sam
Apr 11, 2017 at 9:08

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.

• yeah this is pretty good alright, plotting time is no quicker (than my answer below) but the simplicity is way better.
– Sam
Apr 11, 2017 at 13:21

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) %>%
``````rsll <- spTransform(reduced.shapefile, CRS("+init=epsg:4326"))