[gdal-dev] Geocentric Coordinate System Support

[Ben Discoe]

Rodolfo,

I don't understand.  I am building from a checkout of trunk from GDAL SVN,
which certainly appears to have the Geocentric support in it.  Is this not
the latest code?  Is there some separate 1.9 branch?  My question is how the
geocentric support and vertical datums are supposed to work, e.g. what are
valid arguments to SetVertCS which can produce a valid result when
transforming from geocentric?

-Ben

> -----Original Message-----
> From: Rodolfo Bonnin [mailto:rodolfobonnin at suremptec.com.ar]
> Sent: Tuesday, August 23, 2011 6:12 AM
> To: Ben Discoe
> Subject: Re: [gdal-dev] Geocentric Coordinate System Support
> 
> Hello Ben,
> 
> >From a previous question I asked, I inferred that geocentric will not
> be backported to 1.8 branch,  I think you should jump to version 1.9
> to get that feature. I hope it helps.
> 
> Regards,
> 
> 
> 
> 2011/8/19 Ben Discoe <ben at vterrain.org>:
> > Frank et al.,
> >
> > I am testing the Geocentric support; in particular, the ability to
> transform
> > from Geocentric to a SRS with a vertical datum.
> >
> > Since Geocentric was committed in March, and GDAL 1.8.1 was released in
> > July, i thought that 1.8.1 might contain the Geocentric support.  But,
it
> > does not seem to be there.  The ogr_fromepsg.cpp in 1.8.1 dates from
last
> > December.  So, for testing, i have built GDAL from today's SVN.
> >
> > Basic geocentric support seems to work fine:
> >
> >    OGRSpatialReference srs1, srs2;
> >    srs1.SetWellKnownGeogCS("WGS84");
> >    srs2.importFromEPSG(4978);  // geocentric (ECEF)
> >    OGRCoordinateTransformation *oct =
> > OGRCreateCoordinateTransformation(&srs1, &srs2);
> >
> >    double x, y, z;
> >    x = -122.272778;
> >    y = 37.871667;
> >    z = 0;
> >    oct->Transform(1, &x, &y, &z);
> >
> > This gives (-2691744, -4262401, 3894209) which is a reasonable-looking
> > geocentric value.  However, the first question is:
> >
> > 1. How is the input z (i.e. elevation) being interpreted?  The input is
> only
> > a 2D coordinate system, with no definition of what Z is.  Is it assuming
> > that z is an ellipsoidal height, i.e. in this case relative to the WGS84
> > ellipsoid?
> >
> > Ellipsoidal elevation is not very useful.  To be useful, we must be able
to
> > specify the vertical coordinates of the input, e.g. in order to
interpret
> is
> > as orthometric height, i.e. relative to a geoid.
> >
> > I know that the underlying support for Vertical Datums in PROJ.4 works,
> > using an SRS like:
> >        +proj=longlat +datum=WGS84 +no_defs +geoidgrids=g2009conus.gtx
> >
> > However, how to do this from OGR is not clear.  It might involve
> > OGRSpatialReference::SetVertCS, but i could not find any example
anywhere
> > (even googling) of how to call SetVertCS.  I tried to guess, by working
> > backwards:
> >
> >    srs.importFromProj4("+proj=longlat +datum=WGS84 +no_defs
> > +geoidgrids=g2009conus.gtx");
> >    srs.exportToPrettyWkt(...);
> >
> > Which produces:
> >  COMPD_CS["WGS 84 + Unnamed Vertical Datum",
> >    GEOGCS["WGS 84",
> >        DATUM["WGS_1984",
> >            SPHEROID["WGS 84",6378137,298.257223563,
> >                AUTHORITY["EPSG","7030"]],
> >            TOWGS84[0,0,0,0,0,0,0],
> >            AUTHORITY["EPSG","6326"]],
> >        PRIMEM["Greenwich",0,
> >            AUTHORITY["EPSG","8901"]],
> >        UNIT["degree",0.0174532925199433,
> >            AUTHORITY["EPSG","9108"]],
> >        AUTHORITY["EPSG","4326"]],
> >    VERT_CS["Unnamed",
> >        VERT_DATUM["Unnamed",2005,
> >            EXTENSION["PROJ4_GRIDS","g2009conus.gtx"]],
> >        AXIS["Up",UP]]]
> >
> > 2. How to specify the vertical with OGR (using SetVertCS, or any other
way)
> > such that it will result in PROJ.4 using the necessary GTX file?
> >
> > Thanks,
> > Ben
> >
> >> -----Original Message-----
> >> From: gdal-dev-bounces at lists.osgeo.org [mailto:gdal-dev-
> >> bounces at lists.osgeo.org] On Behalf Of Frank Warmerdam
> >> Sent: Wednesday, March 09, 2011 8:41 AM
> >> To: gdal-dev
> >> Subject: [gdal-dev] Geocentric Coordinate System Support
> >>
> >> Folks,
> >>
> >> Yesterday I introduced Geocentric coordinate system support in
> >> GDAL/OGR.  This includes OGRSpatialReference support for geocentric
> >> coordinate systems, the ability to looking them up from EPSG, and
> >> the ability to translate them to/from PROJ.4 format so they can be
> >> used for transformation.  This is an example of a geocentric coordinate
> >> system looked up from EPSG.  (Note the root is "GEOCCS")
> >>
> >> GEOCCS["WGS 84 (geocentric)",
> >>      DATUM["WGS_1984",
> >>          SPHEROID["WGS 84",6378137,298.257223563,
> >>              AUTHORITY["EPSG","7030"]],
> >>          AUTHORITY["EPSG","6326"]],
> >>      PRIMEM["Greenwich",0,
> >>          AUTHORITY["EPSG","8901"]],
> >>      UNIT["metre",1,
> >>          AUTHORITY["EPSG","9001"]],
> >>      AXIS["Geocentric X",OTHER],
> >>      AXIS["Geocentric Y",OTHER],
> >>      AXIS["Geocentric Z",NORTH],
> >>      AUTHORITY["EPSG","4328"]]
> >>
> >> I think there will be relatively little need for this coordinate
system,
> >> but I would appreciate feedback from anyone using it.
> >>
> >>    http://trac.osgeo.org/gdal/changeset/21916
> >>
> >> Best regards,