How to find the distance from a point to a polygon, via another point, in R

Question

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.

Answer 1

Others users might be interested to know that I managed to solve this.

In order to extend the lines outside of the polygon (problem a):

# Save centre point co-ordinates separately
A.x <- Centre@coords[1]
A.y <- Centre@coords[2]

Using the great answer posted here: https://stackoverflow.com/a/7741655/4312169 I applied the code:

# Generate points outside polygon along line of `C_P_spLines` (I played around with the multiplier at the end to ensure that all the points ended up outside of the polygon, but not too far out)
Out.site.coords <- Sites[,1:3]
for (i in 1:length(Sites$x)) {
Out.site.coords$x[i] <- Sites[i,2] + (Sites[i,2] - A.x) / sqrt((A.x - Sites[i,2])^2 + (A.y - Sites[i,3])^2) * 1
Out.site.coords$y[i] <- Sites[i,3] + (Sites[i,3] - A.y) / sqrt((A.x - Sites[i,2])^2 + (A.y - Sites[i,3])^2) * 1
}

# Make it a SpatialPointsDataFrame
xy <- Out.site.coords[2:3]
Out.sites <- SpatialPointsDataFrame(coords=xy, data=Out.site.coords)
projection(Out.sites) <- CRS('+proj=longlat')
crs(Out.sites) <- "+proj=longlat +datum=WGS84"
class(Out.sites)
Out.sites

# Re-plot with wider extent (based on min and max co-ords from "Out.sites")
plot(wrld_simpl, xlim=c(37,40.5), ylim=c(2.5,6.6), axes=T)
plot(Hull, col="green", add=T)
points(Sites$x, Sites$y, pch=20)
plot(Centre, add=T, col="yellow", pch="+", cex=1.5)
points(Out.sites$x, Out.sites$y, pch=20, col="blue")

# Create SpatialLines object for "centre to out.sites"
Radial_Lines <- list()
for (i in 1:length(Out.site.coords$x)) {
  Radial_Lines[i] <- Lines(list(
    Line(rbind(Centre@coords, Out.site.coords[i,2:3])))
    , ID=Out.site.coords[i,1])
}

# Make it a SpatialLines object
Radial_spLines <- SpatialLines(Radial_Lines, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

I then used gIntersection and some subsetting to extract the outer points from the lines - see answer here: https://gis.stackexchange.com/a/154864/55218