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
]
