Building on the excellent work of Paul Regular, here is a version that should ensure exclusive polygons (i.e. no overlapping).
I've added a new argument fd for fairy dust to address an issue I discovered working with UTM-type coordinates. Basically as I understand the algorithm works by sampling lateral points from the contour lines to determine which side is inside the polygon. The distance of the sample point from the line can create problems if it ends up in e.g. behind another contour. So if your resulting polygons looks wrong try setting fd to values 10^±n until it looks very wrong or about right..
raster2contourPolys <- function(r, levels = NULL, fd = 1) {
## set-up levels
levels <- sort(levels)
plevels <- c(min(values(r)-1, na.rm=TRUE), levels, max(values(r)+1, na.rm=TRUE)) # pad with raster range
llevels <- paste(plevels[-length(plevels)], plevels[-1], sep=" - ")
llevels[1] <- paste("<", min(levels))
llevels[length(llevels)] <- paste(">", max(levels))
## convert raster object to matrix so it can be fed into contourLines
xmin <- extent(r)@xmin
xmax <- extent(r)@xmax
ymin <- extent(r)@ymin
ymax <- extent(r)@ymax
rx <- seq(xmin, xmax, length.out=ncol(r))
ry <- seq(ymin, ymax, length.out=nrow(r))
rz <- t(as.matrix(r))
rz <- rz[,ncol(rz):1] # reshape
## get contour lines and convert to SpatialLinesDataFrame
cat("Converting to contour lines...\n")
cl0 <- contourLines(rx, ry, rz, levels = levels)
cl <- ContourLines2SLDF(cl0)
## extract coordinates to generate overall boundary polygon
xy <- coordinates(r)[which(!is.na(values(r))),]
i <- chull(xy)
b <- xy[c(i,i[1]),]
b <- SpatialPolygons(list(Polygons(list(Polygon(b, hole = FALSE)), "1")))
## add buffer around lines and cut boundary polygon
cat("Converting contour lines to polygons...\n")
bcl <- gBuffer(cl, width = fd*diff(bbox(r)[1,])/3600000) # add small buffer so it cuts bounding poly
cp <- gDifference(b, bcl)
## restructure and make polygon number the ID
polys <- list()
for(j in seq_along(cp@polygons[[1]]@Polygons)) {
polys[[j]] <- Polygons(list(cp@polygons[[1]]@Polygons[[j]]),j)
}
cp <- SpatialPolygons(polys)
cp <- SpatialPolygonsDataFrame(cp, data.frame(id=seq_along(cp)))
# group by elev (replicate ids)
# ids = sapply(slot(cl, "lines"), slot, "ID")
# lens = sapply(1:length(cl), function(i) length(cl[i,]@lines[[1]]@Lines))
## cut the raster by levels
rc <- cut(r, breaks=plevels)
## loop through each polygon, create internal buffer, select points and define overlap with raster
cat("Adding attributes to polygons...\n")
l <- character(length(cp))
for(j in seq_along(cp)) {
p <- cp[cp$id==j,]
bp <- gBuffer(p, width = -max(res(r))) # use a negative buffer to obtain internal points
if(!is.null(bp)) {
xy <- SpatialPoints(coordinates(bp@polygons[[1]]@Polygons[[1]]))[1]
l[j] <- llevels[raster::extract(rc,xy)]
}
else {
xy <- coordinates(gCentroid(p)) # buffer will not be calculated for smaller polygons, so grab centroid
l[j] <- llevels[raster::extract(rc,xy)]
}
}
## assign level to each polygon
cp$level <- factor(l, levels=llevels)
cp$min <- plevels[-length(plevels)][cp$level]
cp$max <- plevels[-1][cp$level]
cp <- cp[!is.na(cp$level),] # discard small polygons that did not capture a raster point
df <- unique(cp@data[,c("level","min","max")]) # to be used after holes are defined
df <- df[order(df$min),]
row.names(df) <- df$level
llevels <- df$level
## define depressions in higher levels (ie holes)
cat("Defining holes...\n")
spolys <- list()
p <- cp[cp$level==llevels[1],] # add deepest layer
p <- gUnaryUnion(p)
spolys[[1]] <- Polygons(p@polygons[[1]]@Polygons, ID=llevels[1])
for(i in seq(length(llevels)-1)) {
p1 <- cp[cp$level==llevels[i+1],] # upper layer
p2 <- cp[cp$level==llevels[i],] # lower layer
x <- numeric(length(p2)) # grab one point from each of the deeper polygons
y <- numeric(length(p2))
id <- numeric(length(p2))
for(j in seq_along(p2)) {
xy <- coordinates(p2@polygons[[j]]@Polygons[[1]])[1,]
x[j] <- xy[1]; y[j] <- xy[2]
id[j] <- as.numeric(p2@polygons[[j]]@ID)
}
xy <- SpatialPointsDataFrame(cbind(x,y), data.frame(id=id))
holes <- over(xy, p1)$id
holes <- xy$id[which(!is.na(holes))]
if(length(holes)>0) {
p2 <- p2[p2$id %in% holes,] # keep the polygons over the shallower polygon
p1 <- gUnaryUnion(p1) # simplify each group of polygons
p2 <- gUnaryUnion(p2)
p <- gDifference(p1, p2) # cut holes in p1
} else { p <- gUnaryUnion(p1) }
spolys[[i+1]] <- Polygons(p@polygons[[1]]@Polygons, ID=llevels[i+1]) # add level
}
cp <- SpatialPolygons(spolys, pO=seq_along(llevels), proj4string=CRS(proj4string(r))) # compile into final object
## make polygons exclusive (i.e. no overlapping)
cpx = gDifference(cp[1,], cp[2,], id=cp[1,]@polygons[[1]]@ID)
for(i in 2:(length(cp)-1)) cpx = spRbind(cpx, gDifference(cp[i,], cp[i+1,], id=cp[i,]@polygons[[1]]@ID))
cp = spRbind(cpx, cp[length(cp),])
## it's a wrap
cp <- SpatialPolygonsDataFrame(cp, df)
cat("Done!")
cp
}
