4

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 doesn't work on all sp objects, unfortunately. can you flesh them out more manually? – mdsumner Jul 15 '15 at 11:10
  • @mdsumner, I've edited the example code so it can now be constructed straight off, rather than using dput() – Andrew Jul 15 '15 at 14:08
  • A radial line sweep--in which the events include the vertices of the polygon and the other "set of points"--is the best method for a one-off solution. (That is, when you won't repeatedly be doing these calculations for a stream of points relative to a fixed central point and fixed polygon.) It requires only a simple modification of the algorithm (and code) at gis.stackexchange.com/a/33449. If you need highly accurate distances, first project everything using an azimuthal equidistant projection originating at the central point. – whuber Jul 15 '15 at 14:30
  • @whuber thanks for your suggestion. Unfortunately, my knowledge of R code and geometry was not good enough to cross-apply it, and I did need to be able to repeat the procedure, but I managed to find a solution (see below). Cheers! – Andrew Jul 17 '15 at 18:53
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

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.