# How to find the coordinates of the intersection points between two SpatialLines objects (or SpatialLines and SpatialPolygon edges)

I am try to find the co-ordinates of points where SpatialLines cross the edge of a SpatialPolygon. I have searched and searched, but can only find results such as `over` which tell me which lines cross but not where, or `gIntersection` which cut the lines down to sections within the polygon.

Here is a subset of my data:

``````# Required packages
library(rgeos)
library(sp)
library(maptools)
library(raster)
library(grDevices)

# Create convex hull
# Points to create hull
Hull_pts <- structure(list(x = c(38.87584, 38.89215, 38.87062, 38.86157, 38.72808, 38.22315, 38.12702, 38.05936, 37.95169, 37.96915, 38.14758, 38.22325, 38.34001, 38.6394, 38.6447),
y = c(4.17092, 4.05521, 3.93639, 3.9168, 3.85601, 3.89487, 4.17066, 4.38951, 4.98351, 4.99706, 5.12187, 5.17153, 5.16003, 4.77422, 4.76607)),
.Names = c("x", "y"),
row.names = c(1574L, 1540L, 1490L, 1482L, 1457L, 1473L, 294L, 1718L, 1131L, 974L, 2838L, 2101L, 111L, 1914L, 1909L),
class = "data.frame")

# Create SpatialLines of Hull
Hull_lines <- list()
for (i in 1:length(Hull_pts\$x)-1) {
Hull_lines[i] <- Lines(list(
Line(rbind(Hull_pts[i,1:2], Hull_pts[i+1,1:2])))
, ID=i)
Hull_lines[15] <- Lines(list(
Line(rbind(Hull_pts[15,1:2], Hull_pts[1,1:2])))
, ID="15")
}
# Make it a SpatialLines object
Hull_spLines <- SpatialLines(Hull_lines, proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

# Create SpatialPolygon of Hull (as alternative)
Hull_poly <- SpatialPolygons(list(
Polygons(list(
Polygon(Hull_pts))
, "ID"))
, proj4string=CRS("+proj=longlat +datum=WGS84"))

# Outer points
Site <- c(seq(1,10,1))
x <- c(38.21467,37.22799, 38.22852, 39.42621, 37.85457, 38.10600, 37.84077, 38.02475, 38.07228, 37.40716)
y <- c(6.327209, 5.691638, 6.237393, 5.382760, 6.042394, 6.078254, 6.007846, 5.973651, 6.044092, 5.746368)
Sites <- data.frame(Site, x, y)
xy <- Sites[2:3]
Sites2 <- SpatialPointsDataFrame(coords=xy, data=Sites)
projection(Sites2) <- CRS('+proj=longlat +datum=WGS84')
crs(Sites2) <- "+proj=longlat +datum=WGS84"

# Create lines
for (i in 1:length(Sites\$x)) {
Line(rbind(c(38.35419, 4.483533), Sites[i,2:3])))
, ID=Sites[i,1])
# Make SpatialLines object

# Plot
plot(Hull_spLines, xlim=c(37,40.5), ylim=c(2.5,6.6), axes=T)
points(Sites2\$x, Sites2\$y, pch=20, col="blue")
``````

What I need to find are the co-ordinates of the points where the radial lines cross the polygon, as achieved here: http://uk.mathworks.com/help/map/ref/polyxpoly.html. I know this is possible in ArcGIS, and it doesn't seem like it should be that difficult.

Does anyone have any ideas?

Other users might be interested to know that I have managed to solve this. Rather than extracting the points directly, I used `gIntersection` in the `rgeos` package to truncate the lines to within the polygon, and then extracted the outer point.

``````# Cut radial lines at polygon edge

# Extract edge points from lines
for (i in 1:length(Sites\$x)){     # i-th list element, 2 and 4 are edge co-ords (1 and 3 are centre point)
Sites\$Edge_x[i] <- coordinates(Centre_edge_lines)[[i]][[1]][2]
Sites\$Edge_y[i] <- coordinates(Centre_edge_lines)[[i]][[1]][4]
}
``````

Thanks for this answer. Just to clarify for other users, if you use coordinates(gIntersection(inputline, inputpoly)) it returns lists of coordinates. Each list is a line segment split by the polygon layer. If you then extract the first and last item from each of these list items, you will get the intersecting points. I wrote up a little generic function to to this. It could probably be simplfied, but it works.

``````intersectedline=gIntersection(SpatialLine, SpatialPoly, byid=T)

intersectpts=function(intersectedline, SpatialPoly){
x=2*rep(NA, length(intersectedline))
y=2*rep(NA, length(intersectedline))
xylast=cbind(x,y)
xy1=cbind(x,y)

for(i in 1:length(intersectedline)){
xy1[i,] <- coordinates(intersectedline)[[i]][[1]][1,]
lastpoint=length(coordinates(intersectedline)[[i]][[1]][,1])
xylast[i,] <- coordinates(intersectedline)[[i]][[1]][lastpoint,]
i=i+1
}

xy=rbind(xy1,xylast)

plot(SpatialPoly)
lines(intersectedline, col="red")
points(xy, col="blue")

return(xy)
}
``````