Coordinate reference systems

heidivanparys commented 3 years ago

At the Danish Agency for Data Supply and Efficiency, we are looking at GeoPackage for distributing (large) datasets. We have been experimenting, and a couple of questions have come up regarding coordinate reference systems.

The implementing rules on the interoperability of spatial data sets and services specify what coordinate reference systems must be used. In the data specifications, that have GML as default encoding, this was proposed implemented as follows (see also http://inspire.ec.europa.eu/documents/data-specifications-template):

These Technical Guidelines propose to use the http URIs provided by the Open Geospatial Consortium as coordinate reference system identifiers (see identifiers for the default CRSs below). These are based on and redirect to the definition in the EPSG Geodetic Parameter Registry (http://www.epsg-registry.org/).

The identifier of the coordinate reference system that is applicable for Denmark in this list is http://www.opengis.net/def/crs/EPSG/0/3044, and that has axis order northing, easting.

In GeoPackage, the specification of the coordinate reference system used is not done via HTTP URIs but by means of an entry in a table called gpkg_spatial_ref_sys, which has amongst others has columns for the definition in Well-Known Text, name of the defining organisation and CRS identifier as assigned by that organisation. See also http://www.geopackage.org/spec/#spatial_ref_sys (that is a link to the not yet approved 1.3 version, but it contains a clarification of the previous versions that may be relevant here).

If I understand the implementing rules correct, they actually do not say anything about axis order. Which means that using EPSG:25832 - which is the same as EPSG:3044 but with a reversed axis order - would be compliant to INSPIRE as well: the datum is the European Terrestrial Reference System 1989, the ellipsoid is GRS 1980 and the projection is Transverse Mercator.

It could look as follows in a GeoPackage file (I left out the other CRSs, obligatory entries for the sake of clarity, and not (yet) making use of the GeoPackage's WKT CRS extension):

srs_name	srs_id	organization	organization_coordsys_id	definition	description
ETRS89 / UTM zone 32N	25832	EPSG	25832	PROJCS["ETRS89 / UTM zone 32N",GEOGCS["ETRS89",DATUM["European_Terrestrial_Reference_System_1989",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6258"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4258"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","25832"]]	ETRS89 / UTM zone 32N with axis order easting, northing

I would expect GeoPackage to be used mainly in GIS-systems, and many tools and users in that world expect coordinates and coordinate reference systems to use "traditional GIS-order". E.g. the tool Spatialite, that can work with the GeoPackage format, only supports that axis order. And community behind the tools intends to keep it that way. See https://www.gaia-gis.it/fossil/libspatialite/tktview/96216ad92d29353a9988f3ba1f8ad6f4349760a6 and https://groups.google.com/g/spatialite-users/c/-sbd9F5Batc

Given the amount of years that was needed to implement proper axis order handling in GML in GDAL, QGIS, etc., I think it would be good to

start now with distributing GeoPackages with the easting, northing order, so a short-term and usable solution
focus on improved handling of axis order by tools in the long-term, so that it wouldn't matter anymore what axis order is used in source data. This is already done in e.g. the PROJ tool, used by many other tools, see also https://proj.org/faq.html#why-is-the-axis-ordering-in-proj-not-consistent

PROJ respects the axis ordering as it was defined by the authority in charge of a given coordinate reference system. This is in accordance to the ISO19111 standard [ISO19111]. Unfortunately most GIS software on the market doesn’t follow this standard. Before version 6, PROJ did not respect the standard either. This causes some problems while the rest of the industry conforms to the standard. PROJ intends to spearhead this effort, hopefully setting a good example for the rest of the geospatial industry.

If support for GeoPackage would be added to the INSPIRE validator, we should think about how to check that a valid coordinate reference system is used. As the WKT is provided, an option could be to look at that, instead of only looking at the identifiers. E.g. for 25832:

projinfo -o WKT1:GDAL,WKT2:2019 EPSG:25832
WKT2:2019 string:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        DATUM["European Terrestrial Reference System 1989",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Europe - 6°E to 12°E and ETRS89 by country"],
        BBOX[38.76,6,83.92,12]],
    ID["EPSG",25832]]

WKT1:GDAL string:
PROJCS["ETRS89 / UTM zone 32N",
    GEOGCS["ETRS89",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6258"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4258"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","25832"]]

Note: this issue was originally posted on the INSPIRE forum: https://inspire.ec.europa.eu/forum/discussion/view/268970/coordinate-reference-systems-for-geopackage

heidivanparys commented 3 years ago

An example of what Spatialite does (even though PROJ version 7 is used behind the scenes) for two points that have CRS EPSG:3044:

select hex(location) as location_hex, aswkt(geomfromgpb(location)) location_wkt, asgml(3, geomfromgpb(location)) location_gml from myfeaturetype;

location_hex	location_wkt	location_gml
47500001E40B00000101000000000000009A0926410000000013925741	POINT(722125 6178892)	722125 6178892</gml:pos></gml:Point>
47500001E40B0000010100000000000000A40926410000000013925741	POINT(722130 6178892)	722130 6178892</gml:pos></gml:Point>

As for the hex and WKT output: the de facto standard for WKB and WKT is (x,y{,z}{,m}) where x is easting or longitude, y is northing or latitude, z is optional elevation, and m is optional measure. Therefore, this output is expected and correct. See also http://www.geopackage.org/spec/#gpb_data_blob_format and https://github.com/opengeospatial/geopackage/issues/540

As for the GML output: the coordinate order is wrong: expected <gml:pos>6178892 722125</gml:pos> and <gml:pos>6178892 722130</gml:pos>.

heidivanparys commented 3 years ago

This issue is not only related to GeoPackage. Closing this one, see INSPIRE-MIF/helpdesk-registry#80 instead.

INSPIRE-MIF / gp-geopackage-encodings

Coordinate reference systems #2