8

I've been musing (and doing prototype coding) for a 'projection free' geographic codebase with your basic point, line & polygon primitives.

Rather than dealing with all the sacrifices that come along with projecting to the plane, however, I'm writing algorithms that work directly on the surface of the ellipsoid.

One of the potential complications is that there are different kinds of "lines" out there:

  • (arcs of) great circles: shortest distance along the (constant-zero-elevation) surface between two points; should correspond exactly to line-of-sight paths.
  • rhumb lines: connect the two points with a path of constant direction -- for example, some state borders follow lines of latitude (which are not great circles).
  • curves: circular arcs (paths of constant distance from a given center-point); Bezier (not sure on correct re-interpretation in the context of a curved surface), etc.

Of the different kinds of paths (including ones I missed), which are important enough that they have an 'exact' representation, vs representing within an error bounds by short segments of a simpler path (e.g. short geodesic arc segments)?

Clarification edits: by 'exact' above, I mean parametric. In other words: computable to any desired accuracy, without a a densification-on-import step.

An edit, much later, to add a citation I've come across that closely parallels my own thoughts on the use of 3D unit vectors as a geographic primitive: A Non-singular Horizontal Position Representation (alt link). Best part? I didn't have to write it all out myself!

Improve this question
6
  • Truly exact representations of any of these objects on any ellipsoid (which is not a sphere) is next to impossible. The geodesics are no longer portions of great circles in general; the rhumb lines will be nasty no matter what; geodesic circular arcs will be particularly messy. Is there really any point to doing this that is worth a couple of orders of magnitude of additional computation for every operation?
    – whuber
    Commented Feb 4, 2011 at 15:46
  • Nothing is exact with that interpretation of the word -- how about "parametric" instead as a better word choice? (Also, as a side note: if I restrict myself to spheroids rather than general ellipsoids, the parametric representations are somewhat less messy.) However, it's true that a lot of those things remain extremely messy/difficult --- hence my question! I'm not interested in a system which destroys existing data quality, but I don't feel a need to represent curves nobody is using either.
    – Dan S.
    Commented Feb 4, 2011 at 17:53
  • I've worked with great circle line segments on a spheroid represented as a pole of rotation (point at lat/long) with a start angle, and end angle. I found it difficult to visualize the math used to manipulate them (quaternions). sciencedirect.com/…
    – Kirk Kuykendall
    Commented Feb 4, 2011 at 19:00
  • @Kirk: An easier-to-work-with representation (opinion here!) is to use orthographic 3D coordinates and represent the start/end points as vectors -- you're still using quaternions (to represent rotations in 3D) but they're much easier to think about.
    – Dan S.
    Commented Feb 4, 2011 at 19:35
  • @Dan - but with 3d coords, you'd need to densify in order to maintain constant elevation, wouldn't you?
    – Kirk Kuykendall
    Commented Feb 4, 2011 at 19:51

3 Answers 3

Reset to default
6

The question concerns what kinds of curves deserve an implicitly exact representation rather than a discretized approximation. The crux of the matter is this: to be successful, the class of curves you support in this manner must be closed under the class of curve- and polygon-creation operations supported in the GIS.

These operations include:

  • Buffering. In this process, you need to construct curves that are parallel to features. ("Parallel" means in the sense of maintaining a fixed distance.) This includes circles and portions thereof (for buffering points), oblique parallels (which are curves equidistant to geodesics on the spheroid, and can reduce to isolated points in special cases), and concentric circles. On the sphere (but not, generally, on an ellipsoid) the oblique parallels are themselves circles.

  • Polygons of influence (Thiessen polygons; Voronoi polygons; Dirichlet cells). To construct the Thiessen polygons for a collection of point features we need to find bisecting lines, which are geodesics (they are straight); but for a collection of other kinds of features, such as points and segments, the boundaries of the Thiessen polygons include portions of parabolas (in the plane). Maybe you don't want to support this...

  • Set-theoretic overlays (intersection, union, difference, complement). These operations do not create any new kinds of curves.

  • Parallel translation and rotation. These are usually not possible to perform exactly on an ellipsoid (because it is not a homogeneous space), but are straightforward on the sphere. On the sphere, these operations do not create new kinds of curves.

The really problematic class of curves you propose consists of the general rhumb lines (loxodromes). Lines of latitude are rhumb lines but (on the sphere at least) they are also circles, so they present no additional problem. But general rhumb lines are complicated beasts: if they are not meridians or parallels, they spiral into one pole or the other. Buffers and parallel translations of rhumb lines will be genuinely new types of curves. You would have to represent these results as broken segments of lines and circles, which would defeat your purpose (and be fairly difficult to compute). Therefore I suggest not trying to support rhumb lines exactly.

In sum, it looks like you can be successful in your program if (a) you work on a spherical model of the earth rather than the more general ellipsoidal ("spheroidal") model and (b) you limit certain constructions such as Thiessen polygons (and medial axes, which are closely related) to collections of points.

Improve this answer
23
  • I really like this way of thinking about it. A few random thoughts in following comments...
    – Dan S.
    Commented Feb 11, 2011 at 18:15
  • Nit: Buffering is fixed-distance, not parallel (except for infinite lines).
    – Dan S.
    Commented Feb 11, 2011 at 18:17
  • On the plane, areas of influence for straight-sided features have straight-line and conic-section curves in their borders -- but I doubt that it is closed. I've actually no idea what curves form the boundary of influence for (conic-section) curved features themselves; maybe some deep thought/research will reveal that they're conic sections too, but I'm skeptical. Being closed in general under these operation seems very unlikely.
    – Dan S.
    Commented Feb 11, 2011 at 18:19
  • Set-theoretic operations: Planar GIS systems typically aren't closed under set complement since it creates an infinitely-sized polygon; sphere/sphereoid/ellipsoid can do better. Note that sometimes intersection points can't (or a difficult) to represent for a lot of curve typess, though, even on the plane.
    – Dan S.
    Commented Feb 11, 2011 at 18:34
  • Translation/rotation: Affine transforms are possible on the plane, although there are other non-affine transform possible that might actually make more sense by making them specifically non-affine. E.g: "move every point due north by 150 meters" is often what is meant by a simple translation in a lot of projections, but of course the distortions of the projection mean that intent is slightly undermined...
    – Dan S.
    Commented Feb 11, 2011 at 18:37
4

Dan,

You may be interested in some of the work I've been doing on geodesics. This is described in this preprint. In particular, note:

  • The direct and inverse geodesics problems may be solved to machine precision. This means about 15 nm for double precision. I can switch to long doubles, add an extra term in the series, and get accuracy of 6 pm. Note in particular, that the solution to the inverse problem always converges (unlike Vincenty's method). The speed is comparable to Vincenty's method (direct solution is somewhat faster, the inverse solution is somewhat slower).

  • I compute the reduced length and geodesic scales. These quantities give the differential properties of geodesics, and allow various geodesic problems (intersections, median lines, etc.) to be solved quickly and accurately using Newton's method. The curvature of the boundaries of buffer regions can be expressed in terms of these quantities. See this note that I sent to the proj.4 mailing list.

  • I define an ellipsoidal gnomonic projection in which geodesics are very nearly straight. This allows problems on the surface of the ellipsoid to be mapped to problems in plane geometry. For example, the intersection of two geodesics can be found exactly by estimating an intersection point, performing a gnomonic projection about that point, re-estimating the intersection and iterating.

  • I give expressions for the area of a geodesic polygon. There's no need to subdivide long edges to get an accurate result.

  • GeographicLib (on sourceforge) implements the algorithms.

  • Finally, I remark that for many purposes geodesics are preferable to any other sort of curve (in particular great ellipses or rhumb lines) because they obey the triangle inequality. This has several consequences:

    • Geodesics and geodesic circles intersect at right angles.
    • The shortest XXX line between a point P and an arbitrary curve C intersects C at right angles only if XXX = geodesic.
    • Geodesics are a natural way to partition data using a quadtree because bounds can be place on the range of distances between a arbitrary point and all the points within a quadnode.
Improve this answer
6
  • Thank you for sharing these ideas; it's a very nice paper containing some great work. Welcome to our site! Some remarks and questions follow. (1) In what sense do "geodesics and geodesic circles intersect at right angles"? You must have some restriction in mind because this is not generally true. (2) The 3D Euclidean distances @Dan refers to will also satisfy the triangle inequality.
    – whuber
    Commented Jul 11, 2011 at 22:09
  • (1) Consider all the geodesics emanating from a point A; this defines one family of curves. Consider next all the geodesic circles centered at A; this defines a second family of curves. Because of the triangle inequality, these two families are orthogonal. This gives one of the basic properties of the azimuthal equidistant projection. (2) Yes you're right, of course. If you're thoroughly into the world of the geometry of surfaces, you recoil at the idea of any properties which depend on how the surface is embedded in 3d space. (See Gauss' "Remarkable Theorem".)
    – cffk
    Commented Jul 12, 2011 at 3:57
  • Actually (1) is more general: in any Riemannian manifold, a geodesic from a point P outside a smooth curve c to a point on c that minimizes the distance between P and c must be orthogonal to c. Your statement about "geodesic circles" immediately follows (provided P is the center: that's the previously unstated restriction). I agree with the sentiment expressed after (2), but we must remember that the objective here is to perform accurate efficient calculations, rather than investigate intrinsic properties of the surface. A well-chosen embedding might facilitate that.
    – whuber
    Commented Jul 12, 2011 at 12:56
  • A belated thanks for this reply. :) I'm too underwater to do more than skim at the moment, but it looks to be a a quite fantastic trove. A quick note on geodesic intersections since you called them out -- mostly for you to review as a check on my own undercaffinated intuition: Exact intersections of spherical geodesics can be found easily by intersecting the planes of the corresponding great circles, and that result carries over to ellipsoids by using an auxiliary sphere -- or am I missing something there?
    – Dan S.
    Commented Jul 20, 2011 at 0:28
  • @Dan I composed an answer to your question. But it seems that comments are limits in length. So see the next answer instead.
    – cffk
    Commented Jul 20, 2011 at 20:50
2

This is the answer to @Dan's question about using the auxiliary sphere to solve intersection problems.

No, the auxiliary sphere doesn't let you solve for intersections directly. The problem is that the mapping from the ellipsoid to the sphere depends on the geodesic (e.g., its azimuth at the equator). Thus the auxiliary sphere is good for solving for a single geodesic but not for solving any problem involving more that one geodesic.

My recommended way to do geodesic intersections and interceptions (shortest path between a point and a geodesic) is to use the gnomonic projection. Geodesics on a sphere map to straight lines in the gnomonic projection and so, provided that your problem is confined to a hemisphere, the gnomonic projection converts these problems to 2d ones.

There is no projection which preserves the straightness of geodesics for an ellipsoid (because its curvature is not constant). However, GeographicLib provides a generalization of the gnomonic projection in which geodesics are very nearly straight. This leads to rapidly converging algorithms for geodesic intersection and interceptions (again provided that the points are all well within a hemisphere). See my answers (with code!) to gpesquero at

https://sourceforge.net/projects/geographiclib/forums/forum/1026621/topic/4085561

Finally, I'd like to point out that I've recently converted GeographicLib's geodesic routines to Javascript, so that you can play around with them in Google Maps. See

http://geographiclib.sourceforge.net/scripts/geod-google.html

http://geographiclib.sourceforge.net/scripts/geod-calc.html

(I didn't convert the gnomonic projection to Javascript yet. That would be reasonably straight forward. I would also convert the azimuthal equidistant projection at the same time, since that's a convenient way to solve another class of geodesic problems involving "median lines".)

ADDENDUM (2014-08-19)

It's also possible to solve for the time of closest approach for two vessels traveling at constant speed along geodesics. Because we know the differential properties of geodesics, it's possible to use Newton's method to get an accurate solution in a few iterations. Code to implement this is posted at

https://sourceforge.net/p/geographiclib/discussion/1026620/thread/33ce09e0

Improve this answer
3
  • I'm quite sure I'm off track, but not sure where; I'd love it if you could help debug my reasoning. (Next comment.) Otherwise: thanks a heap for the helpful code + commentary + link; it's tremendously useful.
    – Dan S.
    Commented Jul 21, 2011 at 0:23
  • Here goes, as concisely as I could make it. My reasoning is expressed in 3D Cartesian, not angular coords: (a) On a sphere, all points in a great circle are coplanar. (b) The transformation to the auxiliary sphere is linear and invertible. (Mistaken thinking?) (c) All points in an elliptical geodesic transform to points along a great circle on the aux. sphere. (d) All points on an elliptical geodesic are coplanar as well, due to (b). Finally, (e): Due to coplanarity, two candidate geodesic intersection points on the ellipsoid can be found by plane intersection.
    – Dan S.
    Commented Jul 21, 2011 at 0:29
  • 1
    @Dan, An ellipsoidal geodesic does not lie in a plane. (If it did, it would have to be a plane curve; and yet we know that, in general, on each circuit of the globe the geodesic falls short by an amount O(f).) The fault in your reasoning is (b) -- the connection between the aux. sphere and the ellipsoid is not linear. The latitude transformation is equivalent to a stretch in the z direction and so is linear. However the longitudes are related by a elliptic integral and this prevents there being any simple linear relationship.
    – cffk
    Commented Jul 21, 2011 at 3:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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