Least Squares in Raster grid (spa)

En este caso, alguien busca como encajar cuadros de diversos tamaños en un raster grid (sus motivos tendrá), de forma que todo el raster con valores tenga esté incluido en algún tipo de cuadro. Deben estar distribuidos en orden de mayor a menor…

La solución se propone para una cuestión vista vista en GIS STACKEXCHANGE: Packing squares polygons in a raster grid{target:“_blank”}

En este caso, no se busca una solución ideal con el menor número de polígonos posibles, si no una solución sistemática. El proceso es relativamente sencillo, se busca encajar los cuadros de mayor tamaño empezando la búsqueda desde la esquina inferior izquierda, de izquierda a derecha y de abajo arriba. Una vez el cuadro tiene espacio para emplazarse, el raster se actualiza y se continúa la búsqueda. Cuando no entran más cuadros grandes se pasa al siguiente tamaño y se repite el proceso.

library(raster)
# Note: ncell start from upper left corner
# Note: polygon start: bottom left corner...check raster val in that order

# Parameters...
dim <- 10
holes = 10

# Create raster with random holes
r <- raster(nrow=dim, ncol=dim, xmn=0, xmx=dim, ymn=0, ymx=dim, crs=NA)
values(r) <- 1
set.seed(123)
nas <- sample(ncell(r), holes)
r[nas] <- NA; r[!is.na(r)] <- 1
par(mfrow=c(1,1))
plot(r)

R-plot

# Set new raster based on original raster (this will be modified during iterations)
nr <- r

# Try sqares backwards (from larger to smaller)
pols <- list() # Store polygons
for(i in dim:1){
    # Set coordinates to explore
    xs <- c(0:(dim))
    ys <- c(1:(dim+1))
    for(iy in seq_along(ys)){
        for (ix in seq_along(xs)){
            cell <- cellFromXY(nr, c(xs[[ix]], ys[[iy]])) # Starting cell
            
            # Check if polygon exceds the raster dims
            polxmax <- xs[[ix]]+i
            polymax <- ys[[iy]]+i
            
            if(polxmax > extent(nr)[2] | polymax > extent(nr)[4]+1){
                cat("Polygon extent exceds the raster size... \n")
                break("Polygon extent exceds the raster size...\n")
            } else if (is.na(nr[cell])){
                cat("Starting cell is NA --> Skip to next starting point \n")
            } else{
                cat("Iteration dim=", i, " | iy=", iy, " | ix=", ix, "\n")
                
                # create moving polygon
                rc <- raster(nrow=i, ncol=i,
                             xmn=xs[[ix]], xmx=xs[[ix]]+i,
                             ymn=ys[[iy]]-1, ymx=ys[[iy]]+i-1, crs=NA)
                pol <- rasterToPolygons(rc, dissolve = TRUE)
                
                ext <- raster::extract(nr, pol) # Extract
                if(NA %in% ext[[1]]){
                    cat("NA in extraction. Raster res ", i, "\n")
                }else{
                    cat("Store pol in pols list \n")
                    pols <- append(pols, list(pol))
                    
                    # update raster (cells to NA)
                    cat("Update raster cells to NAs \n")
                    polr <- rasterize(pol, nr)
                    polr[is.na(polr)] <- 0
                    xraster <- nr+polr
                    nr[xraster==2] <- NA
                    
                    # plot(nr); plot(pol, add = TRUE) #uncoment to show develpment
                }
            }
        }
    }
}

# merge polygons
allpols <- do.call(rbind, pols)

La ejecución del proceso ejecuta también plots para ver el comportamiento de la selección de cuadros.

Para ver el resultado comparativo entre el raster original:

par(mfrow=c(1,2))
plot(r, col = rgb(red = 0, green = 1, blue = 0, alpha = 0.3))
plot(r, col = rgb(red = 0, green = 1, blue = 0, alpha = 0.3))
plot(allpols, col = rgb(red = 1, green = 0, blue = 0, alpha = 0.1), add=TRUE)

R-plot