[Proj] Height and ellipse conversions

[Ben Caradoc-Davies]

Are height adjustments expected when converting between ellipses?

I am trying to understand some unexpected behaviour I found in OGR (via 
the Python bindings) while validating the correctness of 3D conversions 
between ellipses. In the original form I was using OGC WKT definitions, 
but I have boiled the behaviour down to a simple Proj.4 example. I am 
using proj 4.7.0 (proj-bin 4.7.0-2 amd64 debian/sid).

The example:

WGS84 has an ellipsoid with semimajor axis of 6378137 m.
WGS66 has an ellipsoid with semimajor axis of 6378145 m.

So, let's take a point on the WGS84 ellipsoid at the intersection of the 
meridian and the equator, (6378137, 0, 0) in geocentric X, Y, Z metres. 
This point should be 8 m *below* the WGS66 ellipsoid, so I expect that 
if I convert it to WGS66 longlat, the height should be -8.0 metres 
exactly, by definition.

$ echo "6378137 0 0" | cs2cs +proj=geocent +ellps=WGS84 +no_defs +to 
+proj=longlat +ellps=WGS66 +no_defs
0dE	0dN 0.000

What? That can't be right. This height should be -8.000 in WGS66.

$ echo "6378137 0 0" | cs2cs +proj=geocent +ellps=WGS84 +no_defs +to 
+proj=longlat +ellps=WGS84 +no_defs
0dE	0dN 0.000

Yes, I get the same result if the destination is WGS84. This is correct.

$ echo "6378137 0 0" | cs2cs +proj=geocent +ellps=WGS66 +no_defs +to 
+proj=longlat +ellps=WGS84 +no_defs
0dE	0dN -8.000

But I get the *right* answer if I set the *source* ellps to WGS66. What?

$ echo "6378137 0 0" | cs2cs +proj=geocent +ellps=WGS66 +no_defs +to 
+proj=longlat +ellps=WGS66 +no_defs
0dE	0dN -8.000

And the *target* ellps does not seem to matter. This does not seem to 
make sense. In each case, it looks to me like the source ellipse is 
being used as the target. (And what does an ellipse mean for geocent 
anyway?)

Furthermore, no combination of ellps options results in a height 
adjustment when going from longlat to longlat with an ellipse 
adjustment. For example:

$ echo "0dE 0dN 0.0" | cs2cs +proj=longlat +ellps=WGS84 +no_defs +to 
+proj=longlat +ellps=WGS66 +no_defs
0dE	0dN 0.000

If I go straight geocent to geocent, I get different points, It looks 
like the height relative to the ellipse is preserved, but that is a 
different point:

$ echo "6378137 0 0" | cs2cs +proj=geocent +ellps=WGS84 +no_defs +to 
+proj=geocent +ellps=WGS66 +no_defs
6378145.00	0.00 0.00

Is this behaviour expected, or have I misunderstood?

Here is the Python version of the WGS84 geocent to WGS66 lonlat. I got 
the same behaviour with OGC WKT SpatialReference as well as Proj.4 
versions, which led me to the source and here.

import ogr
import osr
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(6378137.0, 0.0, 0.0)
sourceSR = osr.SpatialReference()
targetSR = osr.SpatialReference()
sourceSR.ImportFromProj4("+proj=geocent +ellps=WGS84 +no_defs")
print sourceSR
targetSR.ImportFromProj4("+proj=longlat +ellps=WGS66 +no_defs")
print targetSR
transform = osr.CoordinateTransformation(sourceSR, targetSR)
print point.GetX(), point.GetY(), point.GetZ()
point.Transform(transform)
print point.GetX(), point.GetY(), point.GetZ()

Sample output:

GEOCCS["Geocentric",
     DATUM["unknown",
         SPHEROID["WGS84",6378137,298.257223563]],
     PRIMEM["Greenwich",0]]
GEOGCS["WGS 66",
     DATUM["unknown",
         SPHEROID["WGS66",6378145,298.25]],
     PRIMEM["Greenwich",0],
     UNIT["degree",0.0174532925199433]]
6378137.0 0.0 0.0
0.0 0.0 0.0

Wrong height.

And here is an example of a WGS84 to WGS66 conversion:

import ogr
import osr
# http://spatialreference.org/ref/epsg/4979/prettywkt/
EPSG_4979_WKT = """GEOGCS["WGS 84",
     DATUM["World Geodetic System 1984",
         SPHEROID["WGS 84",6378137.0,298.257223563,
             AUTHORITY["EPSG","7030"]],
         AUTHORITY["EPSG","6326"]],
     PRIMEM["Greenwich",0.0,
         AUTHORITY["EPSG","8901"]],
     UNIT["degree",0.017453292519943295],
     AXIS["Geodetic latitude",NORTH],
     AXIS["Geodetic longitude",EAST],
     AXIS["Ellipsoidal height",UP],
     AUTHORITY["EPSG","4979"]]
"""
# http://spatialreference.org/ref/epsg/4891/prettywkt/
EPSG_4891_WKT = """GEOGCS["WGS 66",
     DATUM["World Geodetic System 1966",
         SPHEROID["NWL 9D",6378145.0,298.25,
             AUTHORITY["EPSG","7025"]],
         AUTHORITY["EPSG","6760"]],
     PRIMEM["Greenwich",0.0,
         AUTHORITY["EPSG","8901"]],
     UNIT["degree",0.017453292519943295],
     AXIS["Geodetic latitude",NORTH],
     AXIS["Geodetic longitude",EAST],
     AXIS["Ellipsoidal height",UP],
     AUTHORITY["EPSG","4891"]]
"""
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(0.0, 0.0, 0.0)
sourceSR = osr.SpatialReference()
targetSR = osr.SpatialReference()
sourceSR.ImportFromWkt(EPSG_4979_WKT)
print sourceSR
targetSR.ImportFromWkt(EPSG_4891_WKT)
print targetSR
transform = osr.CoordinateTransformation(sourceSR, targetSR)
print point.GetX(), point.GetY(), point.GetZ()
point.Transform(transform)
print point.GetX(), point.GetY(), point.GetZ()

The output:

GEOGCS["WGS 84",
     DATUM["World Geodetic System 1984",
         SPHEROID["WGS 84",6378137.0,298.257223563,
             AUTHORITY["EPSG","7030"]],
         AUTHORITY["EPSG","6326"]],
     PRIMEM["Greenwich",0.0,
         AUTHORITY["EPSG","8901"]],
     UNIT["degree",0.017453292519943295],
     AXIS["Geodetic latitude",NORTH],
     AXIS["Geodetic longitude",EAST],
     AXIS["Ellipsoidal height",UP],
     AUTHORITY["EPSG","4979"]]
GEOGCS["WGS 66",
     DATUM["World Geodetic System 1966",
         SPHEROID["NWL 9D",6378145.0,298.25,
             AUTHORITY["EPSG","7025"]],
         AUTHORITY["EPSG","6760"]],
     PRIMEM["Greenwich",0.0,
         AUTHORITY["EPSG","8901"]],
     UNIT["degree",0.017453292519943295],
     AXIS["Geodetic latitude",NORTH],
     AXIS["Geodetic longitude",EAST],
     AXIS["Ellipsoidal height",UP],
     AUTHORITY["EPSG","4891"]]
0.0 0.0 0.0
0.0 0.0 0.0

Note that the heights are the same, which is not expected.

Kind regards,

-- 
Ben Caradoc-Davies <Ben.Caradoc-Davies at csiro.au>
Software Engineer
CSIRO Earth Science and Resource Engineering
Australian Resources Research Centre