Summary

The document describe work related to integration of PROJ 6 with GDAL, which
adds different capabilities: support for CRS WKT 2 version, "late binding"
capabilities for coordinate transformations between CRS, support of time-dimension
for coordinate operations and the use of a unified CRS database.

Motivation

Coordinate systems in GDAL, PROJ, and libgeotiff are missing modern capabilities
and need a thorough refactoring:

The dreaded ad hoc CSV databases in PROJ_LIB and GDAL_DATA are frustrating
for users, pose challenges for developers, and impede interoperability of definitions.

GDAL and PROJ are missing OGC WKT2 support.

PROJ 5.0+ no longer requires datum transformation pivots through WGS84, which
can introduce errors of up to 2m, but the rest of the tools do not take advantage of it.

CSV database

The use of a SQLite-based database for EPSG and other definitions will allow
the projects to add more capability (area-aware validation), transition the
custom peculiar data structures of the projects to something more universally
consumable, and promote definition interoperability between many coordinate
system handling software tools.

WKT2

​OGC WKT2 fixes
longstanding interoperability coordinate system definition discrepancies.
WKT2 contains tools for describing time-dependent coordinate reference systems.
PROJ 5+ is now capable of time-dependent transformations, but GDAL and other
tools do not yet support them.

Several countries are updating their geodetic infrastructure to include
time-dependent coordinate systems. For example, Australia and the United States
are adapting time-dependent coordinate systems in 2020 and 2022, respectively.
The familiar NAD83 and NAVD88 in North America being replaced by NATRF2022 and
NAPGD2022, and the industry WILL have to adapt to these challenges sooner or later.

WGS84 Pivot

PROJ previously required datum transformation that pivoted through WGS84 via a
7-parameter transform. This pivot is a practical solution, but it can introduce
error of about two meters, and many legacy datums cannot be defined in terms of
WGS84. PROJ 5+ now provides the tools to support late-binding through its
​transformation pipeline framework,
but GDAL and the rest of the tools cannot use it yet. Higher accuracy
transformations avoid stepping through WGS84 and eliminates extra transformation
steps with side-car data from a local geodetic authority.

Related work in other libraries

This RFC is the last step in the "gdalbarn" work. Previous steps have consisted
in implementing the related changes in PROJ master per ​PROJ RFC 2
and in libgeotiff master per ​libgeotiff pull request 2.

Proposal

Third-party library requirements

Regarding PROJ, no internal copy of PROJ will be embedded in GDAL master.
It is not doable of supporting older versions of PROJ, as the OGRSpatialReference
class has been largely rewritten to take advantage of functionality
that has been completely moved from GDAL to PROJ: PROJ string import and export,
WKT string import and export, EPSG database exploitation.
To be able to use more easily GDAL master and PROJ master in
complex setups where some GDAL dependencies use a libproj provided by the
system, and where mixing naively PROJ master and this older libproj would result
in runtime crashes, PROJ master can be built with CFLAGS/CXXFLAGS=-DPROJ_RENAME_SYMBOLS
to alias its public symbols, and GDAL will be able to use this custom build.
Note that this is not intended to be used in a long term, since proper packaging
solutions will eventually use PROJ 6 to rebuild all its reverse dependencies.
It should be noted also that PROJ is required at configure / nmake time, that
is the dynamic loading at runtime through dlopen() / LoadLibrary?() is no longer
available.

Regarding libgeotiff, the internal copy in frmts/gtiff/libgeotiff has been
refreshed with the content of upstream libgeotiff master.

All continuous integration systems (Travis-CI and AppVeyor?) have been updated
to build PROJ master as part of the GDAL build.

OGRSpatialReference rewrite

The OGRSpatialReference class is central in GDAL/OGR for all coordinate reference
systems (CRS) manipulations. Up to GDAL 2.4, this class contained mostly a
OGR_SRSNode root node of a WKT 1 representation, and all getters and setters
manipulated this tree representation. As part of this work, the main object
contained internally by OGRSpatialReference is now a PROJ PJ object, and methods
call PROJ C API getters and setters on this PJ object. This enables to be,
mostly (*), representation independent.

WKT1, WKT2, ESRI WKT, PROJ strings import and export is now delegated to PROJ.
The same holds for import of CRS from the EPSG database, that now relies on
proj.db SQLite database. Consequently all the data/*.csv files that contained
CRS related information have been removed from GDAL.
It should be noted that "morphing" from ESRI WKT is now done automatically when
importing WKT.

While general semantics of methods like IsSame?() or FindMatches?() remain the
same, underneath implementations are substantially different, which can lead to
different results than previous GDAL versions in some cases. In the FindMatches?()
case, identification of CRS to EPSG entries is generally improved due to enhanced
query capabilities in the database.

(*) The "mostly" precision is here since it was not
practical to do this rewrite in every place. So for some methods, an internal
WKT1 export is still done. This is the case for methods that take a path to a SRS
node (like "GEOGCS|UNIT") as an argument, or some methods like SetProjection?(),
GetProjParm?(), that expect a OGC WKT1 specific name. Those are thought to be
used mostly be drivers. Changing them to be EPSG names would impact a number of
drivers, some of them little tested regarding SRS support, and which furthermore
mostly support WKT1 representation only.

OGRCoordinateTransformation changes

Since GDAL 2.3 and initial PROJ 5 support, when transforming between two CRS
we still relied on the PROJ.4 string export of the source and target CRS to
create a coordinate operation pipeline. So this limited to "early-binding"
operations, that is using the WGS84 pivot through towgs84 or nadgrids PROJ keywords.
Now PROJ new capabilities to find appropriate coordinate operations between two
CRS is used, offering "late-binding" capabilities to take into account other pivots
than WGS84 or area of uses.

OGRCreateCoordinateOperation() now takes an extra optional arguments to define
options.

One of those options is to define an area of interest that will be
taken into account when searching candidate operations. If several operations
match, the "best" (according to PROJ sorting criterion) will be selected. Note:
it will systematically be used even if later calls to Transform() use coordinates
outside of the initial area of interest.

Another option is the ability to specify the coordinate operation to apply, so as
an override of what GDAL / PROJ would have automatically computed, either
as a PROJ string (generally a +proj=pipeline), or a WKT coordinate
operation/concatenated operation.
Users can typically select a specific coordinate operation by using the new
PROJ projinfo utility that can return the candidate operations from a
source_crs / target_crs tuple.

When no option is specified, GDAL will use PROJ to list all candidate coordinate
operations. For each call to Transform(), it will compute the average coordinate
of the input coordinates and use it to determine the best coordinate operation from
the candidate ones.

The Transform() method now takes an extra argument to contain the coordinate
epoch (generally as a decimal year value) for coordinate operations that are
time-dependent. Related, the transform options of the GDALTransform mechanism
typically used by gdalwarp now accepts a COORDINATE_EPOCH for the same purpose.

Use of OGRSpatialReference in GDAL

Currently GDAL datasets accept and return a WKT 1 string to describe the SRS.
To be more independent of the actual encoding, and for example allowing a
GeoPackage raster dataset to be able to use WKT 2, it is desirable to be
able to attach a SRS that is not dependent of the representation (WKT 1 or WKT 2),
hence using a OGRSpatialReference object instead of a const char* string.

Default WKT version

OGRSpatialReference::exportToWkt() without options will report WKT 1 (with
explicit AXIS nodes. See below "Axis order issues" paragraph) for CRS
compatibles of this representation, and otherwise use WKT2:2018 (typically for
Geographic 3D CRS).

An enhanced version of exportToWkt() accepts options to specify the exact WKT
version used, if multi-line or single-line output must be used, etc.

Alternatively the OSR_WKT_FORMAT configuration option can be used to modify the
WKT version used by exportToWk() (when no explicit version is passed in the
options of exportToWkt())

The gdalinfo, ogrinfo and gdalsrsinfo utililies will default to outputing
WKT2:2018

Axis order issues

This is a recurring pain point. This RFC proposes a new approach (without
pretending to solving it completely) to what was initially done per
RFC 20: OGRSpatialReference Axis Support.
The issue is that CRS official definitions use axis orders that do not conform
to the way raster or vector data is traditionally encoded in GIS applications.
The typical example is the Geographic "WGS 84" definition from EPSG, EPSG:4326,
which uses latitude as the first axis and longitude as the second axis.
RFC 20 decided that by default the AXIS definition would be stripped off from
the WKT when the axis order from the authority did not match the GIS friendly one
(and use a custom EPSGA authority to have WKT with official AXIS elements)

This was technically possible since the WKT 1 grammar makes the AXIS element
definition. However removal of the AXIS definitions was a potential source of
confusion as it was unclear which axis order was actually used. Furthermore,
in WKT2, the AXIS element is compulsory, and the internal PROJ representation
requires also a coordinate system to be defined. So there would have been
two unsatisfactory options:

return patched versions of the official definition with the GIS friendly order,
while still using the official authority code. Practical since we keep the
link with the source code, but a lie since we modify it. Users would not know
whether they must trust the encoded order, or the official order from the
authority.

return patched versions of the official definition with the GIS friendly order,
but without the official authority code. This would be compliant, but we would
lose the link with the authority code.

The solution put forward in this RFC is 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

OAMS_TRADITIONAL_GIS_ORDER or
OAMS_AUTHORITY_COMPLIANT or
OAMS_CUSTOM )

is added to make their job of specifying the axis mapping easier;

OAMS_TRADITIONAL_GIS_ORDER means:

for geographic 2D CRS,

for Latitude NORTH, Longitude EAST (such as EPSG:4326), GetDataAxisToSRSAxisMapping() returns {2,1}, meaning that the data order is longitude, latitude

for Longitude EAST, Latitude NORTH (such as OGC:CRS84), returns {1,2}

for projected CRS,

for EAST, NORTH (ie most projected CRS), return {1,2}

for NORTH, EAST, 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

for all other cases, return {1,2}

OGRCreateCoordinateTransformation() now honors the data axis to srs axis mapping.

Note: contrary to what I indicated in a previous email, gdaltransform behavior
is unchanged, since internally the GDALTransform mechanism forces the GIS friendly
order.

Raster datasets are modified to call SetAxisMappingStrategy?(OAMS_TRADITIONAL_GIS_ORDER) on the
OGRSpatialReference* they return, and assumes it in SetSpatialRef?()
(assumed and unchecked for now)

Vector layers mostly all call SetAxisMappingStrategy?(OAMS_TRADITIONAL_GIS_ORDER) on the
OGRSpatialReference* returned by GetSpatialRef?(). In the case of the GML driver,
if the user defines the INVERT_AXIS_ORDER_IF_LAT_LONG open option, axis
swapping is not done (as previously) and the AUTHORITY_COMPLIANT strategy is
used.
ICreateLayer() when receiving a OGRSpatialReference* may decide (and most will
do it) 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. ogr2ogr is
modified to do the geometry axis swapping in that case.

Related to that change, WKT 1 export now always return the AXIS element,
and EPSG:xxxx thus behaves identically to 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.

Driver changes

Raster drivers that returned / accepted a SRS as a WKT string through the
GetProjectionRef?(), SetProjection?(), GetGCPProjection() and SetGCPs() methods
have been upgraded to use the new virtual methods, in most cases by using the
compatibility layer.

The GDALPamDataset (PAM .aux.xml files) and the GDAL VRT drivers have been
fully upgraded to support the new interfaces, and serialize/deserialize the data
axis to SRS axis mapping values.

The GeoPackage driver now fully supports the official "gpkg_crs_wkt" extension
used to store WKT 2 string definitions in the gpkg_spatial_ref_sys table. The
driver attempts at not using the extension when SRS can be encoded as WKT1 strings,
and will automatically add the "definition_12_063" column to an existing
gpkg_spatial_ref_sys table if a SRS requiring WKT2 (typically a Geographic 3D CRS)
is inserted.

Changes in utilities

gdalinfo and ogrinfo reports the data axis to CRS axis mapping whenever a
CRS is reported. They will also output WKT2_2018 by default, unless "-wkt_format wkt1" is specified.

gdalwarp, ogr2ogr and gdaltransform have gained a -ct switch that can be used
by advanced users to specify a coordinate operation, either as a PROJ string
(generally a +proj=pipeline), or a WKT coordinate operation/concatenated operation,
as explained in the above "OGRCoordinateTransformation changes" paragraph.
Note: the pipeline must take into account the axis order of the CRS, even if the
underlying raster/vector drivers use the "GIS friendly" order.
For example "+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=utm +zone=31 +ellps=WGS84"
when transforming from EPSG:4326 to EPSG:32631.

gdalsrsinfo is enhanced to be able to specify the 2 new supported WKT variants: WKT2_2015 and WKT2_2018.
It will default to outputing WKT2_2018

SWIG binding changes

The enhanced ExportToWkt?() and OGRCoordinateTransformation methods are available
through SWIG bindings. May require additional typemaps for non-Python languages
(particularly for the support of 4D X,Y,Z,time coordinates)

Backward compatibility

This work is intended to be *mostly* backward compatible, yet inevitable differences
will be found. For example the WKT 1 and PROJ string export has been completely
rewritten in PROJ, and so while being hopefully equivalent to what GDAL 2.4 or
earlier generated, this is not strictly identical: number of significant digits,
order of PROJ parameters, rounding, etc etc...

MIGRATION_GUIDE.TXT has been updated to reflect some differences:

OSRImportFromEPSG() takes into account official axis order.

removal of OPTGetProjectionMethods(), OPTGetParameterList() and OPTGetParameterInfo()
No equivalent.

removal of OSRFixup() and OSRFixupOrdering(): no longer needed since objects
constructed are always valid

removal of OSRStripCTParms(). Use OSRExportToWktEx() instead with the
FORMAT=SQSQL option

exportToWkt() outputs AXIS nodes

OSRIsSame(): now takes into account data axis to CRS axis mapping, unless
IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES is set as an option to OSRIsSameEx()

ogr_srs_api.h: SRS_WKT_WGS84 macro is no longer declared by default since
WKT without AXIS is too ambiguous. Preferred remediation: use SRS_WKT_WGS84_LAT_LONG.
Or #define USE_DEPRECATED_SRS_WKT_WGS84 before including ogr_srs_api.h

Out-of-tree raster drivers will be impacted by the introduction of the
new virtual methods GetSpatialRef?(), SetSpatialRef?(), GetGCPSpatialRef() and
SetGCPs(..., const OGRSpatialReference* poSRS), and the removal of their older
equivalents using WKT strings instead of a OGRSpatialReference* instance.

Documentation

New methods have been documented, and documentation of existing methods has
been changed when appropriate during the development. That said, a more
thorough pass will be needed. The tutorials will also have to be updated.

Testing

The autotest suite has been adapted in a number of places since the expected
results have changed for a number of reasons (AXIS node exported in WKT,
differences in WKT and PROJ string generation). New tests have been added for the
new capabilities.

It should be noted that autotest not necessarily checks everything, and issues
have been discovered and fixed through manual testing. The introduction of the
"data axis to CRS axis mapping" concept is also quite error prone, as it requires
setting the OAMS_TRADITIONAL_GIS_ORDER strategy in a lot of different places.

So users and developers are kindly invited to thoroughly test GDAL once this
work has landed in master.

Implementation

While it is provided as a multiple commit for """easier""" review, it will be
probably squashed in a single commit for inclusion in master, as intermediate steps
are not all buildable, due to PROJ symbol renames having occured during the
development, which would break bisectability.