有没有方法可以将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)
raster包（调用rgeos和sp）的解决方案。输出是一个SpatialPolygonsDataFrame，将覆盖网格中的每个值：library('raster')
rr <- raster(t(volcano))
rc <- cut(rr, breaks= 10)
pols <- rasterToPolygons(rc, dissolve=T)
spplot(pols)
这里的讨论会告诉您如何简化（美化）生成的多边形。
raster包的功能。干杯。 - Marc in the box受到这个网站的启发，我创建了一个将等高线转换为填充等高线的函数。它被设置成处理栅格对象并返回一个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()
plot(cp[4,], add=T, col='black') - 这会导致分析和可视化方面的问题（需要特定的绘图顺序，排除透明度等）。在这里，gDifference 可能会对你有所帮助。 - geotheoryfd，用于fairy dust，以解决我在使用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
}
contour(volcano, add=TRUE)已经解决了一部分问题，不是吗？ - thijs van den berghCOLS中的值。为了调试，请尝试一次只绘制一个cont[[k]]，以查看是否有任何问题。 - Carl Witthoftfilled.contour图。有几个问题：1.最低的等高线应该包括到角落的区域 - 我猜测你是对的，这可能是正确编码级别的问题。2.陨石坑的中心应该有一个150的水平面，但它被160覆盖了。 - Marc in the box