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,g2003p01.gtx"
5773,EGM96 geoid height,5171,EGM96 geoid,9001,1,0,6499,9665,"egm96_15.gtx"
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.
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. 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.
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.
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.
-Ben
_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
