[GIS] How to automatically draw lines with a given angle from one single point in QGIS

Question

I need to draw a bundle of lines all originating from the same point and with a 10° angle between each line (e.g. lines bearing at 0,10,20 ..90..110..220..350).
I can't decide a priori the length of each line since I need to intercept the contour of another layer (i.e. coastlines) therefore I guess I will use a really long length.

Is there a way to automatically draw these lines in a shapefile? Do I have to write a specific script or are plugins available? I do not know python but I can work in R if needed.

I've seen similar question have already been asked but always with different softwares or in slightly different ways (e.g. lines from one point to another) so that none of them was useful to me.

Answer 1

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")