1

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.
  • Welcome to GIS SE! As a new user please take the tour to learn about our focused Q&A format. – Midavalo Apr 10 '17 at 15:31
  • do you want to retain at least one polygon per feature? – Sam Apr 11 '17 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 '17 at 9:08
3

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 '17 at 13:21
1

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.

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.