-
Notifications
You must be signed in to change notification settings - Fork 44
OGC WKT Coordinate System Issues
This document is intended to discuss some issues that arise in attempting to use OGC Well Known Text (WKT) descriptions of coordinate reference systems, especially with implementations based on older WKT specifications. It discusses various vendor implementations and issues between three sources of WKT specifications:
- The original "Simple Features" specification (ie. SF-SQL 99-049).
- The newer Coordinate Transformation Services (CT) specification (01-009) which fixes some ambiguities in the original specification (e.g. regarding units of measurement) and defines an extended form of WKT.
- The ISO 19162:2015 specification, also known as "WKT 2".
Note that many issues discussed in this page are resolved in WKT 2. However WKT 1 is still in wide use.
At this time we are aware of at least the following software packages that use some form of WKT internally, or for interchange of coordinate system descriptions:
- Oracle Spatial (WKT is used internally in MDSYS.WKT, loosely SFSQL based).
- ESRI - The Arc8 system's projection engine uses a roughly simple features compatible description for projections.
- Cadcorp - Has the ability to read and write CT 1.0 style WKT. Cadcorp wrote the CT spec.
- OGR/GDAL - Reads/writes WKT as it's internal coordinate system description format. Attempts to support old and new forms as well as forms from ESRI.
- FME - Includes WKT read/write capabilities built on OGR.
- MapGuide - Uses WKT in the SDP data access API. Roughly SF compliant.
- PostGIS - Keeps WKT in the
spatial_ref_systable, but it is up to clients to translate to PROJ.4 format for actual use. Thespatial_ref_systable may be populated using OGR generated translations. - GeoTools - can read and write either the SF or CT style of WKT.
- Apache Spatial Information System (SIS) - can read and write the SF, CT or ISO 19162 style of WKT.
The Simple Feature specification does not list a set of projections and the parameters associated with them. This leads to various selection of parameter names (and sometimes projection names) from different vendors who based their implementation on that oldest specification.
The Coordinate Transformation Services (CT) 1.0 specification lists a few projections and parameters. Please try to adhere to the projection names and parameters listed there when formatting a WKT 1 definition. For projection and parameters not listed in the CT specification, the GeoTIFF Projections List registry provides a list of WKT bindings in usage. That list also tries to relate the projections to the GeoTIFF, EPSG and PROJ.4 formulations where possible.
The ISO 19162 specification tries to fix the projection and parameter names issue by recommending to specify the EPSG (or other authority) code with all projections and parameters. This is possible because EPSG codes exist not only for Coordinate Reference Systems, but also for datum, axes, operation methods, parameters, etc.
Some vendors have different approach for Lambert Conformal Conic projection. In EPSG there is a 1SP and a 2SP form of this. ESRI merges them, and just have different parameters depending on the type.
Another issue is the formulation for Albers projection. While GDAL/OGR used longitude_of_center and latitude_of_center ESRI uses Central_meridian and latitude_of_origin. EPSG defines the "Albers Equal Area" projection with "Longitude of false origin" and "Latitude of false origin" parameters.
ESRI:
PROJECTION["Albers"],
PARAMETER["False_Easting",1000000.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",-126.0],
PARAMETER["Standard_Parallel_1",50.0],
PARAMETER["Standard_Parallel_2",58.5],
PARAMETER["Latitude_Of_Origin",45.0],
OGR:
PROJECTION["Albers"],
PARAMETER["standard_parallel_1",50],
PARAMETER["standard_parallel_2",58.5],
PARAMETER["longitude_of_center",-126],
PARAMETER["latitude_of_center",45],
PARAMETER["false_easting",1000000],
PARAMETER["false_northing",0],
In Simple Features style WKT, the name associated with a datum is the only way to identify the datum. In CT WKT the datum can also have a TOWGS84 parameter indicating it's relationship to WGS84, and an AUTHORITY parameter relating it to EPSG or some other authority space. However, in SF WKT the name itself is the only key.
By convention OGR and Cadcorp have translated the datum names in a particular way from the EPSG database in order to produce comparable names. The rule is to convert all non alphanumeric characters to underscores, then to strip any leading, trailing or repeating underscores. This produces well behaved datum names like "Nouvelle_Triangulation_Francaise".
However, other vendors have done different things. ESRI seems to follow a similar convention but prefixes all datum names with "D_" as well, giving names like "D_WGS_1972". Also they have lots of other differences for reasons that are not clear. For instance for what Cadcorp and OGR call "Nouvelle_Triangulation_Francaise", they call it "D_NTF". Oracle appears to use the raw names without cleanup. So for NTF they use "NTF (Paris meridian)".
The short result of this is that recognizing and comparing datum between different Simple Features implementations require the use of heuristic rules and a table of datum aliases.
It is worthwhile keeping in mind that the BNF grammars for WKT 1 in the SF specs, and the CT spec imply specific orders for most items. For instance the BNF for the PROJCS item in the CT spec is:
<projected cs> =
PROJCS["<name>", <geographic cs>, <projection>, {<parameter>,}* <linear unit> {,<twin axes>}{,<authority>}]
This clearly states that the PROJECTION keyword follows the GEOGCS, followed by the UNIT, AXIS and AUTHORITY items. Providing them out of order is technically a violation of the spec. On the other hand, WKT consumers are encouraged to be flexible on ordering.
The linear PARAMETER values in a PROJCS must be in terms of the linear units for that PROJCS. The main linear units are the false easting and northing type values. Thus, in common cases like a state plane zone in feet, the false easting and northing will also be in feet.
The angular PARAMETER values in a PROJCS must be in terms of the angular units of the GEOGCS. If the GEOGCS is in gradians, for instance, then all the projection angles must also be in gradians.
What units should the prime meridian appear in?
The CT 1.0 specification (7.3.14 PRIMEM) says "The units of the longitude must be inferred from the context. If the PRIMEM clause occurs inside a GEOGCS, then the longitude units will match those of the geographic coordinate system. If the PRIMEM clause occurs inside a GEOCCS, then the units will be in degrees." Note that while the GEOGCS and GEOCCS keywords look very similar, they are not the same (the forth letter is a G in one case and a C in the other case). As a rule of thumb, the prime meridian in WKT 1 shall be expressed in the same units than the latitude and longitude axes of the enclosing coordinate reference system. This rule can be applied in the geographic CRS case (GEOGCS), but can not be applied in the geocentric case (GEOCCS) since the axes are in meters; consequently the angular unit is fixed by the specification to degrees in that later case.
The SF-SQL spec (99-049) does not attempt to address the issue of units of the prime meridian. Existing ESRI EPSG translation to WKT uses degrees for prime meridian, even when the GEOGCS is in gradians as shown in their translation of EPSG 4807:
GEOGCS["GCS_NTF_Paris",
DATUM["D_NTF",
SPHEROID["Clarke_1880_IGN",6378249.2,293.46602]],
PRIMEM["Paris",2.337229166666667],
UNIT["Grad",0.015707963267948967]]
Cadcorp implements according to the CT 1.0 specification as shown in their translation of EPSG 4807:
GEOGCS["NTF (Paris)",
DATUM["Nouvelle_Triangulation_Francaise",
SPHEROID["Clarke 1880 (IGN)",6378249.2,293.466021293627,
AUTHORITY["EPSG",7011]],
TOWGS84[-168,-60,320,0,0,0,0],
AUTHORITY["EPSG",6275]],
PRIMEM["Paris",2.5969213,
AUTHORITY["EPSG",8903]],
UNIT["grad",0.015707963267949,
AUTHORITY["EPSG",9105]],
AXIS["Lat",NORTH],
AXIS["Long",EAST],
AUTHORITY["EPSG",4807]]
Oracle Spatial 8.1.7 uses the following definition for what we assume is supposed to be EPSG 4807. Interestingly it does not bother with using gradians, and it appears that the prime meridian is expressed in radians with very low precision.
GEOGCS [ "Longitude / Latitude (NTF with Paris prime meridian)",
DATUM ["NTF (Paris meridian)",
SPHEROID ["Clarke 1880 (IGN)", 6378249.200000, 293.466021]],
PRIMEM [ "", 0.000649 ],
UNIT ["Decimal Degree", 0.01745329251994330]]
ISO 19162:2015 recommends to explicitly specify the PRIMEM unit for avoiding confusion. However if the unit is omitted, then ISO 19162 confirms the CT 1.0 interpretation (prime meridian in units of the enclosing geographic CRS). In other words, the above "Cadcorp way" is the standard approach for writing WKT 1 definitions that are compatible with WKT 2 / ISO 19162.
Apache SIS follows the "Cadcorp way" by default since it in better agreement with OGC/ISO standards, but provides an option for switching to the "GDAL/ESRI way" for WKT 1 when requested. WKT 2 in always treated in the "Cadcorp way" since it is settled in the ISO 19162 standard.
In EPSG there are two methods of defining the 7 parameter bursa wolf parameters, 9606 (position vector 7-parameter) and 9607 (coordinate frame rotation). The only difference is that the sign of the rotation coefficients is reversed between them.
In GDAL/OGC the TOWGS84 values in WKT are assumed to be done using the sense in 9606 (position vector 7-parameter) and when reading a 9607 the library switches the rotation signs before putting it into a TOWGS84 chunk in WKT.
However, some WKT dump from other implementations are using the 9607 sense. For instance, this item appears to use 9607 values directly without switching the sign.
GEOGCS["DHDN",
DATUM["Deutsche_Hauptdreiecksnetz",
SPHEROID["Bessel 1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],
TOWGS84[582,105,414,-1.04,-0.35,3.08,8.3],
AUTHORITY["EPSG","6314"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["DMSH",0.0174532925199433,AUTHORITY["EPSG","9108"]],
AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4314"]]
The TOWGS84 clause in the 1.0 CT spec just talks about them being the Bursa Wolf transformation parameters (on page 22, 7.3.18). Sections 12.3.15.2 and 12.3.27 are nonspecific as to the handedness of the TOWGS84 rotations. It still not clear of whether TOWGS84 matches EPSG 9606 or EPSG 9607.
ISO 19162:2015 removes any TOWGS84 support. One problem of TOWGS84 is that it encourages the use of WGS 84 as a "pivot" datum, with every datum shift operations defined in terms of that datum. This is not always desirable. For example in the transformation from Martinique 1938 to RGAF09, any operation passing through WGS 84 will inject errors in the results. The use of TOWGS84 as a pivot system is named early-binding approach in EPSG and GIGS guidance notes, as opposed to the recommended late-binding approach. For users who really want the early-binding approach, ISO 19162 replaces TOWGS84 by a new keyword, BOUNDCRS, which resolves the sign issue by specifying explicitly which operation method to use.
Another related question is whether longitudinal projection parameters (ie. central meridian) are relative to the GEOGCS prime meridian or relative to Greenwich. While the simplest approach is to treat all longitudes as relative to Greenwich, it may be that the longitudes were intended to be relative to the prime meridian. However, a review of 7.3.11 (describing PARAMETER) in the CT 1.0 spec provides no support for this opinion, and an inspection of EPSG 25700 in Cadcorp also suggests that the central meridian is relative to Greenwich, not the prime meridian.
PROJCS["Makassar (Jakarta) / NEIEZ",
GEOGCS["Makassar (Jakarta)",
DATUM["Makassar",
SPHEROID["Bessel 1841",6377397.155,299.1528128,
AUTHORITY["EPSG","7004"]],
TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6257"]],
PRIMEM["Jakarta",106.807719444444,
AUTHORITY["EPSG","8908"]],
UNIT["DMSH",0.0174532925199433,
AUTHORITY["EPSG","9108"]],
AXIS["Lat","NORTH"],
AXIS["Long","EAST"],
AUTHORITY["EPSG","4804"]],
PROJECTION["Mercator_1SP",
AUTHORITY["EPSG","9804"]],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",110],
PARAMETER["scale_factor",0.997],
PARAMETER["false_easting",3900000],
PARAMETER["false_northing",900000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["X","EAST"],
AXIS["Y","NORTH"],
AUTHORITY["EPSG","25700"]]
Based on this, GDAL/OGR proceeds on the assumption that while parameters are in the units of the GEOGCS they are not relative the GEOGCS prime meridian.
The specification does not address the precision to which values in WKT should be stored. Some implementations, such as Oracles apparently, use rather limited precision for parameters such as Scale Factor making it difficult to compare coordinate system descriptions or even to get comparable numerical results.
The best practice is to preserve the original precision as specified in the source database, such as EPSG where possible. Given that many systems do not track precision, at least it is advisable to produce values with the equivalent of the C "%.16g" format, maintaining 16 digits of precision, capturing most of the precision of a double precision IEEE floating point value.
- ESRI seems to use Equidistant_Cylindrical for what may be Equirectangular.
Original version Written by Frank Warmerdam
Updated by Martin Desruisseaux on June 4th 2018 with the changes discussed on track 172