5

I am solving the same problem as How to split a polygon into quadrants?, except that I am looking for a solution in R or QGis. So I have a polygon layer shapefile, which is unfortunatelly not exactly rectangular (due to some reprojections etc.). I would like to split the "squares" into 2x2, 4x4 or 8x8 smaller squares, as equal as possible.

enter image description here

Note that I am not looking for tool to generate a new grid, but to split the one I already have.

6
  • In QGIS, I don't think there's a built-in function to "split" polygons into your desired dimensions. If via the interface then you would have to create a new grid, clip your polygon layer onto it and then use Join by location. I'm fairly certain it's possible to do this from R or writing a python script in QGIS if you know your coding but hopefully others can advise.
    – Joseph
    Dec 5, 2014 at 11:12
  • Thanks @Joseph, I don't think this is a good solution, the square net is little bit irregular so each polygon should be divided separately so that we get it split into sub-squares as-equal-as-possible. I will wait for other advices.
    – Tomas
    Dec 5, 2014 at 11:59
  • No problem @Thomas, a final question. Couldn't you Dissolve the polygon layer first to remove all "squares" and then do those steps mentioned? Or must the "squares" remain at their exact location?
    – Joseph
    Dec 5, 2014 at 12:08
  • @Joseph please see the last sentence of my question.
    – Tomas
    Dec 5, 2014 at 12:10
  • Ahh...don't know how I missed that!
    – Joseph
    Dec 5, 2014 at 12:12

1 Answer 1

2

Here I address the question in detail. Here I cover the problem in detail. In resume I propose two ways.

Option 1 - simple way

If you dont need the lines to be coincident among polygons what you can do is this:

library(sf)
library(dplyr)

layer <- read_sf('./data/20220923/A.gpkg')
layer$id <- rownames(layer)
# mapview(layer)
# set number of divisions
div <- 4
ls <- list()
for (i in 1:nrow(layer)){
    x <- st_make_grid(layer[i,], n = div) %>% st_as_sf() # divide pol
    xx <- st_intersection(layer[i,], x) # intersection
    xx$id <- layer[i,]$id 
    ls[[i]] <- xx
}

l2 <- sf::st_as_sf(data.table::rbindlist(ls)) # superfast
# mapview(l2)+layer

Layer looks like this: enter image description here

And then you have it with a 4x4 subdivision example: enter image description here

Option 2 - complex way (better results)

If you need the subdivision lines to be coincident you can follow the process in the link with the steps explained. The function has a paramenter for the number of divisions `div''. Here is the code in case you just what to try it:

# This functions is made to short coords clock wise (implemented inside QuasiRectangularSpliter)
sortcoords <- function(x){
    # sort points to construct lines clock wise lines
    # 1. Upper left corner: from the 2 leftest points, take the upper one
    # 2. Upper right corner: ...
    # 3. Bottom right corner: ...
    # 4. Bottom left corner: ...
    # 5. Close coord = 1
    x <- unique(x)
    x1 <- ((as.data.frame(x) %>% arrange(X))[1:2,] %>% arrange(desc(Y)))[1,]
    x2 <- ((as.data.frame(x) %>% arrange(desc(X)))[1:2,] %>% arrange(desc(Y)))[1,]
    x3 <- ((as.data.frame(x) %>% arrange(desc(X)))[1:2,] %>% arrange(Y))[1,]
    x4 <- ((as.data.frame(x) %>% arrange(X))[1:2,] %>% arrange(Y))[1,]
    x5 <- x1
    do.call(rbind, list(x1, x2, x3, x4, x5)) %>% as.matrix()
}

QuasiRectangularSpliter <- function(layer, div = 4, adjust = 0.5){
    subpolygons <- list()
    for (poly in 1:nrow(layer)){
        # poly <- 2 # funcional
        # poly <- 9 # erratico
        # cat('Creating subpolygons for ', poly, ' of ', nrow(layer), '\n')
        # Create coords
        pol <- layer[poly,]
        xx <- st_cast(pol, "LINESTRING")
        xx <- st_node(xx)

        x <- st_cast(xx, "POINT")
        x <- x %>% st_coordinates()
        x <- sortcoords(x)
        
        # Create lines
        ls <- list()
        for (i in 1:4){
            # i <- 1
            x1 <- x[i, 1]
            x2 <- x[i+1, 1]
            y1 <- x[i, 2]
            y2 <- x[i+1, 2]
            lcoord <- matrix(c(x1,x2,y1,y2), 2)
            l <- st_linestring(lcoord) %>% st_sfc() %>% st_as_sf()
            l$id <- pol$id
            ls[[i]] <- l
        }
        
        l2 <- sf::st_as_sf(data.table::rbindlist(ls)) # superfast
        st_crs(l2) <- st_crs(pol)
        # mapview(l2, color = "red")+layer
        
        # Create point for lines
        lsample <- seq(0, 1, by=1/div) # seq for st_line_sample
        
        ps <- list()
        for(i in 1:4){
            p <- st_line_sample(l2[i,], type = "regular", sample = lsample) %>%
                st_cast("POINT") %>% st_sfc() %>% st_as_sf()
            
            # CORRECTION to be sure about creating totally overlaping lines:
            if(i == 1){ # points in horizontal upper line --> y+1
                st_geometry(p) <- st_geometry(p) + c(0, adjust)
            } else if (i == 2){ # points in vertical right line --> x+1
                st_geometry(p) <- st_geometry(p) + c(adjust, 0)
            } else if (i == 3){ # points in bottom horizontal line --> y-1
                st_geometry(p) <- st_geometry(p) + c(0, -adjust)
            } else { # points in left vertical line --> x-1
                st_geometry(p) <- st_geometry(p) + c(-adjust, 0)
            }
            
            # return
            ps[[i]] <- p
        }
        
        # Create vertical lines
        vlines <- list()
        for(i in 1:(div+1)){
            p1 <- ps[[1]][i,]
            p2 <- ps[[3]][div+2-i, ]
            lcoord <- rbind(st_coordinates(p1), st_coordinates(p2))
            l <- st_linestring(lcoord) %>% st_sfc() %>% st_as_sf()
            l$id <- pol$id
            vlines[[i]] <- l
        }
        
        # Create horizontal lines
        hlines <- list()
        for(i in 1:(div+1)){
            p1 <- ps[[4]][div+2-i,]
            p2 <- ps[[2]][i, ]
            lcoord <- rbind(st_coordinates(p1), st_coordinates(p2))
            l <- st_linestring(lcoord) %>% st_sfc() %>% st_as_sf()
            l$id <- pol$id
            hlines[[i]] <- l
        }
        
        # merge layers
        v2 <- sf::st_as_sf(data.table::rbindlist(vlines)) # superfast
        h2 <- sf::st_as_sf(data.table::rbindlist(hlines)) # superfast
        st_crs(v2) <- st_crs(pol)
        st_crs(h2) <- st_crs(pol)
        mapview(v2)+h2+pol
        
        lgrid <- rbind(v2, h2)
        polygons <- lgrid %>% st_cast("MULTILINESTRING")
        polygons = st_collection_extract(st_polygonize(st_union(polygons))) %>%
            st_as_sf()
        polygons$id = pol$id
        polygons <- st_cast(polygons, "MULTIPOLYGON")
        
        # add polygons to list
        subpolygons[[poly]] <- polygons
        
        # mapview(polygons)+v2+h2+pol
    }
    # merge layers
    mergedpols <- sf::st_as_sf(data.table::rbindlist(subpolygons)) # superfast

    # return
    mergedpols
}

# let's try it:
l4 <- QuasiRectangularSpliter(layer, 4)
mapview(l4)

enter image description here

As you can see here, lines are following the not exactly rectangular shape and are connected among polygons. You can iterate over it or tweak it to make other levels of subdivisions.

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.