[Proj] Height and ellipse conversions

support.mn at elisanet.fi support.mn at elisanet.fi
Tue Jun 11 16:02:44 PDT 2013


Hello,

You might try -->

http://trac.osgeo.org/proj/wiki/VerticalDatums

Janne.

------------------------------------------------------------------------------------------

Ben Caradoc-Davies [Ben.Caradoc-Davies at csiro.au] kirjoitti: 
> 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
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj
> 




More information about the Proj mailing list