I am trying to calculate the distances between a central point, and the edge of a polygon, along straight lines drawn between the central point and another set of points (where necessary extended beyond the outer point to the polygon).
My data look like this:
library(sp)
library(raster)
# Centre point
Centre <- data.frame(x= 38.35419, y= 4.483533)
Centre <- SpatialPointsDataFrame(coords=Centre, data=Centre)
projection(Centre) <- CRS('+proj=longlat')
crs(Centre) <- "+proj=longlat +datum=WGS84"
# Sites
Site <- c(seq(1,10,1))
x <- c(38.29013,37.90987,38.29999,38.66006,38.15978,38.25978,38.15997,38.24062,38.25005,38.00712)
y <- c(5.33006,4.96017,5.23995,4.74010,5.09011,5.09015,5.06016,4.99723,5.06002,4.94634)
Sites <- data.frame(Site, x, y)
xy <- Sites[2:3]
Sites2 <- SpatialPointsDataFrame(coords=xy, data=Sites)
projection(Sites2) <- CRS('+proj=longlat')
crs(Sites2) <- "+proj=longlat +datum=WGS84"
# Polygon
Sr1 = Polygon(cbind(c(37.95169, 37.95169, 38.65000, 38.65000),
c(5.17153, 3.85601, 3.85601, 5.17153)))
Srs1 = Polygons(list(Sr1), "s1")
Hull = SpatialPolygons(list(Srs1), proj4string=CRS("+proj=longlat +datum=WGS84"))
I have managed to calculate the distance from the central point to each outlying point using:
spDistsN1(Sites2, Centre, longlat = TRUE)
Which returns the distance in kilometres.
To visualise the data, I use the code:
library(rgeos)
library(maptools)
data(wrld_simpl)
plot(wrld_simpl, xlim=c(37.5,39.5), ylim=c(4,5), axes=T)
plot(Hull, col='green', axes=T, add=T)
points(Sites2$x, Sites2$y, pch=20)
plot(Centre, add=T, col="yellow", pch="+", cex=1.5)
and have managed to draw lines from the central point to each site:
C_P_lines <- list()
for (i in 1:length(Sites$x)) {
C_P_lines[i] <- Lines(list(
Line(rbind(Centre@data, Sites[i,2:3])))
, ID=Site.coords[i,1])
}
C_P_spLines <- SpatialLines(C_P_lines, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
plot(C_P_spLines, add=T)
What I need to do, effectively, is to a) extend these lines to the edge of the polygon (for sites within the polygon), and then b) extract the intersection points of the lines with the polygon.
a) I have so far been unsuccessful in finding a solution to extending the lines.
b) I need to achieve the results shown here: http://uk.mathworks.com/help/map/ref/polyxpoly.html but in R. I have tried using gIntersection
in the rgeos
package,
Edge_points <- gIntersection(C_P_spLines, Hull, byid=T)
plot(Edge_points, col="red", add=T)
but this just crops the lines to the extent of the polygon.
I'm really new to using R for GIS work.
dput()
– Andrew Jul 15 '15 at 14:08