[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