[gdal-dev] Pre-RFC for OGRSpatialReference use in GDAL and changes in how to deal with axis order issues

[Even Rouault]

Hi,

I call this pre-RFC as I don't have yet any code implementing the following 
ideas (so there might be holes in the analysis), but I wanted to get some 
feedback before investing too much time into this substantial work.

Let's start with the easiest topic:

1) Use of OGRSpatialReference in GDAL

Currently GDAL datasets accept and return a WKT 1 string to describe the SRS. 
To be more independant of the actual encoding, and for example allowing a 
GeoPackage raster dataset to be able to use WKT 2, I propose we switch from a 
string representation to the more abstract representation offered by 
OGRSpatialReference. In the current state of GDAL, OGRSpatialReference is just 
a direct mapping of WKT 1, but in my in-progress gdalbarn branch at
https://github.com/rouault/gdal/tree/gdalbarn , OGRSpatialReference is now 
more or less encapsulating a osgeo::proj::crs object (ie a C++ object 
representing a CRS as modelled by OGC Topic 2 / ISO-19111), that can have 
several representations: OGC WKT1 (GDAL WKT1), WKT2, ESRI WKT1, PROJ strings.

So the proposal would be to change:

GDALDataset:
    virtual const char* GetProjectionRef();

to
    virtual OGRSpatialReference *GetSpatialRef();

    // compatibility layer (conceptually)
    const char* GetProjectionRef() {
		return GetSpatialRef()->ExportToWkt1(); }

Similarly on the set side,

Change
    virtual CPLErr SetProjection(const char* pszWKT);
to
    virtual CPLErr SetSpatialRef(OGRSpatialReference*);

    // compatibility layer (conceptually)
    CPLErr SetProjection(const char* pszWKT) {
		return SetSpatialRef(OGRSpatialReference::importFromWkt(pszWKT)); }

Regarding ownership, I guess we would do similarly to the vector side, that is 
GetSpatialRef() returns an an object that it owns and doesn't increment the 
refcount. And SetSpatialRef() would potentially increment the reference count 
on the passed object or clone it (like ICreateLayer() does)

To ease the transition that will require going through all the ~150 raster 
drivers, I may add in GDALDataset:
  protected:
    OGRSpatialReference* GetSpatialRefFromOldGetProjectionRef()
		// conceptual implementation
		{ return OGRSpatialReference::importFromWkt(OldGetProjectionRef()); }

    virtual const char* OldGetProjectionRef() { return ""; }

Drivers would rename their existing GetProjectionRef() to 
OldGetProjectionRef() and implement a GetSpatialRef() that would just call 
GetSpatialRefFromOldGetProjectionRef(). This is just a technical detail to be 
able to complete the transition of API easier in a semi-automatic way. Over 
time we may remove the hack and implement natively GetSpatialRef(). And 
similarly for the write side.


Now, the harder one:

2) Axis order revamp

Ever recurring headache [1], not sure my below proposal will solve it in a 
definitive way

Let's look at the current situation with a common example

EPSG:4326 returns
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

whereas EPSGA:4326 returns (not the A for Axis)
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4326"]]

In EPSG:4326 expansion, there is no AXIS mention. This is technically 
conformant with OGC WKT1 grammar where AXIS is optional. However, GDAL uses 
this omission to implicitely mean it doesn't honour the axis order that the 
EPSG authority has set for this geographic CRS, which is latitude first, 
longitude second, but instead use longitude, latitude. Note also that EPSGA:
4326 is not really honoured by GDAL: if you use it with gdaltransform, it will 
still expect/return coordinates in longitude, latitude order.

In WKT2, AXIS is a required element, so we cannot hide ourselves behind the 
omission we did in WKT1, and we can't either modify the definition to 
advertize longitude, latitude order, as it would be a great source of 
confusion (using an authority code but putting a contradictory element into it 
is a no-no)

I propose to add a "data axis to SRS axis mapping" concept, which is a bit 
similar to what is done in WCS DescribeCoverage response to explain how the 
SRS axis map to the grid axis of a coverage

Extract from
https://docs.geoserver.org/stable/en/user/extensions/wcs20eo/index.html
for a coverage that uses EPSG:4326
	 <gml:coverageFunction>
        <gml:GridFunction>
          <gml:sequenceRule axisOrder="+2 +1">Linear</gml:sequenceRule>
          <gml:startPoint>0 0</gml:startPoint>
        </gml:GridFunction>
      </gml:coverageFunction>

A similar mapping would be used to define how the 'x' and 'y' components in 
the geotransform matrix or in a OGRGeometry map to the axis defined by the CRS 
definition.

Such mapping would be given by a new method in OGRSpatialReference

    std::vector<int> GetDataAxisToSRSAxisMapping() const

To explain its semantics, imagine that it return 2,-1,3.
That would be interpreted as:
* 2: the first axis of the CRS maps to the second axis of the data
* -1: the second axis of the CRS maps to the first axis of the data, with 
values negated
* 3: the third axis of the CRS maps to the third axis of the data

This would be similar to the PROJ axisswap operation:
	https://proj4.org/operations/conversions/axisswap.html

By default, GetDataAxisToSRSAxisMapping() would return the identity 1,2[,3], 
that is be conformant with the axis order defined by the authority.

As all GDAL and a vast majority of OGR drivers depend on using the "GIS axis 
mapping", a method
	SetDataAxisToSRSAxisMappingStrategy(
			TRADITIONAL_GIS_ORDER or AUTHORITY_COMPLIANT )
would be added to make their job of specifying the axis mapping easier;

TRADITIONAL_GIS_ORDER would mean:
* for geographic 2D CRS,
    - for Latitude NORTH, Longitude EAST (such as EPSG:4326), 
GetDataAxisToSRSAxisMapping() would return {2,1}, meaning that the data order 
is longitude, latitude
    - for Longitude EAST, Latitude NORTH (such as OGC:CRS84), would return 
{1,2}

* for projected CRS,
    - for EAST, NORTH (ie most projected CRS), would return {1,2}
    - for NORTH, EAST, would return {2,1}
    - for North Pole CRS, with East/SOUTH, North/SOUTH, such as EPSG:5041 
("WGS 84 / UPS North (E,N)"), would return {1,2}
    - for North Pole CRS, with northing/SOUTH, easting/SOUTH, such as EPSG:
32661 ("WGS 84 / UPS North (N,E)"), would return {2,1}
    - similarly for South Pole CRS
    - Krovak East North, with EAST, NORTH, such as EPSG:5221 "S-JTSK (Ferro) / 
Krovak", would return {1,2}

	Cases where I'm less sure about what the 'GIS order' is and suspect we 
might not do the right thing currently:
    - Transverse Mercator South Oriented, with WEST, SOUTH, would return {1,2}   
(I assume that this system is not used for rasters)
    - traditional Krovak, with SOUTH, WEST, such as EPSG:2065 "S-JTSK (Ferro) 
/ Krovak", would return {1,2}   (I assume that this system is not used for 
rasters)
    - other cases, return {1,2}

OGRCreateCoordinateTransformation(OGRSpatialReference* source, 
OGRSpatialReference* target) will honour the data axis to srs axis mapping
That means that gdaltransform -s_srs EPSG:4326 -t_srs EPSG:32631 will expect 
latitude, longitude order (so backward incompatible)
(PROJ6 will also be axis order compliant when using EPSG codes)
We may add a -s_axis gis switch (and/or -s_axis long,lat for more explicit) / 
-t_axis syntax to have some easier transition.

Raster datasets will all set 
SetDataAxisToSRSAxisMappingStrategy(TRADITIONAL_GIS_ORDER) on the 
OGRSpatialReference* they return, and will also assume it in SetSpatialRef() 
(will probably be just assumed and unchecked for now)

Vector layers will mostly all call 
SetDataAxisToSRSAxisMappingStrategy(TRADITIONAL_GIS_ORDER) on the 
OGRSpatialReference* returned by GetSpatialRef(). We may decide that some 
drivers like GML no longer do coordinate axis swapping for GML3 and then use 
the AUTHORITY_COMPLIANT strategy.
ICreateLayer() when receiving a OGRSpatialReference* may decide to change the 
axis mapping strategy. That is: if it receives a OGRSpatialReference with 
AUTHORITY_COMPLIANT order, it may decide to switch to TRADITIONAL_GIS_ORDER 
and GetSpatialRef()::GetDataAxisToSRSAxisMapping() will reflect that. Tools 
like ogr2ogr will need to do the geometry axis swapping in that case

I'm midly satisfied with adding the axis mapping in the OGRSpatialReference 
class itself, but I feel it would be overkill to have an extra object to 
capture both (and can't come with a name for it). The coupling of both will 
require careful explanation of what some methods do. For example, I think 
OGRSpatialReference::IsSame() would ignore the data axis mapping setting when 
comparing 2 OGRSpatialReference object.

I would also lean on making WKT 1 export now always return the AXIS element, 
and EPSG:xxxx would then behave like EPSGA:xxxx

So a summary view of this approach is that in the formal SRS definition, we no 
longer do derogations regarding axis order, but we add an additional interface 
to describe how we actually make our match match with the SRS definition.

In case you reached that point, thoughts... ?

Even

[1] Fresh ticket about such confusion:
https://github.com/OSGeo/gdal/issues/1155

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com