WIP: Great circle calculations

[anowacki]

This is a very draft implementation of great circle calculations for Geodesy. Suggestions for improvement very warmly welcomed.

There are a number of things which need to be addressed before this could be merged:

1. Currently, I simply include the source from GeographicLib.jl directly in the tree, since GeographicLib.jl is not registered and I don't know how to depend on it otherwise using the REQUIRE system. This means there are currently no tests for the GeographicLib part. (Eventually, either I will register GeographicLib and we can depend on it, or deprecate that repo in favour of this, refactor the code and move the tests.)
2. The forward calculation ((::GreatCircle)(point, azi, dist)) should probably return a point of the same type as input, but how should this be done when also returning backazimuth, distance and angular distance?
3. Should the calculations return a true struct rather than a named tuple?

[c42f]

Great to see! I think it would be fine to copy the source of your GeographicLib.jl in directly, and to develop it further here. It's not too much of a burden to have the code provided it's got good tests (the tests should definitely come with it). I invited you as a member of JuliaGeo :-)

As for API, how do we name GreatCircle? Is this actually a cache of coefficients related to geodesic calculations? If so it acts more like a set of great circles than a single one, so EllipsoidGeodesics or GreatCircles might make sense? (I wonder whether the differential geometers have a name for this kind of "geodesic bundle" thing. It's probably super obscure though...)

As for the set of function calls we support for use with GreatCircles, we should inspect https://geographiclib.sourceforge.io/html/classGeographicLib_1_1Geodesic.html and figure out how all those many operations will fit into a julia API. Are any of these so important that they should be blessed with the function call syntax? I think this is maybe not the case and we should just make them normal function names instead. This would help clarify what the return type should be as well:

• azimuth(circles, ll1, ll2) - compute the azimuth from ll1 to ll2 (inverse geodesic problem)
• geodesic_distance(circles, ll1, ll2) - compute ellipsoidal distance from ll1 to ll2 (inverse geodesic problem)
• geodesic_destination(circles, p1, azi, dist) (other name geodesic_endpoint?) - compute lat lon position along geodesic (the direct geodesic problem). Ideally this would also used for the direct problem in terms of arc length, but we'd need to use types to distinguish between dist in meters vs ang in radians. We could take a dependency on Unitful I suppose but I'm always a little hesitant there.
• great_circle(circles, ll1, ll2) (other names geodesic_line(circles,ll1,ll2) or geodesic(circles, ll1, ll2)?) - compute the geodesic line as in (GeographicLib::Geodesic::InverseLine)
• great_circle(circles, ll1, azi [,dist]) (GeographicLib::Geodesic::Line / DirectLine)

These are just suggestions. The general theme here is to name functions according to what they do rather than how they do it. This is kind of the opposite of the approach taken in GeographicLib, but I think it will lead to a more intuitive interface.

Naming is hard as always. I think for the set of geodesics and an individual geodesic Geodesics/GeodesicLine or GreatCircles/GreatCircle would make sense. Perhaps the latter as it's less jargon-y, though I think it might also be less correct (GeodesicLine models a line on the ellipsoid with a given starting point and parameterization; there's many different GeodesicLine with different starting points for the same great circle.)

[anowacki]

Thanks very much, Chris! I will take this on board and push something with these suggestions incorporated. Apologies that it will not be immediately.

---

Just to explain how I got started on the initial interface: I took inspiration from the syntax for transformations, where e.g. a UTMfromLLA is not a single transformation from one point to another, but a precomputed cache of possible transformations for a certain datum, zone, etc., which is then called on arguments. That sort of felt idiomatic Geodesy.

As for API, how do we name GreatCircle? Is this actually a cache of coefficients related to geodesic calculations? If so it acts more like a set of great circles than a single one, so EllipsoidGeodesics or GreatCircles might make sense? (I wonder whether the differential geometers have a name for this kind of "geodesic bundle" thing. It's probably super obscure though...)

Having a plural as a package name which provides a type (e.g., CuArrays->CuArray, Measurements->Measurement, Geodesics->Geodesic) is such a classic Julia thing that I think I would get a bit confused by a plural type (Geodesics), personally. I can absolutely see the logic, of course and perhaps others wouldn't get confused. I think naming of the 'cache' object is less crucial if we don't actually call the struct and instead use it as an argument as you suggest. GeodesicCache or something.

As for the set of function calls we support for use with GreatCircles, we should inspect https://geographiclib.sourceforge.io/html/classGeographicLib_1_1Geodesic.html and figure out how all those many operations will fit into a julia API. Are any of these so important that they should be blessed with the function call syntax? I think this is maybe not the case and we should just make them normal function names instead. This would help clarify what the return type should be as well:

• azimuth(circles, ll1, ll2) - compute the azimuth from ll1 to ll2 (inverse geodesic problem)
• geodesic_distance(circles, ll1, ll2) - compute ellipsoidal distance from ll1 to ll2 (inverse geodesic problem)
• geodesic_destination(circles, p1, azi, dist) (other name geodesic_endpoint?) - compute lat lon position along geodesic (the direct geodesic problem). Ideally this would also used for the direct problem in terms of arc length, but we'd need to use types to distinguish between dist in meters vs ang in radians. We could take a dependency on Unitful I suppose but I'm always a little hesitant there.

I can't argue with azimuth, geodesic_distance or geodesic_endpoint, and I actually prefer this form to the struct-calling version I wrote. I personally will usually want an angular/arc distance rather than one in m, and so that would need to be exposed. (arc_distance?)

Degrees are used for coordinates (longitude-latitude) in Geodesy, so it seems consistent to stick to that for angular quantities here, in my opinion, rather than radians. Is the current approach of distinguishing between input arc distance and geodesic distance via keyword argument not helpful here? I agree making the input unitful would not be ideal.

• great_circle(circles, ll1, ll2) (other names geodesic_line(circles,ll1,ll2) or geodesic(circles, ll1, ll2)?) - compute the geodesic line as in (GeographicLib::Geodesic::InverseLine)
• great_circle(circles, ll1, azi [,dist]) (GeographicLib::Geodesic::Line / DirectLine)

Here, you're constructing a cache for a particular 'geodesic line', so the question is then what do you want from it, I suppose. Assuming geodesicline = GeodesicLine(circles, ll1, ll2), you could then have:

• geodesic_endpoint(geodesicline; dist, angle) – calculate endpoint for a line and a distance or angle
• waypoints(geodesicline; npoints, dist, angle) – give either a fixed number of points along the line, or a set of points spaced by dist m or angle°.

Naming is hard as always. I think for the set of geodesics and an individual geodesic Geodesics/GeodesicLine or GreatCircles/GreatCircle would make sense. Perhaps the latter as it's less jargon-y, though I think it might also be less correct (GeodesicLine models a line on the ellipsoid with a given starting point and parameterization; there's many different GeodesicLine with different starting points for the same great circle.)

These paths aren't circles at all because we're on an ellipsoid, naturally, so we should probably avoid GreatCircle entirely and stick with Geodesic/Geodesics. This is my fault for introducing the term here.

I don't really know which is more or less correct, but to me constructing a GeodesicLine is fine for later use. I personally would prefer any object to be constructed using a type-like syntax (i.e., with uppercase starting letter) since that helps me realise that it returns a special type.

Other computed values

One of the reasons that I implemented GeographicLIb.jl the way I did was because each forward or inverse calculation gives you several values (forward and backazimuth, geodesic and arc distance, etc.) Having functions like azimuth is good, but perhaps we should also expose something that just returns everything as well? geodesic?

[anowacki]

It strikes me now that much of this naming problem would go away if we just stored the precomputed cache needed for the geodesic calculations in each Ellipsoid instance, rather than creating a separate type for what is the same thing in the end.

It would require adding extra fields to the struct and ellipsoid creation would be slower, but creating ellipsoids is rare in Geodesy anyway. Built-in datums/ellipsoids would then just be used for the azimuth, geodesic_distance functions:

p1, p2 = LLA(0, 0, 0), LLA(1, 1, 0)
azimuth(p1, p2, wgs84)
geodesic_distance(p1, p2) # datum assumed to be wgs84 like in distance?

A Geodesic could then be a cache object equivalent to a GeographicLib.GeodesicLine.

Thoughts?

(Edit: GeographicLib.Geodesic(1, 1/300) construction takes about 600 ns, 1.9 KB on my machine, and Ellipsoid(a="1", f_inv="300") takes 2 µs, 1.0 KB if that helps inform whether adding it to Ellipsoid is tenable.)

[c42f]

Thanks very much, Chris! I will take this on board and push something with these suggestions incorporated. Apologies that it will not be immediately.

Just to explain how I got started on the initial interface: I took inspiration from the syntax for transformations, where e.g. a UTMfromLLA is not a single transformation from one point to another, but a precomputed cache of possible transformations for a certain datum, zone, etc., which is then called on arguments. That sort of felt idiomatic Geodesy.

Yes, this is the style from CoordinateTransformations.jl. To a certain extent many unusual things can be considered coordinate transformations, but I think it gets confusing when generalized too much.

As for API, how do we name GreatCircle? Is this actually a cache of coefficients related to geodesic calculations? If so it acts more like a set of great circles than a single one, so EllipsoidGeodesics or GreatCircles might make sense? (I wonder whether the differential geometers have a name for this kind of "geodesic bundle" thing. It's probably super obscure though...)

Having a plural as a package name which provides a type (e.g., CuArrays->CuArray, Measurements->Measurement, Geodesics->Geodesic) is such a classic Julia thing that I think I would get a bit confused by a plural type (Geodesics), personally. I can absolutely see the logic, of course and perhaps others wouldn't get confused. I think naming of the 'cache' object is less crucial if we don't actually call the struct and instead use it as an argument as you suggest. GeodesicCache or something.

GeodesicCache would work for me, though I don't entirely love it.

Having the public API based on an Ellipsoid is definitely "the right thing" that the user wants most of the time, I just wasn't sure it can be made efficient. Hrm though having said that, I wonder whether we could convince the compiler to constant-fold the creation of the GeodesicCache:

geodesic_distance(e::Ellipsoid, ll1, ll2) = geodesic_distance(GeodesicCache(e), ll1, ll2)

GeodesicCache(e::Ellipsoid) = some_hopefully_pure_enough_stuff # Maybe with a sprinkling of Base.@pure 

I can't argue with azimuth, geodesic_distance or geodesic_endpoint, and I actually prefer this form to the struct-calling version I wrote. I personally will usually want an angular/arc distance rather than one in m, and so that would need to be exposed. (arc_distance?)

Sure, my list was just a representative example, not meant to be exhaustive.

Degrees are used for coordinates (longitude-latitude) in Geodesy, so it seems consistent to stick to that for angular quantities here, in my opinion, rather than radians. Is the current approach of distinguishing between input arc distance and geodesic distance via keyword argument not helpful here?

Yes of course; degrees for consistency. Keyword args are probably fine for this too, and may also be worth using those for azimuth just for clarity.

These paths aren't circles at all because we're on an ellipsoid, naturally, so we should probably avoid GreatCircle entirely and stick with Geodesic/Geodesics. This is my fault for introducing the term here.

Though technically inaccurate, I don't think it's too bad. Accessibility of terminology is a good thing. I could still go either way here.

One of the reasons that I implemented GeographicLIb.jl the way I did was because each forward or inverse calculation gives you several values (forward and backazimuth, geodesic and arc distance, etc.) Having functions like azimuth is good, but perhaps we should also expose something that just returns everything as well? geodesic?

For sure: we should have some function to return all the values which are byproducts of the calculation, rather than recomputing them multiple times. Perhaps that would be a lower-level function which is named according to the technicalities (geodesic_direct / geodesic_inverse etc). A named tuple return for that seems pretty ok to me.

[c42f]

It strikes me now that much of this naming problem would go away if we just stored the precomputed cache needed for the geodesic calculations in each Ellipsoid instance

To expand, I do really like the idea that the user only need specify the Ellipsoid (or Datum from which the ellipsoid can be extracted). Though I don't think this implies that we should stash coefficients in the Ellipsoid itself. Hopefully we can get the compiler to constant fold some things as mentioned above.

azimuth(p1, p2, wgs84)

I'd write as

azimuth(wgs84, p1, p2)

because the ellipsoid is a context for the calculation and I think it's fair to say that context objects are ususally put first in Julia. (This allows the position of the context to be fixed, where the other arguments may change in type or number.)

geodesic_distance(p1, p2) # datum assumed to be wgs84 like in distance?

Yes, that's rather tempting, perhaps we should allow a default.

A Geodesic could then be a cache object equivalent to a GeographicLib.GeodesicLine.

I'm not sure I quite understand this comment.

[anowacki]

Some initial testing suggests that it'll be hard to get the compiler to do the constant propagation, but I will persevere.

azimuth(p1, p2, wgs84) I'd write as azimuth(wgs84, p1, p2)

Me too, but currently one does distance(p1, p2, wgs84), so consistency probably demands that either the new functions ape the old, or distance gets a change in some future release.

For sure: we should have some function to return all the values which are byproducts of the calculation, rather than recomputing them multiple times. Perhaps that would be a lower-level function which is named according to the technicalities (geodesic_direct / geodesic_inverse etc). A named tuple return for that seems pretty ok to me.

That's a good idea—let's go with that.

A Geodesic could then be a cache object equivalent to a GeographicLib.GeodesicLine.

I'm not sure I quite understand this comment.

I just meant that the singular Geodesic would be opened up as a name for what GeographicLib calls GeodesicLines, in the case that the current (to-be-renamed) GreatCircle is not called directly as we're thinking.

[c42f]

Me too, but currently one does distance(p1, p2, wgs84), so consistency probably demands that either the new functions ape the old, or distance gets a change in some future release.

Oh right. It's so long since I tried actually using this package that I forgot about that 😬. TBH at this stage I'm still paying some attention to Geodesy out of interest rather than need.

I do think it would be completely fine to deprecate that version of distance and swap argument order. But I also can't remember why we had distance originally and looking at it now it seems like a misfeature which could be deprecated entirely.

[bridgewalker]

I am just implementing mm-accuracy great circle calculation myself from original Vincenty/Karney papers. Happy to contribute code if there is still need.

[anowacki]

@bridgewalker This is great to hear. You probably know that I ported Karney's GeographicLib as the basis for this PR (which I clearly haven't advanced for some time). You are welcome to use that, base your work on it, or if you have the time and interest, update this PR or make a new one based off your own work.

In the end, I'm just interested in getting geodesic distances into Geodesy, but other things keep getting in the way. If you were interested in getting this done, I would be delighted for any support you'd like to offer!

[bridgewalker]

Mine is done and tested. Also implemented some other useful formulas. Still needs integrating with Geodesy though.

@bridgewalker This is great to hear. You probably know that I ported Karney's GeographicLib as the basis for this PR (which I clearly haven't advanced for some time). You are welcome to use that, base your work on it, or if you have the time and interest, update this PR or make a new one based off your own work.

In the end, I'm just interested in getting geodesic distances into Geodesy, but other things keep getting in the way. If you were interested in getting this done, I would be delighted for any support you'd like to offer!

[anowacki]

Great! Are you able to share a link to a repo? Are you interested in taking forward this or another PR adding geodesic calculations?

[bridgewalker]

Yes happy to contribute some more. Just working on a fast "Point in polygon" check. But fairly new to githup. So I guess the best way to get my code up is to create a new branch? So I should have a proper file that links for you that links to Geodesy in a few days.

…

On Mon, 17 Aug 2020 at 15:56, Andy Nowacki ***@***.***> wrote: Great! Are you able to share a link to a repo? Are you interested in taking forward this or another PR adding geodesic calculations? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#53 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AQSHVRECPAUZXNCPM7JCFCDSBEZKNANCNFSM4I4OWU7A> .

-- Prof. Dr. Thilo Gross Professor of Biodiversity Theory Universität Oldenburg / HIFMB reallygross.de / biond.org