# How to split a polygon using R

I have an irregular polygon which is a spatial object. I want to divide it into multiple polygons (say, n=25) with equal height. To show how I want to divide it, I've plotted grids on the background using:

``````grid(NA,25)
``````

If it is not possible, can anyone at least tell me how I'd be able to get the coordinates of the intersection points between the gridlines (which is not a spatial object, however) and the polygon's boundary?

• Can't give a full answer right now, but look at the rgeos package for polygon/polygon overlay/intersection etc operations. – Spacedman Aug 1 '13 at 6:56
• Can you share the polygon data either as a shapefile to download or if it is not very large by using `dput` and `structure` as explained in this post? And what is your ultimate aim - are you mapping something? – SlowLearner Aug 3 '13 at 8:42
• well ... I figured it out ... first I had to make a spatial line of my polygon. Then I produced a sequence of straight lines instead of grids. As Spaceman mentioned, I used 'rgeos' package->gIntersection to get the intersecting points between those spatial lines (here: docs.google.com/file/d/0B6GUNg-8d30vQnUxM3JOcGRaVkk/…) – ToNoY Aug 3 '13 at 15:35
• Can you answer your own question with your answer, including the R code you used? In this way it becomes a resource for other people who have the same problem - and that's why this site exists! – SlowLearner Aug 3 '13 at 19:46
• When using `sf` features, the tool `st_make_grid` creates polygons of equal size that covers a shapefile. – user3386170 Jan 16 '19 at 16:42

``````HL <- matrix(c(1,1,2,3,3,2,1, 1,2,4,2,1,0,1), ncol=2, byrow=FALSE)
``````

At first, you have to make a spatial line with these points

``````require(rgdal)
require(sp)
L = Line(HL)
Ls = Lines(list(L), ID = "a")
SL = SpatialLines(list(Ls))
proj4string(SL) = CRS("+proj=utm +zone=46 +datum=WGS84 +units=m +no_defs")
plot(SL)
``````

Know the max and min extent of your data

``````> max(HL)
 4
> min(HL)
 0
``````

Now create vertical and horizontal sequence from the extent with an interval of say, 1 (depends on how you want your grids to be spaced)

``````vp<-seq(0, 4, by= 1)
hp<-seq(0, 4, by= 1)
``````

To make horizontal grid like lines, first create points by combining Xmin and Xmax with all Y values

``````vpXmin<-cbind(rep(min(hp),length(vp)), vp)
vpXmax<-cbind(rep(max(hp),length(vp)), vp)
``````

Now create a loop to combine each row of vpXmin and vpXmax to get the coordinates for each line

``````l <- lapply(1:nrow(vpXmin),function(i) rbind(vpXmax[i,],vpXmin[i,]))
``````

``````L1<-lapply(1:length(l),function(i) Line(l[[i]]))
Ls1<-lapply(1:length(L1),function(i) Lines(list(L1[[i]]), ID='a'))
SL1<-lapply(1:length(Ls1),function(i) SpatialLines(list(Ls1[[i]])))
for (i in 1:length(SL1)) {proj4string(SL1[[i]]) = CRS("+proj=utm +zone=46 +datum=WGS84 +units=m +no_defs")}
``````

Now unlist the spatial lines

``````#  Get the Lines objects which contain multiple 'lines'
ll0 <- lapply( SL1 , function(x) `@`(x , "lines") )
#  Extract the individual 'lines'
ll1 <- lapply( unlist( ll0 ) , function(y) `@`(y,"Lines") )
#  Combine them into a single SpatialLines object and give projection definition
Sll <- SpatialLines( list( Lines( unlist( ll1 ) , ID = "a" ) ),
proj4string=CRS("+proj=utm +zone=46 +datum=WGS84 +units=m +no_defs" ))
``````

Once you have both your grids and the polygon as spatial line objects, getting the intersection points between those is easy

``````require(rgeos)
intrsct <-gIntersection(Sll, SL)
#plot the polygon line
plot(SL, axes=TRUE)
#plot grid lines
#plot intersection points
#the intersection points
View(intrsct)
``````

• Good start! However, there remain some significant problems with this approach. First, having the intersection points does not suffice to split the original polygon. Second, the intersections may contain entire line segments. Third, you have hard-coded some polygon-specific parameters. To see the effects of these problems, apply your code to `HL <- matrix(c(0,0,1,3,2,3,0, 0,4,2,4,2,0,0), ncol=2, byrow=FALSE)`. You will need to set `vp<-seq(0, 4, by=1); hp<-seq(0, 3, by=1)` to obtain non-empty intersections. – whuber Aug 6 '13 at 14:32
• Thanks for the feedback, I know the solution is not universal, but it seems to work in many cases! – ToNoY Aug 7 '13 at 14:40
• Yes, such cases even have names: the polygon can be split into two sides (left and right) each of which is (strictly) monotone with respect to the y-axis. (That means it doesn't ever reverse itself on its path from bottom to top.) Many computational geometry algorithms try to reduce a general situation to one in which monotone curves are processed. All convex polygons have two monotone sides--your procedure works best for them--but there are non-convex polygons with monotone sides, too. – whuber Aug 7 '13 at 14:47

When using sf features, the tool `st_make_grid` creates polygons of equal size that covers a shapefile. Set the spacing between each grid line with `cell size`.

Adapting an exemple from the R Documentation:

``````poly_sf<-st_polygon(list(rbind(c(0,0), c(0.5, 0.2),c(1,0), c(1,1), c(0,0))))
grid<-st_make_grid(poly_sf, cellsize = .1, square = TRUE)
plot(poly_sf, col="red")
If working with `sp` type features, you can convert between `sp` and `sf` with:
``````grid_sp<-as(grid, "Spatial)#converts sf grid polygons into sp features 