Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have recently been looking at airlines web pages that display their routes departing a certain city to all other cities they service. I would like to be able to create similar curved routes between points. Has anyone created scripts or functions that will generate the curved arcs like those displayed in this example?

Flight Paths

In PostGIS, is there an implementation of the ST_MakeLine that would allow you to specify the amount of curve to use when connecting 2 points?

While I am currently using PostGIS and QGIS, I would welcome hearing about other software options that might be able to create the same appearance.

share|improve this question
    
Anybody know of any nice implementations of this? Examples or whatever? –  Mark Boulder Jun 13 '12 at 12:24

4 Answers 4

up vote 17 down vote accepted

Creating great circles could give you the desired effect.

Maybe something like discussed on http://lists.osgeo.org/pipermail/postgis-users/2008-February/018620.html

Update:

I've followed up on this idea in "Visualizing Global Connections". It's a purely PostGIS-based solution using reprojection to create arcs.

SELECT ST_Transform(
  ST_Segmentize(
    ST_MakeLine(
      ST_Transform(a.the_geom, 953027),
      ST_Transform(b.the_geom, 953027)
    ), 
  100000), 
4326)

(The CRS definition for 953027 can be found here: http://spatialreference.org/ref/esri/53027/)

enter image description here

share|improve this answer
2  
I like the idea, although with great circles, the problem you run into is that at shorter distances you will still end up with a generally straight line. I would like to be able to control the amount of arc that I put in the line (ie- arclength = distance * 2). –  RyanDalton Jan 17 '11 at 9:00
    
Here is a good example of the problem with simply using great circles: gc.kls2.com/cgi-bin/… –  RyanDalton Jan 17 '11 at 9:01
1  
After some additional research, I found this post that may be useful in assisting this method. mail-archive.com/postgis-users@postgis.refractions.net/… –  RyanDalton Jan 17 '11 at 19:47
    
For future reader's use, I thought I would just go ahead and link to @underdark's recent blog post that covers this topic. underdark.wordpress.com/2011/08/20/… –  RyanDalton Sep 9 '11 at 4:51

The problem is to figure out how much to bend the arcs to enhance their visual resolution.

Here's one solution (among the many possible). Let's consider all the arcs emanating from a common origin. The arcs get most crowded here. To separate them the best, let's arrange it so they spread out in equally-spaced angles. It's a problem if we draw straight line segments from the origin to the destinations, because there will typically be clusters of destinations in various directions. Let's use our freedom to bend the arcs in order to space the departing angles as evenly as possible.

For simplicity let's use circular arcs on the map. A natural measure of the "bend" in an arc from point y to point x is the difference between its bearing at y and the bearing directly from y to x. Such an arc is a sector of a circle on which y and x both lie; elementary geometry shows that the bending angle equals one-half the included angle in the arc.

To describe an algorithm we need a little more notation. Let y be the point of origin (as projected on the map) and let x_1, x_2, ..., x_n be the destination points. Define a_i to be the bearing from y to x_i, i = 1, 2, ..., n.

As a preliminary step, assume the bearings (all between 0 and 360 degrees) are in ascending order: this requires us to compute the bearings and then sort them; both are straightforward tasks.

Ideally, we would like the bearings of the arcs to equal 360/n, 2*360/n, etc., relative to some starting bearing. The differences between the desired bearings and the actual bearings therefore equal i * 360 / n - a_i plus the starting bearing, a0. The largest difference is the maximum of these n differences and the smallest difference is their minimum. Let's set a0 to be halfway between the max and the min; this is a good candidate for the starting bearing because it minimizes the maximum amount of bending that will occur. Consequently, define

b_i = i * 360 / n - a0 - a_i:

this is the bending to use.

It's a matter of elementary geometry to draw a circular arc from y to x that subtends an angle of 2 b_i, so I'll skip the details and go straight to an example. Here are illustrations of the solutions for 64, 16, and 4 random points placed within a rectangular map

alt text

alt text

alt text

As you can see, the solutions seem to get nicer as the number of destination points increases. The solution for n = 4 shows clearly how the bearings are equally spaced, for in this case the spacing equals 360/4 = 90 degrees and obviously that spacing is exactly achieved.

This solution is not perfect: you can probably identify several arcs that could be manually tweaked to improve the graphic. But it won't do a terrible job and seems to be a really good start.

The algorithm also has the merit of being simple: the most complicated part consists of sorting the destinations according to their bearings.


Coding

I don't know PostGIS, but perhaps the code I used to draw the examples can serve as a guide for implementing this algorithm in PostGIS (or any other GIS).

Consider the following to be pseudocode (but Mathematica will execute it :-). (If this site supported TeX, as the math, stats, and TCS ones do, I could make this a lot more readable.) The notation includes:

  • Variable and function names are case-sensitive.
  • [Alpha] is a lower-case Greek character. ([Pi] has the value you think it ought to have.)
  • x[[i]] is element i of an array x (indexed starting at 1).
  • f[a,b] applies function f to arguments a and b. Functions in proper case, like 'Min' and 'Table', are system-defined; functions with an initial lower case letter, like 'angles' and 'offset', are user-defined. Comments explain any obscure system functions (like 'Arg').
  • Table[f[i], {i, 1, n}] creates the array {f[1], f[2], ..., f[n]}.
  • Circle[o, r, {a, b}] creates an arc of the circle centered at o of radius r from angle a to angle b (both in radians counterclockwise from due east).
  • Ordering[x] returns an array of indexes of the sorted elements of x. x[[ Ordering[x] ]] is the sorted version of x. When y has the same length as x, y[[ Ordering[x] ]] sorts y in parallel with x.

The executable part of the code is mercifully short--less than 20 lines--because over half of it is either declarative overhead or comments.

Draw a map

z is a list of destinations and y is the origin.

circleMap[z_List, y_] := 
Module[{\[Alpha] = angles[y,z], \[Beta], \[Delta], n},
    (* Sort the destinations by bearing *)
    \[Beta] = Ordering[\[Alpha]];
    x = z[[\[Beta] ]]; (* Destinations, sorted by bearing from y *)
    \[Alpha] = \[Alpha][[\[Beta]]]; (* Bearings, in sorted order *)
    \[Delta] = offset[\[Alpha]];
    n = Length[\[Alpha]];
    Graphics[{(* Draw the lines *)
        Gray, Table[circle[y, x[[i]],2 \[Pi] i / n + \[Delta] - \[Alpha][[i]]], 
             {i, 1, Length[\[Alpha]]}],
        (* Draw the destination points *)
        Red, PointSize[0.02], Table[Point[u], {u, x}]
    }]
]

Create a circular arc from point x to point y starting at angle \[Beta] relative to the x-->y bearing.

circle[x_, y_, \[Beta]_] /; -\[Pi] < \[Beta] < \[Pi] := 
Module[{v,  \[Rho], r, o, \[Theta], sign},
    If[\[Beta]==0, Return[Line[{x,y}]]];

    (* Obtain the vector from x to y in polar coordinates. *)
    v = y - x; (* Vector from x to y *)
    \[Rho] = Norm[v]; (* Length of v *)
    \[Theta] = Arg[Complex @@ v]; (* Bearing from x to y *)

    (* Compute the radius and center of the circle.*)
    r = \[Rho] / (2 Sin[\[Beta]]); (* Circle radius, up to sign *)
    If[r < 0, sign = \[Pi], sign = 0];
    o = (x+y)/2 + (r/\[Rho]) Cos[\[Beta]]{v[[2]], -v[[1]]}; (* Circle center *)

    (* Create a sector of the circle. *)
    Circle[o, Abs[r], {\[Pi]/2 - \[Beta] + \[Theta] + sign, \[Pi] /2 + \[Beta] + \[Theta] + sign}]
]

Compute the bearings from an origin to a list of points.

angles[origin_, x_] := Arg[Complex@@(#-origin)] & /@ x;

Compute the midrange of the residuals of a set of bearings.

x is a list of bearings in sorted order. Ideally, x[[i]] ~ 2[Pi]i/n.

offset[x_List] :=
Module[
    {n = Length[x], y},
    (* Compute the residuals. *)
    y = Table[x[[i]] - 2 \[Pi] i / n, {i, 1, n}];
    (* Return their midrange. *)
    (Max[y] + Min[y])/2
]
share|improve this answer
    
I should mention that this solution assumes the destinations more or less surround the origin. When this is not the case, the entire idea (of equally spaced bearings) is not a good one. But it can readily be fixed up by introducing some fake destinations within the angular gaps and later removing those destinations (and their arcs). This process can be automated by computing the average distance between bearings and using that to identify the large gaps, etc. –  whuber Jan 17 '11 at 19:09
    
Nice graphics. I wonder if airlines use an automated tool when they draw up the route maps shown in the back of their in-flight magazine. –  Kirk Kuykendall Jan 17 '11 at 20:14
1  
@Kirk They probably pay someone to do the cartography manually :-). I was inspired by this question to see whether a simple approach could create reasonably good graphics. The answer looks promising. These graphics, by the way, were produced by Mathematica 8 using its Circle and Point primitives and a little vector arithmetic to find the circle centers. –  whuber Jan 17 '11 at 20:33
    
I love the result you showed and me this is the way to go. I'll be honest though, I consider myself technical but I got a little lost in the formula you gave, and how to turn that into PostGIS code therefore becomes nearly impossible. Anyone out there have any ideas on how to translate whuber's concept into workable code? I will try to review and give it a go, but help would be greatly appreciated. –  RyanDalton Jan 27 '11 at 21:21
    
@whuber- Thanks for the updated pseudocode. We'll have to see if we can actually implement it in PostGIS. –  RyanDalton Feb 11 '11 at 0:16

Try ST_CurveToLine

Something like for example:

SELECT ST_CurveToLine('CIRCULARSTRING(1 1,5 3,10 1)'::geometry) as the_geom;

You can visualize this by coping the query into the text box and pressing Map1 at http://www.postgisonline.org/map.php

share|improve this answer
    
+1 Very nice function. –  underdark Jan 17 '11 at 15:58

I would think you'd just want to roll your own polyline with some vector math , http://en.wikipedia.org/wiki/B%C3%A9zier_curve, or if your gis has an ICurve interface.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.