Compound CRS for Gridded Data Extension

[jamesresslernga]

Using a Compound CRS for gridded data.

How does GeoPackage specify a compound CRS in the gpkg_spatial_ref_sys table, which has a singular srs_id and no provision to specify the Vertical CRS?

While some Compound CRS have a unique EPSG code (https://www.opengis.net/def/datum/EPSG/0/9518 for the 4326+3855 compound SRS), most compound CRS do not have a unique code.

Reading the remarks in 2016-17cited below about lack of support for Compound CRS, seems to discourage using a Compound CRS in GeoPackage. Is that view still valid in 2024, or has a method been found in the gpkg tables to indicate a gridded dataset references a Compound CRS, other that the WKT definition itself?

Without using some novel encoding in gpkg_spatial_ref_sys, only stating the horizontal CRS would not inform GeoPackage clients that the CRS is 3D. The WKT2 definition of Compound CRS is defined in OGC 18-010r7 (WKT2). One option is for GeoPackage clients to detect a Compound CRS by examination of the WKT for “COMPOUNDCRS” and extract the VERTCRS . Is using WKT to find VERTCRS a consistent method to determine whether a GeoPackage gridded data set uses a 3D Compound CRS?

Previous Remarks.
In the gridded data extension OGC 17-066r1, Requirement 5 for gpkg-contents includes a remark :
NOTE: Ideally for elevation data the vertical datum for each pyramid of elevation will be specified. However, it is impractical to mandate this for a number of reasons, including the difficulty in testing whether a specific SRS has a valid vertical datum.

In the OGC Elevation Extension interoperability experiment report (OGC 16-094r3), comments considered using a HORIZ CRS and VERT CRS, but this was not implemented in the Gridded data extension. Section 10.3 on CRS in the 16-094r3 report comments on the lack of client support for compound CRS and appears to discourage use of Compound CRS in GeoPackage.

The intent of the 3D CRS requirement was to explicitly specify the horizontal datum, projection, and coordinate system as well as the vertical datum and height type within a single EPSG code. Several participants expressed concern with this approach. The main concern is the lack of 3D CRS implementations in the software used to both create and read OGC GeoPackages and the limited variety of 3D CRSs currently specified in the EPSG CRS registry.
…
While creating a custom WKT to describe a 3D SRS offers a great amount of flexibility to GeoPackage elevation extension producers, this approach may ultimately hinder interoperability. Software clients that currently read WKT strings will need to be expanded to interpret a wide variety of non-standard, custom WKT strings, and software clients that do not currently read the WKT will need to add support. There was a general consensus among IE participants that specifying a 2D CRS plus a vertical CRS for each elevation data tile table is the preferred approach.

Furthermore, in an issue Specifying 3D CRSs · Issue #19 · opengeospatial/geopackage-tiled-gridded-coverage (github.com), Jeff Yutzler concluded:

The concerns stem from pretty poor support for doing any kind of reprojection using 3D SRS (e.g. GeoTools can't do that, impacts on GeoServer and some Java-based desktop applications, amongst others) and the lack of defined EPSG (or other standards body) codes for a 3D equivalent for a lot of 2D codes.

It was suggested that you can just create your own 3D SRS. Sure, but implementations will always be patchy in how well they interpret more "creative" WKT; and we'd also need to look at srs_auth namespacing and likely other unintended consequences. For example, that would break SpatiaLite geopackage support - it doesn't read the WKT.

[desruisseaux]

Assuming that the policy is the same as spatial SQL databases such as PostGIS, the srs_id column of the spatial_ref_sys table is not necessarily EPSG codes. The values can be whatever the geopackage file producer wants. A mapping to EPSG codes can be specified, if desired, in the organization and organization_coordsys_id columns, but it can also be a mapping to IAU codes, or OGC codes, or any other organization. The srs_id values only need to be unique within the database (source: Simple feature access – Part 2: SQL option, §6.1.3).

For compound CRS not defined by EPSG, file producers should be free to use a srs_id value of their choice and declare the CompoundCRS using WKT. This is well covered by the WKT standard, and does not require "creative" WKT. I suggest to encourage users to leverage all standard WKT constructs that they need, without making the GeoPackage standard more complicated with workarounds for software having limited capabilities. The PROJ support of WKT 2 covers a large segment of the market. Java users can use Apache SIS, which has WKT 1and 2 support together with GML, and support the n-dimensional coordinate operations needed for CompoundCRS. They can also use PROJ-JNI, or even better use OGC GeoAPI for making themselves more library-independent.

[jratike80]

The GDAL command used in the comment of #675 does just that and inserts an unique value outside the EPSG code range into gpkg_spatial_ref_sys srs_id=100000. The whole record with all the attributes is

  srs_name (String) = WGS 84 (G1762) + EGM2008 height
  srs_id (Integer64) = 100000
  organization (String) = NONE
  organization_coordsys_id (Integer64) = 100000
  definition (String) = undefined
  description (String) = (null)
  definition_12_063 (String) = COMPOUNDCRS["WGS 84 (G1762) + EGM2008 height",GEOGCRS["WGS 84 (G1762)",DYNAMIC[FRAMEEPOCH[2005]],DATUM["World Geodetic System 1984 (G1762)",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Geodesy. Navigation and positioning using GPS satellite system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",9057]],VERTCRS["EGM2008 height",VDATUM["EGM2008 geoid"],CS[vertical,1],AXIS["gravity-related height (H)",up,LENGTHUNIT["metre",1]],USAGE[SCOPE["Geodesy."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",3855]]]
  epoch (Real) = 2024.123

[desruisseaux]

I do not know if the organization and organization_coordsys_id columns are required to be non null. But if null values are allowed in those two columns, it may be better to leave them empty instead of inserting NONE when there is not really any organization managing the code.

[jratike80]

Both the organization and organization_coordsys_id columns are "not null". Maybe instead of NONE the value could be custom or local but I think that for the client software it would not change the situation. What is important is that there is something that triggers the software to interpret the CRS from the definition_12_063 WKT string (and epoch).

EDIT: The organization NONE is already defined in the GeoPackage standard for the "undefined" coordinate systems, srs_id of -1 and 0

[desruisseaux]

What is important is that there is something that triggers the software to interpret the CRS from the definition_12_063 WKT string (and epoch).

There is nothing special to do. What the software should do is:

• Check if it recognizes the organization value:
  • If yes, use the organization_coordsys_id
  • If no, parse the WKT.

No need to change the Geopackage specification, except maybe write the above as a recommendation. However, I think that organization and organization_coordsys_id should be allowed to be null. Otherwise, a software could think that the dummy value is important and show a warning to the user saying that it does not know that authority and fallback on WKT. If the value is null, then the software knows that the authority does not exist and there is no need to bother the user.

[bosborn]

FWIW, our libraries and software currently attempt to parse the WKT first (definition_12_063 if provided). A lookup using organization/organization_coordsys_id is done as a backup.

Suggesting organization and organization_coordsys_id be modified to nullable would at least be a substantive change, possibly critical.

[desruisseaux]

Some softwares allow users to plugin custom authorities. For example, Apache SIS has this capability and GeoTools too if they didn't removed it. So the use of WKT or authority first is at implementation's choice, but if a software chooses to use organization_coordsys_id first and doesn't recognize the organization value, then it may mean that the user forgot to add a plugin or that some problem occurred during the plugin registration. So it is reasonable for those softwares to report a warning to the users about unrecognized organization value. In such cases, allowing null values would avoid spurious warnings.

[jratike80]

I fear that some programs rely just on srs_id and don't even check the organization_coordsys_id. The foreign data constrainst use srs_id and it is so easy to check the srs_id from the gpkg_contents and the gpkg_geometry_columns and believe that everything is clear if there is "4326" in all places and never consider that the CRS is not necessarily EPSG:4326.

EDIT Actually just this case is clear because it is mandatory that 4326 means EPSG:4326. So imagine any other number instead of 4326.

NONE is already used in the standard as the organization for the special srs_ids 0 and -1.

[desruisseaux]

Yes, taking srs_id as if they were directly EPSG codes is a common mistake observed in the field, at least with SQL databases. But it will be difficult to move to 3D CRS, meteorological data, geospatial in space, and more if we do not fix that misconception. It does not require any change in the standard, "only" explanations.

[jamesresslernga]

From GeoPackage SWG meeting on 13 Feb, the following is recommended for DGIWG Compound CRS for geographic coordinates of elevation coverage.

Set spatial_ref_sys to the following values:

Field Value
srs_name OGC CRS URL of the horizontal CRS, WGS 84 realization or ensemble (4326), http://www.opengis.net/def/crs/EPSG/0/####
srs_id 100001
organization DGIWG
organization_coordsys_id 100001
definition Null – not used
definition_12_063 WKT for compound CRS with horizontal and vertical CRS
epoch Coordinate epoch date (real number YYYY.DDDD, where DDDD is % of 1 year)

A different srs_id number for a compound CRS for a UTM zone may be used.

[jratike80]

When it comes to DGIWG profile, I think this part in the document https://portal.dgiwg.org/files/73489 is probably not what was meant:

If the optional AXIS terms are not present in the WKT, then the following default values are
assumed [OGC 1-009, section 7.3.2].
• Geographic Coordinate Systems: AXIS[“Lon”,EAST],AXIS[“Lat”,NORTH]

The definition for srs_name : http://www.opengis.net/def/crs/EPSG/0/4326 in the same DGIWG standard is as follows

GEOCCS["WGS 84",
  DATUM["WGS_1984",
  SPHEROID["WGS84",6378137,298.257223563]],
  PRIMEM["Greenwich",0],
  UNIT["degree",0.0174532925199433]]

Thus by reading the DGIWG standard literaly that would mean longitude-latitude.
Maybe it would be best to make it totally clear and set the axis mandatory in the DGIWG profile.

AXIS["Geodetic latitude (Lat)",north],
AXIS["Geodetic longitude (Long)",east],

EDIT
DGIWG is using OGC standard 1-009 "Coordinate Transformation Services" as a reference. This is an old standard (but still valid) and it really defines these defaults for the axis order:

If the optional AXIS terms are not present, then the default values are assumed. They are 
Geographic Coordinate Systems: AXIS[“Lon”,EAST],AXIS[“Lat”,NORTH]
Projected Coordinate System: AXIS[“X”,EAST],AXIS[“Y”,NORTH]
Geocentric Coordinate System: AXIS[“X”,OTHER],AXIS[“Y”,EAST],AXIS[“Z”,NORTH]

That does not match with the newer Geographic information — Well-known text representation of coordinate reference systems https://docs.ogc.org/is/18-010r7/18-010r7.html
If someone has contacts with DGIWG they should be asked to correct this misleading piece of information. I have not heard that OGC has plans to review the Coordinate Transformation Services standard.

[jamesresslernga]

The feedback and meeting on 2/13/2024 was very helpful. DGIWG web services technical panel (P5) is recommending adding 2 unique CRS codes for the DGIWG organization for 3D compound CRS for gridded data.

For WGS84 geographic SRS, use the following values.

Field: Value
srs_name: OGC CRS URL of the horizontal CRS, WGS 84 realization or ensemble (4326), http://www.opengis.net/def/crs/EPSG/0/#### 
srs_id: 100001 
organization: DGIWG 
organization_coordsys_id: 100001 
definition: Null – not used
definition_12_063: WKT for compound CRS with horizontal and vertical CRS 
epoch Coordinate epoch date (real number YYYY.DDDD, where DDDD is % of 1 year)

A different srs_id number for a compound CRS for a UTM zone may be used, such as the following.

Field: Value 
srs_name: The OGC CRS URL of the UTM zone http://www.opengis.net/def/crs/EPSG/0/####
srs_id: 100002 
organization: DGIWG 
organization_coordsys_id: 100002 
definition: Null – not used 
definition_12_063: WKT for compound CRS with horizontal and vertical CRS epoch Coordinate 
epoch date (real number YYYY.DDDD, where DDDD is % of 1 year)

This issue may be closed.

[jamesresslernga]

@jratike80 jratike80
Acknowledge the question about coordinate order guidance in DGIWG-126 1.0, section 7.2 that was raised in Feb. 2024. Opening new issue in DGIWG GPKG GitLab to address. That issue is not dependent upon resolution of the compound CRS for this story #676.

[rouault]

Acknowledge the question about coordinate order guidance in section DGIWG-126 1.0, section 7.2 that was raised in Feb. 2024. Opening new issue in DGIWG GPKG GitLab to address.

for clarity, CRS WKT standards in the last 10+ years require explicit AXIS definition in CRS definition, so the concept of default axis order for CRS WKT no longer makes sense. That said, one should keep in mind that the GeoPackage specification explictly overrides axis order for its geometry blob encoding : The axis order in WKB stored in a GeoPackage follows the de facto standard for axis order in WKB and is therefore always (x,y{,z}{,m}) where x is easting or longitude, y is northing or latitude, z is optional elevation, and m is optional measure. This ordering explicitly overrides the axis order as specified in the SRS metadata, applying Case 4 from OGC 08-038r7, Revision to Axis Order Policy and Recommendations