Coordinate transformations version #15

andyferris · 2016-05-13T05:29:54Z

Hi,

I've been doing a complete overhaul of Geodesy. Underneath it all, we have written a new, generic, transformations package called CoordinateTransformations (link).

Two stand-out features are

The ability to compose transformations together into a transformation chain, using ∘ or compose()
UTM/UPS support as a partial, native-Julia port of GeographicLib

These changes make for a pretty massive update to Geodesy. One possible use that we foresee of composed transformations is for doing SLAM (simultaneous localization and mapping) that fuses absolutely georeferenced data and local sensor data. But they can be convenient for simpler things, too. The UTM/UPS addition means less people would need to use Proj4.jl.

This PR also relies on CoordinateTransformations getting into METADATA, which in turn is waiting for Rotations.jl to do the same, but it is a perfect time to do the code review.

Bounds functionality seems to be generic computational geometry stuff, not related to geodesy per se.

* Move ellipsoid stuff into new file datums.jl * Removed XYZ type and all get* methods, since these will be replaced with something easier to use (hopefully!)

Progagated changes to datums through conversion. Removed default datum choice WGS84. Note: now use lower-case wgs84, etc

Also modified macros to work in this case.

* New transformation types like `ECEFfromLLA` that keep track of datum, origin shift for ENU, etc. * Now conversion constructors use `transform()` * `distance()` moved to a new file and made more general

* Introduced UTM and UTMZ data types. Both are 3D point representations, but UTMZ includes the zone and hemisphere information. * Introduced UTM(Z)fromLLA and LLAfromUTM(Z) transformations. * Introduced methods for constructing composed transformations between UTM/UTMZ and many data types. * Added constructor-based conversions for UTMZ. It was decided UTM-based constructors were a little messy, especially for ENU conversions where there are potentially two different zones to consider. Users can do `ENU(UTMZ(utm_origin, zone, hemisphere), point)` and `ENU(origin, UTMZ(utm_point, zone, hemisphere))`, etc for this case. Though, interestingly, we could *assume* the same zone and make a conversion safely (perhaps hemisphere will screw this up?) * Converted Charles Karney's MATLAB code for an accurate 6th-order expansion of the transverse Mercator projection. Unfortunately, this runs ~14 times slower than Proj.4, and will need rewriting.

Used Charles Karney's C++ code from GeograhicLib (MIT licensed). Added some missing distance with UTMZ coordinates.

Now UTM zone 0 is actually the UPS projection, similar to GeographicLib. Code ported from Charles Karney's C++ version (under MIT license).

* UTM now uses a port of C++ GeographicLib and is signficantly faster. * UPS is now implemented with the polar stereographic projection, bugifx * Independently generated test data for UPS and UTM * Improved accuracy for LLAfromECEF transformation (also from GeograhpicLib) * Added many docstrings * Various new methods as necessary * Minor bugfixes and asthetic improvements

* Many more methods for transformations, conversions and distances between UTM points * Some extra UTMZ methods, and UTMZfromUTM/UTMfromUTMZ * Now transverse_mercator_*() will error-out if relative longitude is more than 50 degrees.

* Extensive README for new features * Added CoordinateTransformations as a dependency

andyferris · 2016-05-13T06:03:41Z

@c42f @awbsmith @PaulBellette

c42f · 2016-05-13T07:38:03Z

This is pretty epic, looking forward to checking it out in detail!

andyferris · 2016-05-17T00:36:23Z

Also FYI @cffk - I hope it's OK with you that I ported some of your code to Julia (re-released under MIT license)

andyferris · 2016-05-17T00:43:02Z

OK, now we are one step closer - waiting for CoordinateTransformations to make it to METADATA.

cffk · 2016-05-17T01:04:30Z

@andyferris Yes, this is fine!

andyferris · 2016-05-17T06:11:58Z

Thanks!

PaulBellette · 2016-05-18T01:56:59Z

src/transformations.jl

+            error("Major radius is not positive")
+        end
+        if !(isfinite(f) && f < 1)
+            error("Minor radius is not positive")


Is it at all important that the minor radius is actually minor (e.g. say f = -3) or does that not matter for any of the functions?

Sorry, I'm not sure... we could ask Charles?

Negative f is allowed. So you could change major and minor to equatorial and polar. With the ecef conversions, f needs to be in (-inf,1). For the series approximation, the ellipsoid needs to be close to spherical, i.e., f in (-0.02,0.02), say.

coveralls · 2016-05-18T04:04:23Z

Coverage decreased (-30.04%) to 65.719% when pulling 1b4d116 on FugroRoames:CoordinateTransformations-version into f402145 on JuliaGeo:master.

andyferris · 2016-05-18T04:24:57Z

OK, now CoordinateTransformations is in METADATA, tests pass on Julia-0.4 (but not nightlies - problems upstream?). Any last comments, @c42f?

c42f · 2016-05-18T06:12:44Z

README.md

+```
+
+Here we have used the WGS-84 ellipsoid to calculate the transformation, but other
+datums such as `osgb36`, `nad27` and `grs80` are provided. All transformations


I'd say datums -> ellipsoids here.

c42f · 2016-05-18T06:39:45Z

src/points.jl

+    y::T
+    z::T
+    zone::UInt8
+    hemisphere::Bool # true = north, false = south


Could perhaps change this to isnorth as a shorter and more descriptive name?

c42f · 2016-05-18T07:33:29Z

Sorry for the silence, I've been quite sick all day, but I took a look through to get you unblocked.

General impression is that this seems really great, it's pretty much just what we need.

As a meta-comment, I still don't like the use of the word "datum" here. What you've actually got are transformations which are defined by ellipsoid parameters. A given ellipsoid may be used by many different datums. Eg, GRS80 is not a datum at all, but it's used as the ellipsoid for many datums, eg GDA94. An ECEF->LLA transformation with a given ellipsoid can be used to transform between ECEF and LLA coordinates in datums D1 and D2, provided the same ellipsoid is defined for both D1 and D2.

So I'd suggest sidestepping the issue of datum entirely in this iteration and just passing around ellipsoids instead. In principle datums are tags which should be attached to points, or used in a lookup table to construct a transformation between datums (see eg, pj_datums.c from proj4).

andyferris · 2016-05-19T01:40:03Z

@c42f Currently, what is passed around are type-tags:

export wgs84

immutable WGS84; end
const wgs84 = WGS84()
ellipsoid(::WGS84) = wgs84_el # itself a pre-compiled constant

So the "tag" was a thing that can get ellipsoids on demand, potentially even other information that helps you do the transform. I was really hoping it would encapsulate a datum, and be more useful in the future to do more general things.

Though, you do make a good point about grs80... it's not a datum but it is a tag that gives us access to enough information to do the implemented transformations. Not sure if that means we should rename all the fields of the types, but I'd be open to suggestions.

andyferris · 2016-05-19T02:30:36Z

Actually, @c42f, I did have a few more crazy ideas to run past you:

Can we/should we export ° = pi/180 and use radians in the LLA types? How would data normally be on disk? Degrees or radians?
We could also define immutable Hemisphere; isnorth::Bool; end and export const north = Hemisphere(true) and const south = Hemisphere(false), and overload show() accordingly.
Abusing juxtaposition and with (2) we can have 35 north == 35*north == UTMZone(35,true) or similar. (alternatively, abuse subtraction with zone = 25-south).

These changes could potentially make interaction/input/output a bit easier - what is displayed in a user-friendly way by show could become real syntax!

coveralls · 2016-05-19T04:58:11Z

Coverage decreased (-23.5%) to 72.305% when pulling 13e5407 on FugroRoames:CoordinateTransformations-version into f402145 on JuliaGeo:master.

coveralls · 2016-05-19T05:19:23Z

Coverage decreased (-19.3%) to 76.451% when pulling b0c3441 on FugroRoames:CoordinateTransformations-version into f402145 on JuliaGeo:master.

c42f · 2016-05-19T06:12:53Z

Sure, I saw what you did with the datum tags. I guess it can make sense to define the LLA->ECEF transformation for a given datum, though it's only the ellipsoid which is really required. The question then becomes: should it be the role of Geodesy to maintain a datum->ellipsoid mapping, given that there's any number of datums, but only relatively few ellipsoids. Given that Geodesy isn't enforcing that the datum of the point types matches with the transformations, this suggests we should supply a set of ellipsoids.

I went and looked up the datums provided. We have:

wgs84 - originally used the GRS80 ellipsoid, but now uses WGS84
osgb36 - Airy 1830 ellipsoid
nad27 - Clarke 1866
grs80 - not a datum, but used by NAD83, GDA94, etc etc

Uh oh, so now there's a conundrum - I doubt people care so much about the ellipsoid as the datum, so using ellipsoid names will probably just confuse people more...

I'll think about the other ideas, just wanted to mention the above as I've got a few other things to look at thisarvo.

c42f · 2016-05-19T06:19:53Z

BTW, in the current design I notice you don't actually use the datum tag in any downstream types, so having the tag as an actual type becomes mostly unnecessary.

cffk · 2016-05-19T11:48:12Z

"Can we/should we export ° = pi/180 and use radians in the LLA types?
How would data normally be on disk? Degrees or radians?"

My experience with GeographicLib is that it's preferable to keep
latitudes, longitudes, and azimuth stored as degrees. Then special
values (lat = 90° for the north pole, longitude boundaries for UTM,
azimuth = 90° for due east, etc.) are all represented exactly. In
addition range reduction for longitudes and azimuths is exact. A
counterexample is GeoTrans where the code for assiging UTM zones is full
of epsilons to patch up the inaccuracy of converting radians to degrees.

c42f · 2016-05-20T00:04:59Z

special values are all represented exactly

Ah, that's very interesting. I did wonder about this choice in GeographicLib. Radians are very attractive for mathematical reasons; degrees rather less but somewhat compensated for by being more familiar. It hadn't occurred to me that degrees have some attractive numerical properties as well. (Side note - in my personal experience of using GeographicLib, I've had data from an external source in radians anyway, so I had to convert to and from degrees to pass the data. Ah well.)

@andyferris - in any case, I'm not much of a fan of storing things differently from how they appear when printed, this is cute but it's just asking for confusion IMO. If we really wanted to support both, I think we'd have to make it a type parameter - mightn't be too bad but would introduce a little extra complexity.

Plus fixed uncovered bugs and some typos

coveralls · 2016-05-22T23:23:44Z

Coverage decreased (-13.4%) to 82.314% when pulling 828b45e on FugroRoames:CoordinateTransformations-version into f402145 on JuliaGeo:master.

coveralls · 2016-05-23T01:17:09Z

Coverage decreased (-13.4%) to 82.314% when pulling a5f367c on FugroRoames:CoordinateTransformations-version into f402145 on JuliaGeo:master.

andyferris · 2016-05-23T01:25:09Z

OK, going to merge this and tag a new Geodesy in METADATA.

andyferris · 2016-05-23T01:38:27Z

Thanks @cffk for the advice, too. I guess there is no one-size-fits-all approach in any case, and whether you care about points extremely near a UTM boundary appearing in the right zone is probably application-dependent, but I do appreciate that you have implemented the most correct way of doing things.

Andy Ferris added 12 commits April 11, 2016 12:05

Removed generic bounds functions in bounds.jl

f8480d6

Bounds functionality seems to be generic computational geometry stuff, not related to geodesy per se.

Began refactoring datums

d033f6b

* Move ellipsoid stuff into new file datums.jl * Removed XYZ type and all get* methods, since these will be replaced with something easier to use (hopefully!)

Refactored datums and introduced LatLon

76e0c7c

Progagated changes to datums through conversion. Removed default datum choice WGS84. Note: now use lower-case wgs84, etc

Switched testsuite to BaseTestNext

8a03061

Also modified macros to work in this case.

Added GRS80 datum

a25f9a2

Implemented AbstractTransformation interface

3f1ab43

* New transformation types like `ECEFfromLLA` that keep track of datum, origin shift for ENU, etc. * Now conversion constructors use `transform()` * `distance()` moved to a new file and made more general

Upgraded LLAfromECEF with more accurate transform

13aecc3

Used Charles Karney's C++ code from GeograhicLib (MIT licensed). Added some missing distance with UTMZ coordinates.

Added polar stereographic projection

7aecb68

Now UTM zone 0 is actually the UPS projection, similar to GeographicLib. Code ported from Charles Karney's C++ version (under MIT license).

More consistent UTM and UTMZ behaviour

2bee24d

* Many more methods for transformations, conversions and distances between UTM points * Some extra UTMZ methods, and UTMZfromUTM/UTMfromUTMZ * Now transverse_mercator_*() will error-out if relative longitude is more than 50 degrees.

README and REQUIRE updates

09d4963

* Extensive README for new features * Added CoordinateTransformations as a dependency

andyferris mentioned this pull request May 17, 2016

Register CoordinateTransformations JuliaLang/METADATA.jl#5205

Merged

PaulBellette reviewed May 18, 2016
View reviewed changes

Typo and travis kick

1b4d116

c42f reviewed May 18, 2016
View reviewed changes

andyferris force-pushed the CoordinateTransformations-version branch from 13e5407 to b0c3441 Compare May 19, 2016 05:16

andyferris force-pushed the CoordinateTransformations-version branch from b0c3441 to 409b24c Compare May 19, 2016 06:02

More complete tests coverage

828b45e

Plus fixed uncovered bugs and some typos

andyferris force-pushed the CoordinateTransformations-version branch from 409b24c to 828b45e Compare May 22, 2016 23:21

Changed hemisphere to isnorth, +docstrings,tests

c2910b0

andyferris force-pushed the CoordinateTransformations-version branch from a5f367c to c2910b0 Compare May 23, 2016 01:24

andyferris merged commit caa03df into JuliaGeo:master May 23, 2016

andyferris deleted the CoordinateTransformations-version branch May 23, 2016 01:39

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Coordinate transformations version #15

Coordinate transformations version #15

andyferris commented May 13, 2016

andyferris commented May 13, 2016

c42f commented May 13, 2016

andyferris commented May 17, 2016

andyferris commented May 17, 2016

cffk commented May 17, 2016

andyferris commented May 17, 2016

PaulBellette May 18, 2016

andyferris May 18, 2016

cffk May 18, 2016

coveralls commented May 18, 2016

andyferris commented May 18, 2016

c42f May 18, 2016

c42f May 18, 2016

c42f commented May 18, 2016

andyferris commented May 19, 2016

andyferris commented May 19, 2016

coveralls commented May 19, 2016

coveralls commented May 19, 2016

c42f commented May 19, 2016

c42f commented May 19, 2016

cffk commented May 19, 2016

c42f commented May 20, 2016

coveralls commented May 22, 2016

coveralls commented May 23, 2016

andyferris commented May 23, 2016

andyferris commented May 23, 2016

Coordinate transformations version #15

Coordinate transformations version #15

Conversation

andyferris commented May 13, 2016

andyferris commented May 13, 2016

c42f commented May 13, 2016

andyferris commented May 17, 2016

andyferris commented May 17, 2016

cffk commented May 17, 2016

andyferris commented May 17, 2016

PaulBellette May 18, 2016

Choose a reason for hiding this comment

andyferris May 18, 2016

Choose a reason for hiding this comment

cffk May 18, 2016

Choose a reason for hiding this comment

coveralls commented May 18, 2016

andyferris commented May 18, 2016

c42f May 18, 2016

Choose a reason for hiding this comment

c42f May 18, 2016

Choose a reason for hiding this comment

c42f commented May 18, 2016

andyferris commented May 19, 2016

andyferris commented May 19, 2016

coveralls commented May 19, 2016

coveralls commented May 19, 2016

c42f commented May 19, 2016

c42f commented May 19, 2016

cffk commented May 19, 2016

c42f commented May 20, 2016

coveralls commented May 22, 2016

coveralls commented May 23, 2016

andyferris commented May 23, 2016

andyferris commented May 23, 2016