# [Proj] Height and ellipse conversions

Tue Jun 11 04:37:36 EST 2013

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

\$ 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)
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)
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,

--