# 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. 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? 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/…) 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! Aug 3 '13 at 19:46
• When using sf features, the tool st_make_grid creates polygons of equal size that covers a shapefile. 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)
[1] 4
> min(HL)
[1] 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. 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! 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. 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")