9

In the green area of the raster I have to fit as many square polygons as possible, of three different sizes:

  • 3x3 cell size
  • 2x2 cell size
  • 1x1 cell size

With the fact that the approach should be greedy, so first fitting the biggest polygons, then the middle sized, and then if there is any space left, the smallest size. It does not have to be an exact approach but at least an approximation.

Considering this is just a test area for the algorithm, which should eventually perform this search on bigger real case rasters, where basically we have objects of three different sizes and are trying to place them in the allowed zones of the raster.

Having said that, I am aware it's not at all an easy task. I am using Qgis with R scripting to try to formulate this search.

Raster image from R

EDIT: By moving a 3x3 window square and suming up the cell values, it can be determined where are the possible locations for the squares:

raster matrix

So, wherever there is a 9 means the square can fit, but the squares cannot overlap so now the question is how to find the best combination of 9 cells that do not overlap

9
  • Scan from top left to bottom right. When you find a 3x3 green spot, fill it. Continue. That gets you a number of 3x3 fills in linear time, but is it optimal? Can we prove that? Is there an optimal solution and what's the complexity of finding it? I don't think the complexity can be that bad since there's only a finite number of cells in the raster...
    – Spacedman
    Commented Aug 24, 2017 at 11:38
  • I now have code that does the above, ie scan and find space for squares. But I'm not sure it constitutes an answer unless I know what you really want is an optimum solution (if such exists...)
    – Spacedman
    Commented Aug 24, 2017 at 13:54
  • Is a square what we see that looks like pixels ? In that case, the square area seems to be of 200x200 units, no ? 3x3 suare area would make a great number of squares...
    – gisnside
    Commented Aug 24, 2017 at 15:58
  • @Spacedman Yes exactlly what I managed to do, by moving the 3x3 square and summing up the values of the cells. If the sum is 9, that means the square can fit. After fitting all possible 3x3, mark those locations as 9, repeat the process with 2x2, and same for 1x1. But depending on where you start the search, different solutions are possible and how to prove the optimality
    – Dragana
    Commented Aug 24, 2017 at 20:54
  • 1
    Please share your code! If it’s big, perhaps a link to a gist?
    – Simbamangu
    Commented Aug 29, 2017 at 5:51

1 Answer 1

3

I've created a simple algorithm that does that. The approach is the following:

  1. Create the larger polygon size and check cell by cell. If the polygon exceeds the extent of the raster, go to the next size. Sizes are checked backwards from larger to smaller ones.
  2. Go through cells from starting point (bottom left corner) and check if it fits within the raster (all values inside polygon must be not NA). If this happens the raster is updated so there will be no other polygon that overlaps.
  3. Set next polygon size and explore again to find gaps with values.

The process has still room to be optimized as it doesn't look for the least number of polygons solution. This could be achieved adding randomization to the polygon search or creating different starting points. That could complicate things...

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)

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

# Try squares 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 exceeds the raster dims
            polxmax <- xs[[ix]]+i
            polymax <- ys[[iy]]+i
            
            if(polxmax > extent(nr)[2] | polymax > extent(nr)[4]+1){
                cat("Polygon extent exceeds the raster size... \n")
                break("Polygon extent exceeds 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)
                }
            }
        }
    }
}

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

Show results:

# plot with r
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)

enter image description here enter image description here

20220328 - Update

In case you just want squares of size 3 as max, change the first for loop by a parameter:

...
squaremaxsize = 3

for(i in squaremaxsize:1){
...

That way you will find a new result...

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.