# How to split a polygon into quadrants in R/Qgis?

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.

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

• 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. 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. 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? Dec 5, 2014 at 12:08
• @Joseph please see the last sentence of my question. Dec 5, 2014 at 12:10
• Ahh...don't know how I missed that! Dec 5, 2014 at 12:12

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\$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:

And then you have it with a 4x4 subdivision example:

# 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")

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)
``````

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.