[gdal-dev] Vertical and geocentric coordinate support in OGR/PROJ4

Even Rouault even.rouault at mines-paris.org
Wed Aug 24 14:38:04 EDT 2011 

Le mercredi 24 août 2011 02:48:44, Ben Discoe a écrit :
> Hi folks,
> 
> The issues of geocentric CS and vertical CS are distinct, but they are
> related in an important way: Geocentric coordinates are not very useful to
> the GIS world unless they can be converted to orthometric ("sea level")
> heights, which involves a vertical CS.
> 
> Currently, the GDAL/OGR/PROJ4 stack has limited support for geocentric &
> vertical CS, but the limitation of what works and what doesn't is
> difficult to guess from the documentation.  I'd like to contribute my
> understanding to date.
> 
> GDAL is actually doing quite a lot right, despite some major historical
> issues it has to deal with.  As far as i can tell:
> 
> 1. PROJ4 was never meant to handle vertical CS; as described on
> http://trac.osgeo.org/proj/wiki/VerticalDatums it's a late addition.
> 
> 2. The EPSG codes are a messy, incomplete patchwork of 2D and 3D coordinate
> systems, in which vertical CS also seems to be an awkward, late addition.
> 
> 3. OGC's WKT encoding (and hence OGR's) also wasn't quite designed to
> handle vertical CS; it can store a VERT_CS node, but that apparently isn't
> sufficient to actually define the transformation.
> 
> Understandably, as a result of the above history, the OGRSpatialReference
> implementation of VERT_DATUM nodes is somewhat awkward and ad-hoc.
> 
> The simplest way to create an OGRSpatialReference with a vertical CS is to
> use an import function, for example:
> 
> A. Working backwards from PROJ4:
>   srs.importFromProj4("+proj=longlat +datum=WGS84 +no_defs
> +geoidgrids=g2009conus.gtx");
> 
> or
> 
> B. Example of a compound, Paris-based CS from EPSG:
>   srs.importFromEPSG(7400);
> 
> The first example produces an SRS that works (assuming you have
> g2009conus.gtx installed on your PROJ_LIB folder).  The WKT has a section
> that looks like: VERT_CS["NAVD88 height",
>         VERT_DATUM["North American Vertical Datum 1988",2005,
>             EXTENSION["PROJ4_GRIDS","g2009conus.gtx"]],
>         AXIS["Up",UP]]
> 
> The PROJ4_GRIDS EXTENSION is what makes PROJ4 work.
> 
> The second example will not actually function.  It creates WKT like this:
> COMPD_CS["NTF (Paris) + NGF IGN69 height",
>     GEOGCS[...],
>     VERT_CS["NGF IGN69 height",
>         VERT_DATUM["Nivellement General de la France - IGN69",2005,
>             AUTHORITY["EPSG","5119"]],
>         AXIS["Up",UP],
>         UNIT["metre",1,
>             AUTHORITY["EPSG","9001"]],
>         AUTHORITY["EPSG","5720"]],
>     AUTHORITY["EPSG","7400"]]
> 
> It lacks a PROJ4_GRIDS, because in the GDAL data file for composite CS,
> compdcs.csv, the entry is: 7400,"NTF (Paris) + NGF IGN69
> height",4807,5720,1,0
> 
> This means the vertical CS is 5720, which appears in vertcs.csv as:
>    5720,NGF IGN69 height,5119,Nivellement General de la France -
> IGN69,9001,1,0,6499,,
> 
> Note there is no .gtx file present for 5720.  There are in fact only 3
> entries in vertcs.csv which contain a .gtx file, and hence only 3 that
> will produce valid results:
> 
>   3855,EGM2008 geoid height,1027,EGM2008
> geoid,9001,1,0,6499,9665,"egm08_25.gtx" 5703,NAVD88 height,5103,North
> American Vertical Datum
> 1988,9001,1,0,6499,9665,"g2003conus.gtx,g2003alaska.gtx,g2003h01.gtx,g2003
> p01.gtx" 5773,EGM96 geoid height,5171,EGM96
> geoid,9001,1,0,6499,9665,"egm96_15.gtx"

My understanding (based on http://trac.osgeo.org/geotiff/changeset/1893 ) is 
that the GTX files are hand-added in a vertcs.override.csv file and merged with 
the data that come from the original EPSG database. So it is "just" a matter 
of finding the appropriate grid, adding new entries in vertcs.override.csv and 
then running the script that regenerates the CSV files ...

> 
> Of those three (3855, 5703, 5773), there is not a single entry in
> compdcs.csv which uses them.  Hence, no compound CS set by EPSG code will
> actually work.

Actually, there's another way of building a compound CS from EPSG codes you 
probably didn't notice. You can use a string like 
"EPSG:horizontal_code+vertical_code" as a valid input to SetFromUserInput().

For example if you want EPSG 4326 with the EGM 96 geoid, you would do 
"EPSG:4326+5773".

Here's the output of testepsg utility :

$ testepsg "EPSG:4326+5773"
Validate Succeeds.
WKT[EPSG:4326+5773] =
COMPD_CS["WGS 84 + VERT_CS",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    VERT_CS["EGM96 geoid height",
        VERT_DATUM["EGM96 geoid",2005,
            AUTHORITY["EPSG","5171"],
            EXTENSION["PROJ4_GRIDS","egm96_15.gtx"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Up",UP],
        AUTHORITY["EPSG","5773"]]]

[...]

PROJ.4 rendering of [EPSG:4326+5773] = +proj=longlat +datum=WGS84 
+geoidgrids=egm96_15.gtx +no_defs 


> 
> The third way to create an OGRSpatialReference with a vertical CS is like
> this:
> 
>   srs4.SetVertCS("NAVD88 height", "North American Vertical Datum 1988");
>   srs4.SetExtension("VERT_DATUM", "PROJ4_GRIDS", "g2009conus.gtx");
> 
> AFAICT, the arguments to SetVertCS() are not a matter of functionality;
> they can be anything, but it seems good practice to put human-readable
> descriptions in there, ideally taken from a standard text like vertcs.csv.
>  (This example is a bit confusing, because the vertical datum is named
> 1988, but the latest version of it is 2009).  The main purpose of calling
> SetVertCS is to promote the CS to a compound CS with a VERT_DATUM node. 

Yes, that's syntaxic sugar. Technically SetCompoundCS() would be the most 
"straightforward" way of building a compound CS from an horizontal and a 
vertical CS. SetVertCS() has indeed that extra magic to save a step :

" * This method will ensure a VERT_CS node is created if needed.  If the
 * existing coordinate system is GEOGCS or PROJCS rooted, then it will be
 * turned into a COMPD_CS."

> The next call (SetExtension) can then decorate the VERT_DATUM with the
> 'extension' to make the SRS actually work, i.e. allow it to exportToProj4
> correctly so that OGRCreateCoordinateTransformation will work correctly.

Yes, the extension is needed if you want to make use of the vertical CS with 
proj4. The text descriptions attached to the COMPD_CS, VERT_CS and VERT_DATUM 
nodes are ignored in the proj4 world.

> 
> Here is a complete, functional example of converting a geocentric value to
> a geographic coordinate with orthometric height:
> 
>   OGRSpatialReference srs1, srs2;
>   srs1.importFromEPSG(4978);		// 4978 is the EPSG code for geocentric
> (ECEF) srs2.SetWellKnownGeogCS("WGS84");
>   srs2.SetVertCS("NAVD88 height", "North American Vertical Datum 1988");
>   srs2.SetExtension( "VERT_DATUM", "PROJ4_GRIDS", "g2009conus.gtx");
> 
>   double x, y, z;
>   x = -2691744.9992481042;	// A point on the ellipsoid
>   y = -4262401.8609118778;
>   z =  3894209.4515372305;
> 
>   // ecef -> geoid
>   OGRCoordinateTransformation *oct =
> OGRCreateCoordinateTransformation(&srs1, &srs2); oct->Transform(1, &x, &y,
> &z);
>   printf(" lon,lat,height %lf,%lf,%lf\n", x, y, z);
> 
> If it works, it prints the correct value:
> lon,lat,height -122.272778,37.871667,32.269611
> 
> If the VertCS/Extension is omitted, or if PROJ cannot find the .gtx file,
> or anything else goes wrong, you simply get the ellipsoidal height (close
> to 0.0).
> 
> In conclusion, although the vertical CS support in OGR/PROJ4 is a bit
> limited (and neither the limitations nor supported cases are documented),
> it can be made to work.  Hopefully, this email will prove useful as
> documentation for anyone else needing to do this transformation.

Perhaps you would be willing to initiate a Trac wiki page with your findings ?

> 
> As an aside, this whole discussion is purely about what works in the Open
> stack; i have absolutely no idea what sort of WKT (COMPD_CS / VERT_CS /
> VERT_DATUM) would work in the e.g. ESRI universe.

There seems indeed to have differences on how this is handled in the ESRI 
world. I see this comment in OGRSpatialReference::importFromWkt() 
implementation :

/* -------------------------------------------------------------------- */
/*      The following seems to try and detect and unconsumed            */
/*      VERTCS[] coordinate system definition (ESRI style) and to       */
/*      import and attach it to the existing root.  Likely we will      */
/*      need to extend this somewhat to bring it into an acceptable     */
/*      OGRSpatialReference organization at some point.                 */
/* -------------------------------------------------------------------- */

> 
> -Ben
> 
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev

More information about the gdal-dev mailing list