Skip to content

Instantly share code, notes, and snippets.

@geotheory
Created February 20, 2012 15:28
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save geotheory/1869673 to your computer and use it in GitHub Desktop.
Save geotheory/1869673 to your computer and use it in GitHub Desktop.
Simplifying spatial polygons in R {rgeos}
library(rgeos)
data(wrld_simpl)
plot(wrld_simpl, col='white', bg='grey9', border=NA)
limx <- c(-6,1)
limy <- c(49,60)
plot(wrld_simpl, xlim=limx, ylim=limy, col='white', bg='grey9', border=NA)
# Apply RDP algorithm with distance variable = 0.4 degrees
wrld_RDP <- gSimplify(wrld_simpl, 0.4, topologyPreserve=TRUE)
plot(wrld_RDP, xlim=limx, ylim=limy, col='white', bg='grey9', border=NA)
# SOME USEFUL TOOLS FOR ANALYSING THE RESULTING OBJECT
library(maptools)
str(wrld_simpl, max.level = 2)
str(wrld_simpl@data)
wrld_simpl$NAME
# Looking at id=207 (UK)
wrld_RDP@polygons[[207]]@ID # country code (ISO 3166-1 alpha-3)
wrld_simpl[[3]][207] # Country code
wrld_simpl[[5]][207] # Country name
length(wrld_RDP@polygons[[207]]@Polygons) # count sub-polygons within parent
wrld_RDP@polygons[[207]]@Polygons[[1]]@coords # list vertex coords of a child
wrld_RDP@polygons[[207]]@Polygons[[1]]@area # return its area (ie. sq degrees)
# Count child polygons per country id
subpolys <- list("")
for(i in 1:length(wrld_RDP)){
subpolys[[i]] <- length(wrld_RDP@polygons[[i]]@Polygons)
}
subpolys <- as.numeric(polys)
subpolys[1:3] # vector of child polygon counts per parent
# Useful method (thanks to 'The Praise of Insects' blog) for removing small polygons
plot(wrld_RDP, xlim=c(-6,-4), ylim=c(55,59), col='white', bg='grey9')
# Matrix of polygon areas
areas <- lapply(wrld_RDP@polygons, function(x) sapply(x@Polygons, function(y) y@area))
areas[1:3]
bigpolys <- lapply(areas, function(x) which(x > 0.1))
for(i in 1:length(bigpolys)){
if(length(bigpolys[[i]]) >= 1 && bigpolys[[i]][1] >= 1){
wrld_RDP@polygons[[i]]@Polygons <- wrld_RDP@polygons[[i]]@Polygons[bigpolys[[i]]]
wrld_RDP@polygons[[i]]@plotOrder <- 1:length(wrld_RDP@polygons[[i]]@Polygons)
}
}
plot(wrld_RDP, xlim=c(-6,-4), ylim=c(55,59), col='white', bg='grey9')
@ikashnitsky
Copy link

Line 39: Error: object 'polys' not found

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment