Segement a SpatialPolygons or SpatialPolygonsDataFrame

SegmentSpatialPolygons.r

## package for handling spatial data:
library(sp)
## package that contains simple world map data
library(maptools)
## package for modifying spatial data:
library(rgeos)
## package used for setting up clipping areas
library(raster)

## function that segments the polygons of SpatialPolygons(DataFrame) object
## 'object' by a factor of 'segments', where 'segments' should be an
## integer value >= 1
SPSegmentation <- function(object, segments) {

  ## apply to each element of the SpatialPolygons object:
  object@polygons <- lapply(object@polygons, function(pol) {

    ## apply to each polygon of each element of each SpatialPolygons object:
    pol@Polygons <- lapply(pol@Polygons, function(p) {

      ## apply linear interpolation to the coordinates:
      p@coords <- apply(p@coords, 2, function(x) {
        approx(x, n = (nrow(p@coords) - 1)*segments + 1)$y
      })
      return(p)

    })
    return(pol)

  })
  return(object)
}

data(wrld_simpl)

## first clip the world data to an area useful for our example...

## start be defining the clipping area:
clip <- as(extent(-10, 15, 50, 60), "SpatialPolygons")
## add the proper projection to the data:
proj4string(clip) <- proj4string(wrld_simpl)
## and clip it:
wrld_clip <- gIntersection(wrld_simpl, clip, byid=TRUE)

## setup a simple rectangle:
simple_rect <- as(extent(-8, 13, 51, 59), "SpatialPolygons")
## add the proper projection to the data:
proj4string(simple_rect) <- proj4string(wrld_simpl)

## transform the clipped world map to UTM:
wrld_clip_trans   <- spTransform(wrld_clip,   CRS("+proj=utm +zone=31 +datum=WGS84"))
## transform the simple rectangle without segmentation:
simple_rect_trans <- spTransform(simple_rect, CRS("+proj=utm +zone=31 +datum=WGS84"))

## now first segment the simple rectangle (increasing the
## number of vertices by a factor of 10)
segmented_rect <- SPSegmentation(simple_rect, 10)
## transform the segmented rectangle to UTM:
segmented_rect_trans <- spTransform(segmented_rect, CRS("+proj=utm +zone=31 +datum=WGS84"))

## plot everything:
png("segmentation-demo.png", 1200, 400, res = 300, type = "cairo")
par(mfcol=c(1,3), mar = c(0, 0, 2, 0))
plot(wrld_clip, col = "lightgreen", main = "original projection")
plot(simple_rect, add = T, border = "red")
plot(wrld_clip_trans, col = "lightgreen", main = "new projection\nno segmentation")
plot(simple_rect_trans, add = T, border = "red")
plot(wrld_clip_trans, col = "lightgreen", main = "new projection\nsegmentation")
plot(segmented_rect_trans, add = T, border = "red")
dev.off()