如何将轮廓线变成填充轮廓？

Question

有谁知道将contourLines 多边形的输出转换为填充轮廓的方法，就像filled.contours 一样。是否有顺序必须如何绘制多边形才能查看所有可用级别？这是一个不起作用的示例 sn-p 代码：

#typical plot
filled.contour(volcano, color.palette = terrain.colors)

#try
cont <- contourLines(volcano)
fun <- function(x) x$level
LEVS <- sort(unique(unlist(lapply(cont, fun))))
COLS <- terrain.colors(length(LEVS))
contour(volcano)
for(i in seq(cont)){
    COLNUM <- match(cont[[i]]$level, LEVS)
    polygon(cont[[i]], col=COLS[COLNUM], border="NA")
}
contour(volcano, add=TRUE)

Answer 1

使用raster 包（调用rgeos 和sp）的解决方案。输出是一个SpatialPolygonsDataFrame，它将覆盖网格中的每个值：

library('raster')
rr <- raster(t(volcano))
rc <- cut(rr, breaks= 10)
pols <- rasterToPolygons(rc, dissolve=T)
spplot(pols)

Here's a discussion 将向您展示如何简化（“美化”）生成的多边形。

Answer 2

感谢this 网站的一些启发，我设计了一个将等高线转换为填充等高线的功能。它设置为处理栅格对象并返回 SpatialPolygonsDataFrame。

raster2contourPolys <- function(r, levels = NULL) {

  ## set-up levels
  levels <- sort(levels)
  plevels <- c(min(values(r), na.rm=TRUE), levels, max(values(r), 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")
  cl <- contourLines(rx,ry,rz,levels=levels) 
  cl <- ContourLines2SLDF(cl)

  ## 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 = 0.0001) # 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)))

  ## 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[extract(rc,xy)]
    } 
    else { 
      xy <- coordinates(gCentroid(p)) # buffer will not be calculated for smaller polygons, so grab centroid
      l[j] <- llevels[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
  cp <- SpatialPolygonsDataFrame(cp, df)
  cat("Done!")
  cp

}

它可能存在一些效率低下的问题，但它在我使用测深数据进行的测试中运行良好。下面是一个使用火山数据的例子：

r <- raster(t(volcano))
l <- seq(100,200,by=10)
cp <- raster2contourPolys(r, levels=l)
cols <- terrain.colors(length(cp))
plot(cp, col=cols, border=cols, axes=TRUE, xaxs="i", yaxs="i")
contour(r, levels=l, add=TRUE)
box()

Answer 3

在 Paul Regular 的出色工作的基础上，这里有一个版本应该确保独占多边形（即没有重叠）。

我为 fairy dust 添加了一个新参数 fd，以解决我发现使用 UTM 类型坐标的问题。基本上，据我了解，该算法通过从轮廓线中采样横向点来确定多边形内部的哪一侧。如果样本点与线的距离最终出现在例如在另一个轮廓后面。因此，如果您生成的多边形看起来错误，请尝试将 fd 设置为值 10^±n，直到它看起来非常错误或几乎正确..

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
}