4

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.

  • It's been ages since I've done this, but a point, and angle (bearing) and a length is also know as a vector, so at point 10,10, 45 deg from North, draw a 100 KM long line. So that might help find a solution, also it's common that a line can be unbounded, so a line from 10,10 to 10,15 can be infinitely long, but it's defined by that small segment. – MichaelStoner Jun 4 '14 at 20:10
  • Two questions: (1) Do you actually want to get the line segments that connect the point to the nearest intersection with the contour (or do you really want long line segments that cross the contour); and (2) are you still open to an R-based solution? – Josh O'Brien Jun 6 '14 at 16:24
  • I prefer to have a longer line that allow for possible more interseption (e.g. if the first interseption is a small island infront of the coast I probably use the second interseption). yes I'm alwasy open to R solutions! – Filippo Jun 6 '14 at 20:21
3

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

enter image description here

0

You could do this using the Azimuth and Distance plugin.

From what I remember this worked well in QGIS 1.8.

  • Azimut and Distance is somehow close but I can just draw a line at the time... still better than nothing thoug. – Filippo Jun 4 '14 at 21:37
  • If I enter more than one angle in the 'Next vertex' box then 1) if 'Survey type' is set on 'boundary' I obtain a line with as many segments as angles, 2) if 'Survey type' is set on 'polar/radial' I obtain a single line corresponding to the last angle entered. Am I missing something here? – Filippo Jun 4 '14 at 21:44
  • Oh you're right, it's not working properly. Best to create your own line memory layer and use this, use polar/radial and press 'clear list' after drawing a line. This will keep the starting vertex intact. I've never had this problem before. – nhopton Jun 5 '14 at 9:42

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.