Here is how you could do this sort of thing in R. You'll almost certainly want to tweak the function to behave just how you'd like it. (My function, for instance, computes the closest point of intersection with a line and returns the great circle segment terminating at that point. If you want to use the second intersection, or the intersection with a polygon of area greater than x square kilometers, or just want the function to return the locations of all intersections, this will at least have gotten you started.)
Without further ado, here's the code:
## (0) Load required packages
library(maptools) ## For ContourLines2SLDF(), used for example data
library(sp) ## For spDists() and CRS()
library(rgeos) ## For gIntersection()
library(geosphere) ## For gcIntermediate(), yielding paths along great circles
## (1) Define a couple of functions ...
## ... first a utility function that converts a two-column matrix of
## x-y coordinates (like that returned by geosphere::gcIntermediate)
## to a SpatialLines object (as needed by, e.g., rgeos::gIntersection) ...
xyMatToSL <- function(mat, prj=CRS("+proj=longlat")) {
id <- paste0(sample(letters, 10), collapse = "")
linesList <- list(Lines(Line(mat), ID = id))
SpatialLines(linesList, proj4string = prj)
}
## ... and then a main function that computes a segment of a great
## circle beginning from a given focal point at a given initial
## bearing and travels until its first intersection with a line OR
## until it has gone maxDist meters.
focalToFirstGCIntersection <-
function(focalPt, lineLayer, initBearing=80, maxDist=1e6) {
## Draw a long line from pt at an initial bearing
destPt <- destPoint(focalPt, b = initBearing, d = maxDist)
gcSegment <- xyMatToSL(gcIntermediate(focalPt, destPt, n=100, addStartEnd=TRUE))
## Find nearest points of intersection or, if none, return point at given distance
candidatePts <- gIntersection(gcSegment, lineLayer)
destPt <- if(is.null(candidatePts)) {
destPt} else {candidatePts[which.min(spDists(focalPt, candidatePts))]}
## Construct a line from point to it
xyMatToSL(gcIntermediate(focalPt, destPt, addStartEnd=TRUE))
}
And here is a simple example of its application to a made up data set:
## (3) An example
## Example data: a point and some lines
prj <- CRS("+proj=longlat")
egPoint <- readWKT("POINT(.5 .5)", p4s=prj)
egLines <- ContourLines2SLDF(contourLines(volcano,nlevels=3), proj4string=prj)
## Compute great circle segments of maximum length 40 km.
thetas <- seq(10,360,by=20)
outLines <- lapply(thetas, function(X) {
focalToFirstGCIntersection(egPoint, egLines, initBearing=X, maxDist=40000)
})
outLines <- do.call(rbind, outLines)
## Plot results
plot(egLines)
plot(outLines, col = "blue", lwd = 1.5, add = TRUE)
plot(egPoint, pch = 16, col = "gold", add = TRUE)
## One way to output the lines, here as a shapefile in R's working directory
## df <- data.frame(row.names=names(outLines), theta=thetas)
## outLines <- SpatialLinesDataFrame(outLines, df)
## writeOGR(outLines, getwd(), "outLines", driver="ESRI Shapefile")