Vertical and geocentric coordinate support in OGR/PROJ4

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

Vertical and geocentric coordinate support in OGR/PROJ4

Ben Discoe
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

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Vertical and geocentric coordinate support in OGR/PROJ4

support.mn
Hello,

What I have noticed is that most GPS sensors disagree about
the geoid (MSL) height but the reference ellipsoid is always the
same between different GPS units (what ever make they might
be). It is also true that GNSS sensor can calculate reference
ellipsoids without any large effort. So I would presume that
the best vertical reference would be some ellipsoid but hardly
any geoid even if there are areas on the earth where the geoid
is well known and calculated.

But sure if you need horizontal datums you also need vertical
datums. The question is more what is the reference? I would
say that it is "towgs84" both horizontal and vertical in proj4.

But even airplanes do not care about that since they use flight
levels and they are always referenced to the current local
pressure and that changes all the time when the weather is
changing.  (All "FL" planes go lower or higher when the pressure
changes)

"A Flight Level (FL) is a standard nominal altitude of an aircraft,
in hundreds of feet. This altitude is calculated from the
International standard pressure datum of 1013.25 hPa
(29.92 inHg), the average sea-level pressure, and therefore
is not necessarily the same as the aircraft's true altitude either
above mean sea level or above ground level." (or any reference
ellipsoid)

http://en.wikipedia.org/wiki/Flight_level

Regards: Janne.

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

Ben Discoe [[hidden email]] kirjoitti:

> 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
>
> _______________________________________________
> Proj mailing list
> [hidden email]
> http://lists.maptools.org/mailman/listinfo/proj
>

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Vertical and geocentric coordinate support in OGR/PROJ4

Mikael Rittri
In reply to this post by Ben Discoe
Thanks, Ben, for the information.

We have some customers who need to convert between
height-above-geoid and height-above-WGS84-ellipsoid,
but they don't need much accuracy: a couple of
meters error is okay. My advice has been to use
some representation of EGM2008, as the most detailed
global geoid model.

I would like to learn more about the maximal (or typical)
deviation between EGM2008 and other geoid models. For
Europe, my impression is that most of the national height
systems differ less than 0.5 m from the common European
height system EVRF2007 (see "Map Projections for Europe",
page 101, http://www.ec-gis.org/document.cfm?id=425&db=document ).
A curious exception is a difference of 2.3 m for
Belgium, which uses (or has used) a height system
based on mean low tide of some kind, instead of
mean sea level. I think Ireland has some similar
height system.

But I am not sure about the difference between EVRF2007
and EGM2008. I have looked at a few sample points where
there was a difference of up to about 0.4 m. I haven't
found any conversion method between EVRF2007 and EGM2008,
or any estimate of the differences. Does anyone know?

For the rest of the world, I am even less sure of
things. I have heard that the true mean sea level
can differ by a couple of meters from the geoid,
since mean sea level is not an equipotential surface
(due to variations in water temperature, salinity,
and whatever). But since national heights can be
measured from low tide - or maybe from something else -
I guess one can expect even larger deviations.  

Regards,

Mikael Rittri
Carmenta
Sweden
http://www.carmenta.com 

-----Original Message-----
From: [hidden email] [mailto:[hidden email]] On Behalf Of Ben Discoe
Sent: den 24 augusti 2011 02:49
To: [hidden email]; [hidden email]
Subject: [Proj] Vertical and geocentric coordinate support in OGR/PROJ4

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

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: Vertical and geocentric coordinate support in OGR/PROJ4

Greg Troxel

Mikael Rittri <[hidden email]> writes:

> We have some customers who need to convert between
> height-above-geoid and height-above-WGS84-ellipsoid,
> but they don't need much accuracy: a couple of
> meters error is okay. My advice has been to use
> some representation of EGM2008, as the most detailed
> global geoid model.

That sounds reasonable.

> I would like to learn more about the maximal (or typical)
> deviation between EGM2008 and other geoid models. For
> Europe, my impression is that most of the national height
> systems differ less than 0.5 m from the common European
> height system EVRF2007 (see "Map Projections for Europe",
> page 101, http://www.ec-gis.org/document.cfm?id=425&db=document ).
> A curious exception is a difference of 2.3 m for
> Belgium, which uses (or has used) a height system
> based on mean low tide of some kind, instead of
> mean sea level. I think Ireland has some similar
> height system.
NGS has some podcasts/lectures available about height, and they are well
worth listening to:

  http://www.ngs.noaa.gov/corbin/class_description/NSRS_Modernization_0411.shtml

(and see related items).

One thing that wasn't clear to me at first is that 'geoid model' can
mean several things.  The most straightforward is about the distance of
a particular equipotential above the ellipsoid.  But one can also make a
model that gives the difference between orthometic heights in some datum
and ellipsoid heights.  Such a model implicitly includes distortions in
the vertical datum.  But if what you really want is to obtain orthometic
heights from GPS, that's what you want.

This page lists some of the US geoid models.  Note that if you compare
the values at a point, GEOID09 and EGM2008 will differ not only because
of only GEOID09 matching NAVD88, but because NAD83 ellipsoid heights are
different from WGS84 ellipsoid heights.   That might be the case in
Europe if there is a modern regional datum distinct from ITRF/WGS84.
Near me (42N 71W), GEOID09 is -27.3m and EGM2008 is -28.1m.  But, this
doesn't mean that the models differ by 0.8m -- because the NAVD88 zero
surface is not necessarily the WGS84 geoid.

http://www.ngs.noaa.gov/GEOID/

> But I am not sure about the difference between EVRF2007
> and EGM2008. I have looked at a few sample points where
> there was a difference of up to about 0.4 m. I haven't
> found any conversion method between EVRF2007 and EGM2008,
> or any estimate of the differences. Does anyone know?

From a quick scan of EVRF2007, it appears that it is conceptually
similar to NAVD88, in that it's based on leveling and has a single fixed
point.  (But it's apparently dynamic height rather than orthometric.)
So a geoid model that relates the EVRF2007 0 height surface to the WGS84
ellipsoid would not necessarily have the same values as EGM2008.

http://www.bkg.bund.de/nn_164706/geodIS/EVRS/EN/EVRF2007/evrf2007__node.html__nnn=true

> For the rest of the world, I am even less sure of
> things. I have heard that the true mean sea level
> can differ by a couple of meters from the geoid,
> since mean sea level is not an equipotential surface
> (due to variations in water temperature, salinity,
> and whatever). But since national heights can be

and persistent/average winds.

> measured from low tide - or maybe from something else -
> I guess one can expect even larger deviations.  

I think the key question is what people who say they want "height above
geoid" really want.  I would think people that don't have a significant
understanding of geodesy probably want "orthometric height in the same
datum as my other data", which is a reasonable thing to want.


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj

attachment0 (200 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Vertical and geocentric coordinate support in OGR/PROJ4

Noel Zinn (cc)
In reply to this post by Mikael Rittri
Thanks for the EVRF2007 link, Mikael.  I have already thanked Ben for his
post that charts a path between gravity-based verticals and geocentric
coordinates (ECEF), which provide geodetically rigorous, scalable
visualization in 3D without cartographic distortion.  The two (gravity and
ECEF) are inextricably linked.  It's where geospatial is bound to go and Ben
is forging a path.

Since my endeavor is a .com and not a .org, I am not adverse to using
commercial products that speed my work.  At the risk of being off-topic or
inappropriate, I will say that one of them is Blue Marble Geographic
Calculator, in which I have no financial interest, and which provides Proj.4
functionality plus much else at a price.  Blue Marble excels at vertical
transformations by supporting the following gravity models simply and
easily:

Australia - AUSGEOID 98, AusGeoid2009
Canadian Vertical Datum CVGD28
Colombia - GEOCOL 04
Denmark - DVR90
France and Corsica - RAC09, RAF09
Great Britain - OSTN02
Iberia - IGM 95, IGG2005
Japan - Japan Height Datum via GSIGEO2005
Local Geodetic Datum Ellipsoid Height
The Netherlands - NLGEO2004
South Africa - SAGEOID2010
United States - NAVD88, NGVD29 via Geoid 96, Geoid 99, Geoid 03, or Geoid
2009
Worldwide - EGM96, EGM08, OSU91A
Also supports local offset height models for vertical datum transformations

My point is not to tout a commercial product but to reaffirm the thrust of
Ben's post.  I'd say the .com world is ahead of the .org world in providing
ease of access to the vertical transformations that are required to bring
disparate parts of the world into a common, distortion-free reference
system.

Noel Zinn, Principal, Hydrometronics LLC
+1-832-539-1472 (office), +1-281-221-0051 (cell)
[hidden email] (email)
http://www.hydrometronics.com (website)

-----Original Message-----
From: Mikael Rittri
Sent: Wednesday, August 24, 2011 8:13 AM
To: PROJ.4 and general Projections Discussions
Subject: Re: [Proj] Vertical and geocentric coordinate support in OGR/PROJ4

Thanks, Ben, for the information.

We have some customers who need to convert between
height-above-geoid and height-above-WGS84-ellipsoid,
but they don't need much accuracy: a couple of
meters error is okay. My advice has been to use
some representation of EGM2008, as the most detailed
global geoid model.

I would like to learn more about the maximal (or typical)
deviation between EGM2008 and other geoid models. For
Europe, my impression is that most of the national height
systems differ less than 0.5 m from the common European
height system EVRF2007 (see "Map Projections for Europe",
page 101, http://www.ec-gis.org/document.cfm?id=425&db=document ).
A curious exception is a difference of 2.3 m for
Belgium, which uses (or has used) a height system
based on mean low tide of some kind, instead of
mean sea level. I think Ireland has some similar
height system.

But I am not sure about the difference between EVRF2007
and EGM2008. I have looked at a few sample points where
there was a difference of up to about 0.4 m. I haven't
found any conversion method between EVRF2007 and EGM2008,
or any estimate of the differences. Does anyone know?

For the rest of the world, I am even less sure of
things. I have heard that the true mean sea level
can differ by a couple of meters from the geoid,
since mean sea level is not an equipotential surface
(due to variations in water temperature, salinity,
and whatever). But since national heights can be
measured from low tide - or maybe from something else -
I guess one can expect even larger deviations.

Regards,

Mikael Rittri
Carmenta
Sweden
http://www.carmenta.com

-----Original Message-----
From: [hidden email]
[mailto:[hidden email]] On Behalf Of Ben Discoe
Sent: den 24 augusti 2011 02:49
To: [hidden email]; [hidden email]
Subject: [Proj] Vertical and geocentric coordinate support in OGR/PROJ4

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

_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj 


_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj