6

I have an irregular polygon which is a spatial object.

enter image description here

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?

|improve this question
  • 1
    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
  • 2
    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
  • 5
    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 at 16:42
2

Say your polygon is made up of the following points:

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,]))

To make spatial objects from your gridded lines, follow these steps

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(Sll, lty=2, add=TRUE)
#plot intersection points
plot(intrsct, add=TRUE, pch=16, col='red')
#the intersection points
View(intrsct)

Here's the final outcome: https://docs.google.com/file/d/0B6GUNg-8d30vX29wQ0gzTU5pZkE/edit?usp=sharing

|improve this answer
  • 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
1

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")
plot(grid, add=T)

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
grid_sf<-st_sf(grid)#convert sp back to sf

enter image description here

|improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.