Splitting a Polygon into Quadrants
This is the QUESTION IN GIS STACK EXCHANGE that it is not totally solved in our favorite forum. The aim is to create a sub-grid of 2x2, 4x4, etc within each polygon. This can be used recursively in case more sub-grids are needed. The goal is not to create a new grid but to adapt a sub-grid to the boundaries of the original polygons.
The process will be addressed in two different ways:
- The Option 1: is a simple one and easier to understand in case you are in a hurry and you don’t need totally aligned lines among polygons, but just a sub-grid within each one.
- The Option 2: is more complex but it creates the sub-grid adapting the lines to the boundaries shift.
The process is just valid for a layer composed by quadrilateral polygons (not necesarily rectangular) as you can see in the main image.
Option 1: quite simple process
The following layer has been created based in the specifications from the question –> “not exactly rectangular”. It was made in QGIS and looks like this:
library(mapview)
library(sf)
library(dplyr)
library(htmlwidgets)
library(rgeos)
layer <- read_sf('../../data/20220923/A.gpkg')
layer$id <- rownames(layer)
One possible solution is to make a grid within each of the polygons. If you apply it directly to the whole extension of polygons the results could not be as good as you are looking for. So do it individually instead:
# 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
x$id <- layer[i,]$id
ls[[i]] <- x
}
l2 <- sf::st_as_sf(data.table::rbindlist(ls)) # superfast
This tools helps you to get the most identical parts within the boundaries, but as you can see, the boundaries of each polygon present not a good result. What you can do to get the perfect internal boundaries is to add other step in the previous function: an intersection between the original polygon and the internal grid. You have to notice that the grid tool result is going to contain totally the original polygon.
Check it with the applied correction:
xx <- st_intersection(layer[16,], l2)
So applying this to the whole function we have:
# set number of divisions
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
If you dont care about the unconnected lines among polygons, this can be fine. If you are looking for a perfect result we encourage you to get deeper in geometry tools in the following section.
Option 2: hard process
As we have checked it out in the last process, we found a solution but it looks no so perfect as it should, because polygon nodes are in different places in each of the polygons. The hard process is done creating the nodes previously by dividing the periferal lines from the polygons and then creating the grid from those lines.
We are goint to perform this polygon by polygon in order to get the most equalized subpolygons within it. First select one polygon to develop the algorithm. Later on, it will be applyed to the whole layer.
pol <- layer[1,]
# Create coords
xx <- st_cast(pol, "LINESTRING")
xx <- st_node(xx)
x <- st_cast(xx, "POINT")
x <- x %>% st_coordinates()
# 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)
Now that we have the separated lines we have to get the internal points and match them to created the “vertical” and “horizontal” lines.
div <- 8
lsample <- seq(0, 1, by=1/div) # seq for st_line_sample
adjust <- 1 # modifier of internal coords to get totally overlaping lines later
ps <- list()
for(i in 1:4){
# i <- 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
}
Points are sampled over lines in a clockwise order. With this in mind,
we will create vertical and horizontal lines to intersect them later to
create the grid. Remember, lines will be separated by div + 1 points.
To create vertical lines, we will start from left to right. The first line will we coincident with the poligon edge one. The same will happen to the last one. So these lines are already “created”. These are the line 4 and the line 2. We important part is to create the internal ones. To do so, we need to create create the lines by paring the correct points from up and down. The first point from the upper line goes with the last point from the bottom line and so on. For the horizontal lines, the process is the simetrical one.
vlines <- list()
for(i in 1:(div+1)){
# i <- 2
if (i == 1){
vlines[[i]] <- ls[[4]]
} else if (i == div+1){
vlines[[i]] <- ls[[2]]
} else{
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
}
}
hlines <- list()
for(i in 1:(div+1)){
# i <- 2
if (i == 1){
hlines[[i]] <- ls [[1]]
} else if (i == div+1){
hlines[[i]] <- ls[[3]]
} else{
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
}
}
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
You can test it with other divisions.
Now, we have to create a grid with all the intersections of lines:
lgrid <- rbind(v2, h2)
# cast to polygons
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")
Now that we have develop the algorith for each polygon we can apply it to the whole layer. Take into account that all the subpolygons created are identified with the id of the main polygon, but they are not sorted in any way. You could do that by another algorith, for examplo sorting them from left to right and from up to down and then assigning a new feature id for each square.
Extra step This function has been create to create lines ordered
clock wise. This is needed as polygons are not always constructed in the
same order by their defined coords. This is implemented in the
QuasiRectangularSpliter function show below.
# This has been implemented later
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()
}
Here is the simplyfied function to apply to your layer:
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
}
Split the original grid into 4x4:
l4 <- QuasiRectangularSpliter(layer, 4)
# mapview(l4)
And 8x8:
l8 <- QuasiRectangularSpliter(layer, 8)
# mapview(l8)
Hope it helps¡
