Changes between Version 1 and Version 2 of rfc73_proj6_wkt2_srsbarn


Ignore:
Timestamp:
Jan 9, 2019, 7:19:42 AM (5 years ago)
Author:
Even Rouault
Comment:

Initial version of RFC 73

Legend:

Unmodified
Added
Removed
Modified
  • rfc73_proj6_wkt2_srsbarn

    v1 v2  
    11= RFC 73: Integration of PROJ6 for WKT2, late binding capabilities, time-support and unified CRS database =
     2
     3|| Author:  || Even Rouault ||
     4|| Contact: || even.rouault @ spatialys.com ||
     5|| Started: || 2018-Jan-08 ||
     6|| Status:  || Development, integration target: GDAL 2.5 ||
     7
     8== Summary ==
     9
     10The document describe work related to integration of PROJ 6 with GDAL, which
     11adds different capabilities: support for CRS WKT 2 version, "late binding"
     12capabilities for coordinate transformations between CRS, support of time-dimension
     13for coordinate operations and the use of a unified CRS database.
     14
     15== Motivation ==
     16
     17The motivations are those exposed in https://gdalbarn.com/#why , which are
     18copied here
     19
     20Coordinate systems in GDAL, PROJ, and libgeotiff are missing modern capabilities
     21and need a thorough refactoring:
     22* The dreaded ad hoc CSV databases in PROJ_LIB and GDAL_DATA are frustrating
     23  for users, pose challenges for developers, and impede interoperability of definitions.
     24* GDAL and PROJ are missing OGC WKT2 support.
     25* PROJ 5.0+ no longer requires datum transformation pivots through WGS84, which
     26  can introduce errors of up to 2m, but the rest of the tools do not take advantage of it.
     27
     28=== CSV database ===
     29
     30The use of a SQLite-based database for EPSG and other definitions will allow
     31the projects to add more capability (area-aware validation), transition the
     32custom peculiar data structures of the projects to something more universally
     33consumable, and promote definition interoperability between many coordinate
     34system handling software tools.
     35
     36=== WKT2 ===
     37
     38[http://docs.opengeospatial.org/is/12-063r5/12-063r5.html OGC WKT2] fixes
     39longstanding interoperability coordinate system definition discrepancies.
     40WKT2 contains tools for describing time-dependent coordinate reference systems.
     41PROJ 5+ is now capable of time-dependent transformations, but GDAL and other
     42tools do not yet support them.
     43
     44Several countries are updating their geodetic infrastructure to include
     45time-dependent coordinate systems. For example, Australia and the United States
     46are adapting time-dependent coordinate systems in 2020 and 2022, respectively.
     47The familiar NAD83 and NAVD88 in North America being replaced by NATRF2022 and
     48NAPGD2022, and the industry WILL have to adapt to these challenges sooner or later.
     49
     50=== WGS84 Pivot ===
     51
     52PROJ previously required datum transformation that pivoted through WGS84 via a
     537-parameter transform. This pivot is a practical solution, but it can introduce
     54error of about two meters, and many legacy datums cannot be defined in terms of
     55WGS84. PROJ 5+ now provides the tools to support late-binding through its
     56[https://proj4.org/usage/transformation.html#geodetic-transformation transformation pipeline framework],
     57but GDAL and the rest of the tools cannot use it yet. Higher accuracy
     58transformations avoid stepping through WGS84 and eliminates extra transformation
     59steps with side-car data from a local geodetic authority.
     60
     61== Related work in other libraries ==
     62
     63This RFC is the last step in the "gdalbarn" work. Previous steps have consisted
     64in implementing the related changes in PROJ master per [https://proj4.org/community/rfc/rfc-2.html PROJ RFC 2]
     65and in libgeotiff master per [https://github.com/OSGeo/libgeotiff/pull/2 libgeotiff pull request 2].
     66
     67== Proposal ==
     68
     69=== Third-party library requirements ===
     70
     71GDAL master (future 2.5.0) will require PROJ master (future PROJ 6.0) and
     72libgeotiff master (future libgeotiff 1.5 or 2.0) for build and execution.
     73
     74Regarding PROJ, no internal copy of PROJ will be embedded in GDAL master.
     75It is not doable of supporting older versions of PROJ, as the OGRSpatialReference
     76class has been largely rewritten to take advantage of functionality
     77that has been completely moved from GDAL to PROJ: PROJ string import and export,
     78WKT string import and export, EPSG database exploitation.
     79To be able to use more easily GDAL master and PROJ master in
     80complex setups where some GDAL dependencies use a libproj provided by the
     81system, and where mixing naively PROJ master and this older libproj would result
     82in runtime crashes, PROJ master can be built with CFLAGS/CXXFLAGS=-DPROJ_RENAME_SYMBOLS
     83to alias its public symbols, and GDAL will be able to use this custom build.
     84Note that this is not intended to be used in a long term, since proper packaging
     85solutions will eventually use PROJ 6 to rebuild all its reverse dependencies.
     86It should be noted also that PROJ is required at configure / nmake time, that
     87is the dynamic loading at runtime through dlopen() / LoadLibrary() is no longer
     88available.
     89
     90Regarding libgeotiff, the internal copy in frmts/gtiff/libgeotiff has been
     91refreshed with the content of upstream libgeotiff master.
     92
     93All continuous integration systems (Travis-CI and AppVeyor) have been updated
     94to build PROJ master as part of the GDAL build.
     95
     96=== OGRSpatialReference rewrite ===
     97
     98The OGRSpatialReference class is central in GDAL/OGR for all coordinate reference
     99systems (CRS) manipulations. Up to GDAL 2.4, this class contained mostly a
     100OGR_SRSNode root node of a WKT 1 representation, and all getters and setters
     101manipulated this tree representation. As part of this work, the main object
     102contained internally by OGRSpatialReference is now a PROJ PJ object, and methods
     103call PROJ C API getters and setters on this PJ object. This enables to be,
     104mostly (*), representation independent.
     105
     106WKT1, WKT2, ESRI WKT, PROJ strings import and export is now delegated to PROJ.
     107The same holds for import of CRS from the EPSG database, that now relies on
     108proj.db SQLite database. Consequently all the data/*.csv files that contained
     109CRS related information have been removed from GDAL.
     110It should be noted that "morphing" from ESRI WKT is now done automatically when
     111importing WKT.
     112
     113While general semantics of methods like IsSame() or FindMatches() remain the
     114same, underneath implementations are substantially different, which can lead to
     115different results than previous GDAL versions in some cases. In the FindMatches()
     116case, identification of CRS to EPSG entries is generally improved due to enhanced
     117query capabilities in the database.
     118
     119(*) The "mostly" precision is here since it was not
     120practical to do this rewrite in every place. So for some methods, an internal
     121WKT1 export is still done. This is the case for methods that take a path to a SRS
     122node (like "GEOGCS|UNIT") as an argument, or some methods like SetProjection(),
     123GetProjParm(), that expect a OGC WKT1 specific name. Those are thought to be
     124used mostly be drivers. Changing them to be EPSG names would impact a number of
     125drivers, some of them little tested regarding SRS support, and which furthermore
     126mostly support WKT1 representation only.
     127
     128=== OGRCoordinateTransformation changes ===
     129
     130Since GDAL 2.3 and initial PROJ 5 support, when transforming between two CRS
     131we still relied on the PROJ.4 string export of the source and target CRS to
     132create a coordinate operation pipeline. So this limited to "early-binding"
     133operations, that is using the WGS84 pivot through towgs84 or nadgrids PROJ keywords.
     134Now PROJ new capabilities to find appropriate coordinate operations between two
     135CRS is used, offering "late-binding" capabilities to take into account other pivots
     136than WGS84 or area of uses.
     137
     138OGRCreateCoordinateOperation() now takes an extra optional arguments to define
     139options.
     140
     141One of those options is to define an area of interest that will be
     142taken into account when searching candidate operations. If several operations
     143match, the "best" (according to PROJ sorting criterion) will be selected. Note:
     144it will systematically be used even if later calls to Transform() use coordinates
     145outside of the initial area of interest.
     146
     147Another option is the ability to specify the coordinate operation to apply, so as
     148an override of what GDAL / PROJ would have automatically computed, either
     149as a PROJ string (generally a +proj=pipeline), or a WKT coordinate
     150operation/concatenated operation.
     151Users can typically select a specific coordinate operation by using the new
     152PROJ projinfo utility that can return the candidate operations from a
     153source_crs / target_crs tuple.
     154
     155When no option is specified, GDAL will use PROJ to list all candidate coordinate
     156operations. For each call to Transform(), it will compute the average coordinate
     157of the input coordinates and use it to determine the best coordinate operation from
     158the candidate ones.
     159
     160The Transform() method now takes an extra argument to contain the coordinate
     161epoch (generally as a decimal year value) for coordinate operations that are
     162time-dependent. Related, the transform options of the GDALTransform mechanism
     163typically used by gdalwarp now accepts a COORDINATE_EPOCH for the same purpose.
     164
     165=== Use of OGRSpatialReference in GDAL ===
     166
     167Currently GDAL datasets accept and return a WKT 1 string to describe the SRS.
     168To be more independent of the actual encoding, and for example allowing a
     169GeoPackage raster dataset to be able to use WKT 2, it is desirable to be
     170able to attach a SRS that is not dependent of the representation (WKT 1 or WKT 2),
     171hence using a OGRSpatialReference object instead of a const char* string.
     172
     173The following new methods are added in GDALDataset:
     174* virtual const OGRSpatialReference* GetSpatialRef() const;
     175* virtual CPLErr SetSpatialRef(const OGRSpatialReference*);
     176* virtual const OGRSpatialReference* GetGCPSpatialRef() const;
     177* virtual CPLErr SetGCPs(int nGCPCount, const GDAL_GCP *pasGCPList, const OGRSpatialReference*);
     178
     179To ease the transition, the following non virtual methods are added in
     180GDALDataset:
     181* const OGRSpatialReference* GetSpatialRefFromOldGetProjectionRef() const;
     182* CPLErr OldSetProjectionFromSetSpatialRef(const OGRSpatialReference* poSRS);
     183* const OGRSpatialReference* GetGCPSpatialRefFromOldGetGCPProjection() const;
     184* CPLErr OldSetGCPsFromNew( int nGCPCount, const GDAL_GCP *pasGCPList,
     185                              const OGRSpatialReference * poGCP_SRS );
     186
     187and the previous GetProjectionRef(), SetProjection(), GetGCPProjection() and SetGCPs()
     188are available as projected virtual methods, prefixed by an underscore
     189
     190This way to convert an existing driver, it is a matter of renaming its
     191GetProjectionRef() method as _GetProjectionRef(), and adding:
     192{{{
     193const OGRSpatialReference* GetSpatialRef() const override {
     194    return GetSpatialRefFromOldGetProjectionRef();
     195}
     196}}}
     197
     198=== Default WKT version ===
     199
     200OGRSpatialReference::exportToWkt() without options will report WKT 1 (with
     201explicit AXIS nodes. See below "Axis order issues" paragraph) for CRS
     202compatibles of this representation, and otherwise use WKT 2:2015 (typically for
     203Geographic 3D CRS)
     204
     205An enhanced version of exportToWkt() accepts options to specify the exact WKT
     206version used, if multi-line or single-line output must be used, etc.
     207
     208=== Axis order issues ===
     209
     210This is a recurring pain point. This RFC proposes a new approach (without
     211pretending to solving it completely) to what was initially done per
     212[rfc20_srs_axes RFC 20: OGRSpatialReference Axis Support].
     213The issue is that CRS official definitions use axis orders that do not conform
     214to the way raster or vector data is traditionally encoded in GIS applications.
     215The typical example is the Geographic "WGS 84" definition from EPSG, EPSG:4326,
     216which uses latitude as the first axis and longitude as the second axis.
     217RFC 20 decided that by default the AXIS definition would be stripped off from
     218the WKT when the axis order from the authority did not match the GIS friendly one
     219(and use a custom EPSGA authority to have WKT with official AXIS elements)
     220
     221This was technically possible since the WKT 1 grammar makes the AXIS element
     222definition. However removal of the AXIS definitions was a potential source of
     223confusion as it was unclear which axis order was actually used. Furthermore,
     224in WKT2, the AXIS element is compulsory, and the internal PROJ representation
     225requires also a coordinate system to be defined. So there would have been
     226two unsatisfactory options:
     227- return patched versions of the official definition with the GIS friendly order,
     228  while still using the official authority code. Practical since we keep the
     229  link with the source code, but a lie since we modify it. Users would not know
     230  whether they must trust the encoded order, or the official order from the
     231  authority.
     232- return patched versions of the official definition with the GIS friendly order,
     233  but without the official authority code. This would be compliant, but we would
     234  lose the link with the authority code.
     235
     236The solution put forward in this RFC is to add a "data axis to SRS axis mapping"
     237concept, which is a bit similar to what is done in WCS DescribeCoverage response
     238to explain how the
     239SRS axis map to the grid axis of a coverage
     240
     241Extract from
     242https://docs.geoserver.org/stable/en/user/extensions/wcs20eo/index.html
     243for a coverage that uses EPSG:4326
     244{{{
     245      <gml:coverageFunction>
     246        <gml:GridFunction>
     247          <gml:sequenceRule axisOrder="+2 +1">Linear</gml:sequenceRule>
     248          <gml:startPoint>0 0</gml:startPoint>
     249        </gml:GridFunction>
     250      </gml:coverageFunction>
     251}}}
     252
     253A similar mapping is added to define how the 'x' and 'y' components in
     254the geotransform matrix or in a OGRGeometry map to the axis defined by the CRS
     255definition.
     256
     257Such mapping is given by a new method in OGRSpatialReference
     258
     259    const std::vector<int>& GetDataAxisToSRSAxisMapping() const
     260
     261To explain its semantics, imagine that it return 2,-1,3.
     262That is interpreted as:
     263* 2: the first axis of the CRS maps to the second axis of the data
     264* -1: the second axis of the CRS maps to the first axis of the data, with values negated
     265* 3: the third axis of the CRS maps to the third axis of the data
     266
     267This is similar to the PROJ axisswap operation:
     268        https://proj4.org/operations/conversions/axisswap.html
     269
     270By default, on a newly create OGRSpatialReference object,
     271GetDataAxisToSRSAxisMapping() returns the identity 1,2[,3],
     272that is, conform to the axis order defined by the authority.
     273
     274As all GDAL and a vast majority of OGR drivers depend on using the "GIS axis
     275mapping", a method
     276        SetAxisMappingStrategy(
     277            OAMS_TRADITIONAL_GIS_ORDER or
     278            OAMS_AUTHORITY_COMPLIANT or
     279            OAMS_CUSTOM )
     280is added to make their job of specifying the axis mapping easier;
     281
     282OAMS_TRADITIONAL_GIS_ORDER means:
     283* for geographic 2D CRS,
     284  - for Latitude NORTH, Longitude EAST (such as EPSG:4326), GetDataAxisToSRSAxisMapping() returns {2,1}, meaning that the data order is longitude, latitude
     285  - for Longitude EAST, Latitude NORTH (such as OGC:CRS84), returns {1,2}
     286
     287* for projected CRS,
     288  - for EAST, NORTH (ie most projected CRS), return {1,2}
     289  - for NORTH, EAST, return {2,1}
     290  - for North Pole CRS, with East/SOUTH, North/SOUTH, such as EPSG:5041 ("WGS 84 / UPS North (E,N)"), would return {1,2}
     291  - for North Pole CRS, with northing/SOUTH, easting/SOUTH, such as EPSG:32661 ("WGS 84 / UPS North (N,E)"), would return {2,1}
     292  - similarly for South Pole CRS
     293  - for all other cases, return {1,2}
     294
     295OGRCreateCoordinateTransformation() now honors the data axis to srs axis mapping.
     296
     297Note: contrary to what I indicated in a previous email, gdaltransform behavior
     298is unchanged, since internally the GDALTransform mechanism forces the GIS friendly
     299order.
     300
     301Raster datasets are modified to call SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) on the
     302OGRSpatialReference* they return, and assumes it in SetSpatialRef()
     303(assumed and unchecked for now)
     304
     305Vector layers mostly all call SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER) on the
     306OGRSpatialReference* returned by GetSpatialRef(). In the case of the GML driver,
     307if the user defines the INVERT_AXIS_ORDER_IF_LAT_LONG open option, axis
     308swapping is not done (as previously) and the AUTHORITY_COMPLIANT strategy is
     309used.
     310ICreateLayer() when receiving a OGRSpatialReference* may decide (and most will
     311do it) to change the  axis mapping strategy. That is: if it receives a
     312OGRSpatialReference with
     313AUTHORITY_COMPLIANT order, it may decide to switch to TRADITIONAL_GIS_ORDER
     314and GetSpatialRef()::GetDataAxisToSRSAxisMapping() will reflect that. ogr2ogr is
     315modified to do the geometry axis swapping in that case.
     316
     317Related to that change, WKT 1 export now always return the AXIS element,
     318and EPSG:xxxx thus behaves identically to EPSGA:xxxx
     319
     320So a summary view of this approach is that in the formal SRS definition, we no
     321longer do derogations regarding axis order, but we add an additional interface
     322to describe how we actually make our match match with the SRS definition.
     323
     324=== Driver changes ===
     325
     326Raster drivers that returned / accepted a SRS as a WKT string through the
     327GetProjectionRef(), SetProjection(), GetGCPProjection() and SetGCPs() methods
     328have been upgraded to use the new virtual methods, in most cases by using the
     329compatibility layer.
     330
     331The GDALPamDataset (PAM .aux.xml files) and the GDAL VRT drivers have been
     332fully upgraded to support the new interfaces, and serialize/deserialize the data
     333axis to SRS axis mapping values.
     334
     335The GeoPackage driver now fully supports the official "gpkg_crs_wkt" extension
     336used to store WKT 2 string definitions in the gpkg_spatial_ref_sys table. The
     337driver attempts at not using the extension when SRS can be encoded as WKT1 strings,
     338and will automatically add the "definition_12_063" column to an existing
     339gpkg_spatial_ref_sys table if a SRS requiring WKT2 (typically a Geographic 3D CRS)
     340is inserted.
     341
     342=== Changes in utilities ===
     343
     344* gdalinfo and ogrinfo reports the data axis to CRS axis mapping whenever a
     345  CRS is reported. For example:
     346
     347{{{
     348Driver: GTiff/GeoTIFF
     349Files: out.tif
     350Size is 20, 20
     351Coordinate System is:
     352GEOGCS["WGS 84",
     353    DATUM["WGS_1984",
     354        SPHEROID["WGS 84",6378137,298.257223563,
     355            AUTHORITY["EPSG","7030"]],
     356        AUTHORITY["EPSG","6326"]],
     357    PRIMEM["Greenwich",0],
     358    UNIT["degree",0.0174532925199433,
     359        AUTHORITY["EPSG","9122"]],
     360    AXIS["Latitude",NORTH],
     361    AXIS["Longitude",EAST],
     362    AUTHORITY["EPSG","4326"]]
     363Data axis to CRS axis mapping: 2,1 <-- here
     364Origin = (2.000000000000000,49.000000000000000)
     365Pixel Size = (0.100000000000000,-0.100000000000000)
     366}}}
     367
     368* gdalwarp, ogr2ogr and gdaltransform have gained a -ct switch that can be used
     369  by advanced users to specify a coordinate operation, either as a PROJ string
     370  (generally a +proj=pipeline), or a WKT coordinate operation/concatenated operation,
     371  as explained in the above "OGRCoordinateTransformation changes" paragraph.
     372  Note: the pipeline must take into account the axis order of the CRS, even if the
     373  underlying raster/vector drivers use the "GIS friendly" order.
     374  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"
     375  when transforming from EPSG:4326 to EPSG:32631.
     376
     377* gdalsrsinfo is enhanced to be able to specify the 2 new supported WKT variants: WKT2_2015 and WKT2_2018
     378
     379=== SWIG binding changes ===
     380
     381The enhanced ExportToWkt() and OGRCoordinateTransformation methods are available
     382through SWIG bindings. May require additional typemaps for non-Python languages
     383(particularly for the support of 4D X,Y,Z,time coordinates)
     384
     385== Backward compatibility ==
     386
     387This work is intended to be *mostly* backward compatible, yet inevitable differences
     388will be found. For example the WKT 1 and PROJ string export has been completely
     389rewritten in PROJ, and so while being hopefully equivalent to what GDAL 2.4 or
     390earlier generated, this is not strictly identical: number of significant digits,
     391order of PROJ parameters, rounding, etc etc...
     392
     393MIGRATION_GUIDE.TXT has been updated to reflect some differences:
     394- OSRImportFromEPSG() takes into account official axis order.
     395- removal of OPTGetProjectionMethods(), OPTGetParameterList() and OPTGetParameterInfo()
     396  No equivalent.
     397- removal of OSRFixup() and OSRFixupOrdering(): no longer needed since objects
     398  constructed are always valid
     399- removal of OSRStripCTParms(). Use OSRExportToWktEx() instead with the
     400  FORMAT=SQSQL option
     401- exportToWkt() outputs AXIS nodes
     402- OSRIsSame(): now takes into account data axis to CRS axis mapping, unless
     403  IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES is set as an option to OSRIsSameEx()
     404- ogr_srs_api.h: SRS_WKT_WGS84 macro is no longer declared by default since
     405  WKT without AXIS is too ambiguous. Preferred remediation: use SRS_WKT_WGS84_LAT_LONG.
     406  Or #define USE_DEPRECATED_SRS_WKT_WGS84 before including ogr_srs_api.h
     407
     408Out-of-tree raster drivers will be impacted by the introduction of the
     409new virtual methods GetSpatialRef(), SetSpatialRef(), GetGCPSpatialRef() and
     410SetGCPs(..., const OGRSpatialReference* poSRS), and the removal of their older
     411equivalents using WKT strings instead of a OGRSpatialReference* instance.
     412
     413== Documentation ==
     414
     415New methods have been documented, and documentation of existing methods has
     416been changed when appropriate during the development. That said, a more
     417thorough pass will be needed. The tutorials will also have to be updated.
     418
     419== Testing ==
     420
     421The autotest suite has been adapted in a number of places since the expected
     422results have changed for a number of reasons (AXIS node exported in WKT,
     423differences in WKT and PROJ string generation). New tests have been added for the
     424new capabilities.
     425
     426It should be noted that autotest not necessarily checks everything, and issues
     427have been discovered and fixed through manual testing. The introduction of the
     428"data axis to CRS axis mapping" concept is also quite error prone, as it requires
     429setting the OAMS_TRADITIONAL_GIS_ORDER strategy in a lot of different places.
     430
     431So users and developers are kindly invited to thoroughly test GDAL once this
     432work has landed in master.
     433
     434== Implementation ==
     435
     436Done by Even Rouault, [http://www.spatialys.com Spatialys].
     437Available per [https://github.com/OSGeo/gdal/pull/1185 PR 1185]
     438Funded through [https://gdalbarn.com/ gdalbarn] sponsoring.
     439
     440While it is provided as a multiple commit for """easier""" review, it will be
     441probably squashed in a single commit for inclusion in master, as intermediate steps
     442are not all buildable, due to PROJ symbol renames having occured during the
     443development, which would break bisectability.
     444
     445== Voting history ==
     446
     447TBD