Come si trasformano le curve di livello in contorni riempiti?

Question

Qualcuno sa di un modo per trasformare l'uscita dei poligoni contourLines per tracciare come contorni riempiti, come con filled.contours. Esiste un ordine su come devono essere tracciati i poligoni per poter vedere tutti i livelli disponibili? Ecco un esempio frammento di codice che non funziona:

#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

Una soluzione che utilizza il pacchetto raster (che chiama rgeos e sp). L'uscita è un SpatialPolygonsDataFrame che coprirà ogni valore in griglia:

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

Here's a discussion che vi mostrerà come semplificare ('abbellire') i poligoni risultanti.

Answer 2

Grazie a qualche ispirazione da this sito, ho lavorato una funzione per convertire le linee di contorno ai contorni pieni. È configurato per elaborare un oggetto raster e restituire 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 

} 

Probabilmente detiene diversi inefficienze, ma ha funzionato bene nei test che ho condotto utilizzando i dati batimetrici. Ecco un esempio utilizzando i dati del vulcano:

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

Costruire per l'ottimo lavoro di Paolo regolare, qui è una versione che dovrebbe garantire poligoni esclusivi (cioè senza sovrapposizione).

Ho una nuova discussione fd per Fairy Dust per risolvere un problema che ho scoperto a lavorare con le coordinate UTM-type. Fondamentalmente, come ho capito, l'algoritmo funziona campionando i punti laterali dalle linee di livello per determinare quale lato è all'interno del poligono. La distanza del punto di campionamento dalla linea può creare problemi se finisce, ad es. dietro un altro contorno. Quindi, se i vostri poligoni risultante è prova sbagliato impostazione fd ai valori 10^± n fino a che non sembra molto sbagliato o circa la destra ..

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 
}