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
dput
andstructure
as explained in this post? And what is your ultimate aim - are you mapping something? – SlowLearner Aug 3 '13 at 8:42sf
features, the toolst_make_grid
creates polygons of equal size that covers a shapefile. – user3386170 Jan 16 '19 at 16:42