libLAS RFC 1: Vertical Coordinate Systems

Author

Frank Warmerdam

Contact

warmerdam@pobox.com

Date

1/6/10

Status

Implemented

Version

libLAS 1.3

Contents

• libLAS RFC 1: Vertical Coordinate Systems

  • Vertical CS in GeoTIFF

  • Vertical CS in OGC WKT

  • Vertical CS in PROJ.4

  • Sample Data

  • Proposed LASSpatialReference Interfaces

  • Proposed Changes to las2las

  • Changes in GDAL

  • Changes in libgeotiff

This page holds resources to handling of vertical coordinate systems (datums and units) in LAS files. Traditionally LAS files have not been rigorous with regard to recording the vertical datum and vertical units. In actual use I (Frank Warmerdam) have been unable to find any existing LAS files (January 2010) with vertical datum information set, though some do set vertical units information.

With the support of !LizardTech, I have been asked to extend libLAS to support reading and writing vertical coordinate system information in LAS files, and to the extent practical to establish some reference information on “best practices” with regard to describing vertical coordinate systems in LAS files.

Vertical CS in GeoTIFF

The planned approach is to store vertical coordinate system in GeoTIFF tags. The GeoTIFF specification has always had vertical coordinate system support, tied to the EPSG codes for vertical coordinate systems and datums. However, traditionally, the vertical coordinate system support in GeoTIFF has also not been widely used. So it is hoped that formalization of GeoTIFF tags for vertical datums in LAS format would also apply for GeoTIFF in general and any other geotiff oriented formats (GeoJP2 for instance). It is anticipated that the GeoTIFF representation of a LAS file with a vertical datum set might look something like this:

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsPoint
      ProjectedCSTypeGeoKey (Short,1): PCS_WGS84_UTM_zone_15N
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      VerticalCSTypeGeoKey (Short,1): Unknown-5701
      VerticalCitationGeoKey (Ascii,7): "Newlyn"
      VerticalDatumGeoKey (Short,1): Unknown-5101
      VerticalUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.

Note that VerticalCSTypeGeoKey 5701 is “ODN height”, a vertical coordinate system from the EPSG coordinate reference system table. The vertical coordinate system 5701 is based on vertical datum 5101 (Ordnance Datum Newlyn) from the EPSG datum table. The 5701 vertical coordinate system also specifies use of meters as the vertical linear measure.

So there are four GeoTIFF tags related to vertical coordinate system information:

• VerticalCSTypeGeoKey: an EPSG coordinate system code of type vertical (normally 5600-5799).

• VerticalCitationGeoKey: a description of the vertical coordinate system, often the EPSG user name for the VerticalCS.

• VerticalDatumGeoKey: An EPSG vertical datum code (normally 5000 to 5299).

• VerticalUnitsGeoKey: An EPSG units code (normally 9000 to 9199, 9001=meter)

Interesting things to note:

• There does not appear to be a VerticalUnitSizeGeoKey allowing user defined vertical units similarly to the ProjLinearSizeGeoKey, so it is not appropriate to use KvUserDefined for VerticalUnitsGeoKey.

• The libgeotiff epsg_vertcs.inc file seems to imply that the vertical datum codes in the range 5000 to 5299 should be considered vertical coordinate systems which they are not in the EPSG data model. It is not clear whether this is something that changed in the EPSG data model after the initiation of the GeoTIFF specification, or if it was an error by the GeoTIFF/libgeotiff authors. This mismatch is why Unknown-5101 above is not reported as Newlyn directly. Spec issue reported at http://trac.osgeo.org/geotiff/ticket/24

• The EPSG and OGC models include an EPSG code for a composite coordinate system. The composite coordinate system is a collection of a horizontal coordinate system (like PCS_WGS84_UTM_zone_15N) and a vertical coordinate system (like 5701 - ODN height). An example is [http://spatialreference.org/ref/epsg/7405/html EPSG:7405] - OSGB36 / British National Grid + ODN. There is no way to represent this in a GeoTIFF tag.

Vertical CS in OGC WKT

As noted above, EPSG and OGC have a concept of a compound coordinate system. An example from the CT 1.0 specification (OGC 09-001) looks like this:

COMPD_CS["OSGB36 / British National Grid + ODN",
    PROJCS["OSGB 1936 / British National Grid",
        GEOGCS["OSGB 1936",
            DATUM["OSGB_1936",
                SPHEROID["Airy 1830",6377563.396,299.3249646,
                    AUTHORITY["EPSG","7001"]],
                TOWGS84[375,-111,431,0,0,0,0],
                AUTHORITY["EPSG","6277"]],
            PRIMEM["Greenwich",0,
                AUTHORITY["EPSG","8901"]],
            UNIT["DMSH",0.0174532925199433,
                AUTHORITY["EPSG","9108"]],
            AXIS["Lat",NORTH],
            AXIS["Long",EAST],
            AUTHORITY["EPSG","4277"]],
        PROJECTION["Transverse_Mercator"],
        PARAMETER["latitude_of_origin",49],
        PARAMETER["central_meridian",-2],
        PARAMETER["scale_factor",0.999601272],
        PARAMETER["false_easting",400000],
        PARAMETER["false_northing",-100000],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["E",EAST],
        AXIS["N",NORTH],
        AUTHORITY["EPSG","27700"]],
    VERT_CS["Newlyn",
        VERT_DATUM["Ordnance Datum Newlyn",2005,
            AUTHORITY["EPSG","5101"]],
        UNIT["metre",1,
            AUTHORITY["EPSG","9001"]],
        AXIS["Up",UP],
        AUTHORITY["EPSG","5701"]],
    AUTHORITY["EPSG","7405"]]

The VERT_CS portion is analogous to the information we store in GeoTIFF file about vertical coordinate systems, though we have no mechanism to hold the overall name or authority code of the COMPD_CS wrapper.

Interesting things to note:

• The VERT_DATUM keyword includes a name and a datum type code (2005 in this case). The set from the specification is specified in on page 58, section 12.3.8 which is roughly captured in [wiki:OGCVerticalDatumTypes OGC Vertical Datum Types].

• For reference, the same compound coordinate system is handled similarly in [http://spatialreference.org/ref/epsg/7405/html Geotools]

• GDAL/OGR, used to process libLAS OGC WKT does not currently have meaningful support for vertical or compound coordinate systems.

Vertical CS in PROJ.4

PROJ.4 currently has no support for vertical coordinate systems, nor any way of representing them.

Sample Data

A standard sample file has been prepared demonstrating good practice setting a vertical cs. The geotiff tags look like:

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      ProjectedCSTypeGeoKey (Short,1): PCS_WGS84_UTM_zone_17N
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      VerticalCSTypeGeoKey (Short,1): Unknown-5703
      VerticalCitationGeoKey (Ascii,38): "North American Vertical Datum of 1988"
      VerticalDatumGeoKey (Short,1): Unknown-5103
      VerticalUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.

The test file is available at: http://hg.liblas.org/main/raw-file/tip/test/data/srs_vertcs.las

Proposed LASSpatialReference Interfaces

The LASSpatialReference class is used to manipulate the spatial references which are considered to include the vertical coordinate system. It is proposed to add the following interfaces.

LASSpatialReference::SetVerticalCS

This is a fairly low level method to set the vertical coordinate system information on the spatial reference (setting the appropriate geokeys), and updating the VLR.

/// Sets the vertical coordinate system using geotiff key values.
/// This operation should normally be done after setting the horizontal
/// portion of the coordinate system with something like SetWKT(),
/// SetProj4(), SetGTIF() or SetFromUserInput()
/// \param verticalCSType - An EPSG vertical coordinate system code,
/// normally in the range 5600 to 5799, or -1 if one is not available.
/// \param citation - a textual description of the vertical coordinate
/// system or an empty string if nothing is available.
/// \param verticalDatum - the EPSG vertical datum code, often in the
/// range 5100 to 5299 - implied by verticalCSType if that is provided, or
/// -1 if no value is available.
/// \param verticalUnits - the EPSG vertical units code, often 9001 for Metre.

void LASSpatialReference::SetVerticalCS( int verticalCSType,
                                        std::string const& citation,
                                        int verticalDatum,
                                        int verticalUnits )

It is mainly intended for those wishing to extend coordinate system information from another source (WKT, PROJ.4, an existing las file) with vertical cs data in a convenient fashion. This will be used to implement -a_vertcs in las2las.

LASSpatialReference::GetWKT()

This method will be modified to support read with or without the compound coordinate system.

enum WKTModeFlag
{
    eHorizontalOnly = 1,
    eCompoundOK = 2
};

/// Returns the OGC WKT describing Spatial Reference System.
/// If GDAL is linked, it uses GDAL's operations and methods to determine
/// the WKT.  If GDAL is not linked, no WKT is returned.
/// \param mode_flag May be eHorizontalOnly indicating the WKT will not
/// include vertical coordinate system info (the default), or
/// eCompoundOK indicating the the returned WKT may be a compound
/// coordinate system if there is vertical coordinate system info
/// available.
std::string GetWKT( WKTModeFlag mode_flag = eHorizontalOnly ) const;

Note, that while eHorizontalOnly is the default, this is just to avoid compatibility problems. eCompoundOK is the general version of the function and returns a COMPD_CS coordinate system with the vertical coordinate system if there is vertical information available, otherwise is just returns the normal horizontal coordinate system.

Internally this is implemented by improvements to GTIFGetOGISDefn() in GDAL.

LASSpatialReference::SetWKT()

This existing method will be changed to support setting vertical CS information if the passed in WKT is a compound coordinate system. The interface would not change, just the breadth of functionality. The actual implementation of this functionality is in GTIFSetFromOGISDefn() in GDAL.

Proposed Changes to las2las

It is proposed to add two new options to las2las

• -a_srs <coordinate system info>: This would override the coordinate system information as part of the las2las conversion, but without actually doing any geometric transformation - analygous to -a_srs in tools like gdal_translate or ogr2ogr.

• -a_vertcs verticalCSType citation verticalDatum verticalUnits: This option would allow overriding/setting the vertical coordinate system information during las2las processing. Mostly this would be used to set vertical coordinate system information when it was missing. The arguments are the epsg code numbers except for citation which is a string (empty to not set a citation).

The goal is to make las2las more useful as a “coordinate system fixing up” tool, including a fairly convenient ability to set the vertical coordinate system information on files that were missing it.

Changes in GDAL

The following changes have been made in GDAL to better support libLAS vertical coordinate system handling:

• SetFromUserInput() was modified to recognise COMPD_CS as a valid root so a compound coordinate system can be set this way. Mainly useful for las2las.

• GTIFGetOGISDefn() was extended to build a COMPD_CS coordinate system if there is vertical coordinate system info.

• GTIFSetFromOGISDefn() was extended to write vertical cs tags if there is vertical cs information in the WKT.

• OGRSpatialReference was extended with the StripVertical() method to remove COMPD_CS wrappers from a compound coordinate system (used by GetWKT() to strip off vertical stuff).

These changes will be in GDAL 1.7.

Changes in libgeotiff

The following changes have been made in libgeotiff to better support libLAS vertical coordinate system handling:

• GTIFValueName() now knows how to translate !VerticalUnitsGeoKey.

• GTIFNew() now reads existing geotiff file tags in a way that support update in place within the GTIF structure.

• libLASSortKeys() was rewritten to be less fragile in the face of malformed sets of geokeys, such as those common in LAS files with several “0” keys at the end.

These changes will be in libgeotiff 1.3.0. They are mostly useful for las2las.