[geotk] VerticalCRS transformation

[Martin Desruisseaux]

Hello Aaron

Le 09/03/11 05:50, Aaron Braeckel a écrit :
> Caused by: org.opengis.referencing.operation.OperationNotFoundException: No 
> transformation available from system "VerticalDatum[Geoidal]" to 
> "VerticalDatum[EGM2008 geoid]".
>     at 
> org.geotoolkit.referencing.operation.DefaultCoordinateOperationFactory.createOperationStep(DefaultCoordinateOperationFactory.java:697)
>     at 
> org.geotoolkit.referencing.operation.DefaultCoordinateOperationFactory.createOperation(DefaultCoordinateOperationFactory.java:273)
>     at 
> org.geotoolkit.referencing.operation.DefaultCoordinateOperationFactory.createOperationStep(DefaultCoordinateOperationFactory.java:1306)
>     ... 26 more
>
> I may be making the wrong assumption, but it seems as if these two 
> VerticalDatums should be transformable.  The only thing I've noticed that may 
> be a factor is the fact that EPSG 3855 references a 2008 geoid and the 
> DefaultVerticalDatum.GEOIDAL references a 2005 geoid.  Any thoughts on why 
> these two VerticalCRSs could not be transformed?

Actually the 2005 number in DefaultVerticalDatum.GEOIDAL is not a year, but is 
an enum from a code list defined by OGC 01-009 specification. Vertical datum 
types are:

    * 2000 - Other.
    * 2001 - Orthometric.
    * 2002 - Ellipsoidal.
    * 2003 - Barometric Altitude.
    * 2004 - Normal.
    * 2005 - Geoid Model Derived.
    * 2006 - Depth.

So the 2005 code list just said that this vertical datum is a geoid model, but 
doesn't said which one (I admit this is confusing, but those numerical codes are 
not defined by us...). The "EGM2008 geoid" datum however seems to describe more 
accurately a specific model.

Geotoolkit.org supports the transformation of *ellipsoidal* heights. The 
transformation of geoidal heights is more complicated: we basically need to 
transform the geoidal height to an ellipsoidal height, apply the projection (if 
any) then transform back the ellipsoidal height into a geoidal height. This is 
because all the projection mathematics are defined relative to ellipsoids, not 
geoids.

Transformations between geoidal and ellipsoidal heights require a geoid model. 
There is one defined in the geotk-referencing3D module (EarthGravitationalModel 
<http://www.geotoolkit.org/apidocs/org/geotoolkit/referencing/operation/transform/EarthGravitationalModel.html>), 
which supports only the /WGS 84 Earth Gravity Model/ at this time. This is a 
port of Fortran code that use spherical harmonics for doing its computation. The 
"EGM2008 geoid" is not yet implemented. An other current limitation is that 
while the mathematics are implemented as MathTransform implementation, they are 
not yet automatically recognized by CoordinateOperationFactory (see GEOTK-139 
<http://jira.geotoolkit.org/browse/GEOTK-139>).

As a workaround, maybe you could try some of the following?

    *

      For all vertical CRS, use ellipsoidal height instead than geoidal height.
      Geotk should work fine for all ellipsoidal heights. Inconvenient is that
      height values are less realist (an ellipsoidal height of -10 metres for
      instance may actually be 200 metres above the sea level). This is actually
      the approach I would recommand if you can easily change the vertical CRS.

    *

      Drop the vertical CRS from the GridCoverage2D which is read or from the
      target CRS to which the coverage is resampled. You can use the
      CRS.getHorizontalCRS(...) and CRS.getTemporalCRS(...) convenience methods
      for extracting the other CRS components.


     Martin