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.

  • Anybody know of any nice implementations of this? Examples or whatever? – Mark Boulder Jun 13 '12 at 12:24
up vote 26 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


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_Transform(a.the_geom, 953027),
      ST_Transform(b.the_geom, 953027)

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

enter image description here

  • 4
    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
  • 1
    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
  • That's great!! Used in my project to draw lines between user checkins and venue locations, picked up from Forsquare – Lorenzo Barbagli Nov 23 '15 at 10:14

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.


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] :=
    {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
  • 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

  • I ended up trying this out to curve a set of "two point" linestrings. – Brent Edwards Apr 7 '15 at 18:40

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.

I ended up trying this out to curve a set of "two point" linestrings using the ST_CurveToLine function as suggested by @Nicklas Avén.

I passed the following 3 coordinate sets to the ST_OffsetCurve function :

  1. Start of the original line
  2. Midpoint of a line offset parallel to the original line
  3. End of the original line

I used the ST_OffsetCurve function to calculate the offset - 1/10th of the length of the original line in my example.

Here is the SQL I used to generate the curved lines from the original straight lines :

    ST_CurveToLine('CIRCULARSTRING(' || st_x(st_startpoint(the_geom)) || ' ' || st_y(st_startpoint(the_geom)) || ', ' || st_x(st_centroid(ST_OffsetCurve(the_geom, st_length(the_geom)/10, 'quad_segs=4 join=bevel'))) || ' ' || st_y(st_centroid(ST_OffsetCurve(the_geom, st_length(the_geom)/10, 'quad_segs=4 join=bevel'))) || ', ' || st_x(st_endpoint(the_geom)) || ' ' ||  st_y(st_endpoint(the_geom)) || ')') AS the_curved_geom
  • Really useful, but for some reason the result doesn't respect my srid. Any idea why? – DMS02 Oct 26 '16 at 14:13
  • Could you provide some more detail - srid of input geometry, output srid is missing, is different, errors generated (what application(s) - QGIS, PostgreSQL). – Brent Edwards Oct 27 '16 at 14:52
  • The table where I want to insert the resulting curved lines has a enforce_srid_geom constraint. When I execute the query I get an error saying this query violates that constraint. With a table without that constraint it works but then when adding it to QGIS it is listed with srid 0. My query: INSERT INTO test (the_curved_geom) SELECT [here your SQL] FROM lines – DMS02 Oct 28 '16 at 6:47
  • Try running the postgis.net/docs/ST_GeometryType.html and postgis.net/docs/ST_SRID.html functions on the geometry column (the_curved_geom) and check if there are conflicts with your test table and enforce_srid_geom. If so, you can transform the geometry/srid as needed or modify your test table/constraint. – Brent Edwards Oct 28 '16 at 12:38

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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