# Curved point-to-point "route maps" [closed]

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?

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?

I am currently using PostGIS and QGIS.

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

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

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
]
``````

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 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
``````

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.