怎样才能把轮廓线到充满轮廓？(How does one turn contour lines int

有谁知道一个办法把输出contourLines多边形，以绘制作为填充的轮廓，与filled.contours 。 是否有怎样的多边形必须然后才能看到所有可用的水平绘制的订单？ 下面是代码不起作用的例子片断：

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

这里有一个讨论会告诉你如何简化（“美化”）所产生的多边形。

Answer 2:

由于从一些启示这个网站，我的工作了一个功能轮廓线转换为填充的轮廓。 它的建立来处理光栅对象，并返回一个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:

保罗定期的优秀工作的基础上，这里是应该确保专属的多边形（即没有重叠）的一个版本。

我添加了一个新的说法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
}

文章来源: How does one turn contour lines into filled contours?