1--- 2title: "Migration to PROJ6/GDAL3" 3author: "Roger Bivand" 4output: 5 html_document: 6 toc: true 7 toc_float: 8 collapsed: false 9 smooth_scroll: false 10 toc_depth: 2 11bibliography: PROJ.bib 12link-citations: yes 13vignette: > 14 %\VignetteIndexEntry{Migration to PROJ6/GDAL3} 15 %\VignetteEngine{knitr::rmarkdown} 16 %\VignetteEncoding{UTF-8} 17--- 18 19# Status 22 April 2020 20 21[R spatial follows GDAL and PROJ development](https://r-spatial.org/r/2020/03/17/wkt.html) describes the main points of the status now. **sp** uses WKT2 CRS with PROJ6+/GDAL3+. **sp** and **rgdal** read, write, and project and transform objects using PROJ strings before PROJ6/GDAL3 and when **raster** calls `rgdal::rawTransform()` with PROJ6+/GDAL3+. **sp** and **rgdal** read, write, and project and transform objects using WKT2_2019 strings from PROJ6/GDAL3. The mechanism used is described in the R-spatial blog. 22 23This rolling RFC (not many comments) depicts the state of play from October 2019 to April 2020, involving a lot of reverse dependency testing of almost 1000 CRAN packages. When the development versions of **sp** (>= 1.4-2) and **rgdal** (>= 1.5-8) are submitted to CRAN, a fair number of reverse dependencies will be broken (as with recent **sf** releases). Some maintainers sensibly find that fixing the PROJ6/GDAL3 transition is sufficiently invasive to make it sensible to re-base to **sf** from **sp**. Others are regretably unresponsive so far, many find it hard to check on PROJ6+/GDAL3. 24 25For further links (in addition to the blog post), see (https://github.com/r-spatial/discuss/issues/28), and **sf** github issues: https://github.com/r-spatial/sf/issues/1146, https://github.com/r-spatial/sf/issues/1187, https://github.com/r-spatial/sf/issues/1231, 26https://github.com/r-spatial/sf/issues/1328 and many others. 27 28We've established that we should have preferred WKT over PROJ strings starting eight years ago: https://lists.osgeo.org/pipermail/gdal-dev/2012-November/034558.html (read the last paragraph), and possibly even earlier. 29 30 31# Migration to PROJ6/GDAL3 32 33## PROJ 34 35Because so much open source (and other) software uses the PROJ library and framework, many are affected when PROJ upgrades. Until very recently, PROJ has been seen as very reliable, and the changes taking place now are intended to confirm and reinforce this reliability. Before PROJ 5 (PROJ 6 is out now, PROJ 7 is coming early in 2020), the `+datum=` tag was used, perhaps with `+towgs84=` with three or seven coefficients, and possibly `+nadgrids=` where datum transformation grids were available. However, transformations from one projection to another first inversed to longitude-latitude in WGS84, then projected on to the target projection. 36 37 38 39'Fast-forward 35 years and PROJ.4 is everywhere: It provides coordinate handling for almost every geospatial program, open or closed source. Today, we see a drastical increase in the need for high accuracy GNSS coordinate handling, especially in the agricultural and construction engineering sectors. This need for geodetic-accuracy transformations is not satisfied by "classic PROJ.4". But with the ubiquity of PROJ.4, we can provide these transformations "everywhere", just by implementing them as part of PROJ.4' [@evers+knudsen17]. 40 41 42 43Following the introduction of geodetic modules and pipelines in PROJ 5 [@knudsen+evers17; @evers+knudsen17], PROJ 6 moves further. Changes in the legacy PROJ representation and WGS84 transformation hub have been coordinated through the [GDAL barn raising](https://gdalbarn.com/) initiative. Crucially WGS84 often ceases to be the pivot for moving between datums. A new OGC WKT is coming, and an SQLite EPSG file database has replaced CSV files. SRS will begin to support 3D by default, adding time too as SRS change. See also [PROJ migration notes](https://proj.org/development/migration.html). 44 45There are very useful postings on the PROJ mailing list from Martin Desruisseaux, first [proposing clarifications](https://lists.osgeo.org/pipermail/proj/2019-July/008748.html) and a [follow-up](https://lists.osgeo.org/pipermail/proj/2019-August/008750.html) including a summary: 46 47> * "Early binding" ≈ hub transformation technique. 48 49> * "Late binding" ≈ hub transformation technique NOT used, replaced by 50a more complex technique consisting in searching parameters in the 51EPSG database after the transformation context (source, target, 52epoch, area of interest) is known. 53 54> * The problem of hub transformation technique is independent of WGS84. 55It is caused by the fact that transformations to/from the hub are 56approximate. Any other hub we could invent in replacement of WGS84 57will have the same problem, unless we can invent a hub for which 58transformations are exact (I think that if such hub existed, we 59would have already heard about it). 60 61> The solution proposed by ISO 19111 (in my understanding) is: 62 63> * Forget about hub (WGS84 or other), unless the simplicity of 64early-binding is considered more important than accuracy. 65 66> * Associating a CRS to a coordinate set (geometry or raster) is no 67longer sufficient. A {CRS, epoch} tuple must be associated. ISO 6819111 calls this tuple "Coordinate metadata". From a programmatic 69API point of view, this means that getCoordinateReferenceSystem() 70method in Geometry objects (for instance) needs to be replaced by a 71getCoordinateMetadata() method. 72 73For users of North American spatial data, this [ESRI news item](https://www.esri.com/about/newsroom/arcuser/moving-from-static-spatial-reference-systems-in-2022/) gives a broad-brush picture of some of the motivations and oncoming changes. 74 75## Impacts of GDAL barnraising on **sp** workflows 76 77The following examples will contrast the behaviour of PROJ4/GDAL2 (similar to PROJ5/GDAL2, and PROJ4/GDAL1) and PROJ6+/GDAL3+. In particular, the behaviour of the `exportToProj4()` function in GDAL's `OGRSpatialReference` class has changed: 78 79> Warning 80 Use of this function is discouraged. Its behaviour in GDAL >= 3 / PROJ >= 6 is significantly different from earlier versions. In particular +datum will only encode WGS84, NAD27 and NAD83, and +towgs84/+nadgrids terms will be missing most of the time. PROJ strings to encode CRS should be considered as a a legacy solution. Using a AUTHORITY:CODE or WKT representation is the recommended way. 81 82(https://gdal.org/doxygen/classOGRSpatialReference.html#a271b3de4caf844135b0c61e634860f2b); see also (https://github.com/r-spatial/sf/issues/1187) and links therein to (https://github.com/r-spatial/discuss/issues/28) and (https://github.com/r-spatial/sf/issues/1146). 83 84This function is used for both raster and vector data read through GDAL to provide the PROJ 4 string used to specify the coordinate reference system of **sp** `"Spatial"` objects using the `"CRS"` (S4, new style) class. Such classes cannot be modified without making it impossible for users to load serialised objects, such as **sp** RDS objects from GADM [for example for Norway](https://gadm.org/download_country_v3.html). 85 86My fork of **sp** (https://github.com/rsbivand/sp) is currently at 1.3-3 or higher, and contains extra code conditionally using the development version of **rgdal** on the R-Forge repository at 1.5-1 or higher. The commented out blocks marked PROJ4/GDAL2 were generated on a Windows platform with **sp** 1.3-2 and **rgdal** 1.4-7, using PROJ4/GDAL2. The other commented out blocks were **sp** fork and **rgdal** development version, revision 886. The uncommented output is what the package build platform put there. 87 88 89```{r} 90library(sp) 91packageVersion("sp") 92``` 93 94```{r} 95## [1] '1.3.3' 96``` 97 98```{r} 99## PROJ4/GDAL2 [1] '1.3.2' 100``` 101 102The `"CRS"` class definition is unchanged going forward to maintain backward compatibility. 103 104```{r} 105getClass("CRS") 106``` 107 108```{r} 109rgdal::rgdal_extSoftVersion() 110``` 111 112```{r} 113## GDAL GDAL_with_GEOS PROJ sp 114## "3.0.2" "TRUE" "6.2.1" "1.3-3" 115``` 116 117```{r} 118## PROJ4/GDAL2 GDAL GDAL_with_GEOS PROJ.4 sp 119## PROJ4/GDAL2 "2.2.3" "TRUE" "4.9.3" "1.3-1" 120``` 121 122```{r} 123packageVersion("rgdal") 124``` 125 126```{r} 127## [1] '1.5.1' 128``` 129 130```{r} 131## PROJ4/GDAL2 [1] '1.4.7' 132``` 133 134The changes introduced affect `CRS()` when **rgdal** is available; if **rgdal** is not available, the `"CRS"` object just contains a lightly checked PROJ-style string. Because some terms are deprecated or defunct from PROJ6/GDAL3, we need to be careful. GDAL's `exportToProj4()` uses the PROJ `proj_as_proj_string()` function in its new API to return the PROJ string. Terms which are deprecated or defunct may be omitted. In the PROJ6+/GDAL3+ case, `CRS()` calls `rgdal::checkCRSArgs_ng()`, a new generation function replacing the legacy `rgdal::checkCRSArgs()` function. 135 136It passes on the input PROJ-style string to `rgdal::showSRID()`, which is a many-to-many converter. It can take PROJ-style strings, WKT strings, urn-style strings and EPSG-style strings, and converts to WKT (many types) and PROJ. The PROJ-to-PROJ conversion is equivalent to PROJ4/GDAL2 behaviour, using `importFromProj4()` and friends to instantiate an SRS object in GDAL, and `exportToProj4()` and `exportToWkt()` to emit strings. If the output string appears to be missing a specification term implied by the input, a warning is given; the warnings are at present overly cautious. 137 138```{r} 139(crs <- CRS("+proj=longlat +ellps=WGS84")) 140``` 141 142```{r} 143## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): 144## Discarded datum Unknown_based_on_WGS84_ellipsoid in CRS definition 145 146## CRS arguments: +proj=longlat +ellps=WGS84 +no_defs 147``` 148 149In `rgdal::checkCRSArgs_ng()`, `rgdal::showSRID()` may be called several times. If the passed `"CRS"` object only has a non-NA PROJ-style string, this is used to populate the SRS object, as in this case. In addition to emitting a checked PROJ string, a WKT2 string is also returned (WKT2_2018), and this string is assigned as a `comment()` to the `"CRS"` object. This representation is far more robust than the PROJ-style string, giving authorities and table look-up ID values. (*WKT comment strings are reported here split across lines, because some find right-scrolling unpretty; the real format is as a character string on one line only*). 150 151```{r} 152cat(strwrap(gsub(",", ", ", (comment(crs)))), sep="\n") 153``` 154 155```{r} 156## [1] "GEOGCRS[\"unknown\", DATUM[\"Unknown based on WGS84 ellipsoid\", 157## ELLIPSOID[\"WGS 84\", 6378137, 298.257223563, LENGTHUNIT[\"metre\", 1], 158## ID[\"EPSG\", 7030]]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\", 159## 0.0174532925199433], ID[\"EPSG\", 8901]], CS[ellipsoidal, 2], AXIS[\"longitude\", 160## east, ORDER[1], ANGLEUNIT[\"degree\", 0.0174532925199433, ID[\"EPSG\", 9122]]], 161## AXIS[\"latitude\", north, ORDER[2], ANGLEUNIT[\"degree\", 0.0174532925199433, 162## ID[\"EPSG\", 9122]]]]" 163``` 164 165For the PROJ4/GDAL2 (and similar) cases, if **rgdal** is available, `CRS()` calls `rgdal::checkCRSArgs()`, calling `RGDAL_checkCRSArgs()`, a compiled function, which calls `pj_init_plus()` in PROJ to check the validity of the string and possibly expand terms. If this succeeds, `pj_get_def()` is used to return the PROJ string. Both of these functions are part of the deprecated PROJ API that is still accessible in PROJ 6, but will soon be made defunct. 166 167```{r} 168## PROJ4/GDAL2 CRS arguments: +proj=longlat +ellps=WGS84 169``` 170 171A number of further examples will be given here, including the case of one of three `+datum=` values that are still acknowledged in GDAL3/PROJ6: `WGS84`, `NAD27` and `NAD83`. 172 173```{r} 174(crs <- CRS("+proj=longlat +datum=WGS84")) 175``` 176 177```{r} 178## CRS arguments: +proj=longlat +datum=WGS84 +no_defs 179``` 180 181```{r} 182cat(strwrap(gsub(",", ", ", (comment(crs)))), sep="\n") 183``` 184 185```{r} 186## [1] "GEOGCRS[\"unknown\", DATUM[\"World Geodetic System 1984\", 187## ELLIPSOID[\"WGS 84\", 6378137, 298.257223563, LENGTHUNIT[\"metre\", 1]], 188## ID[\"EPSG\", 6326]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\", 189## 0.0174532925199433], ID[\"EPSG\", 8901]], CS[ellipsoidal, 2], AXIS[\"longitude\", 190## east, ORDER[1], ANGLEUNIT[\"degree\", 0.0174532925199433, ID[\"EPSG\", 9122]]], 191## AXIS[\"latitude\", north, ORDER[2], ANGLEUNIT[\"degree\", 0.0174532925199433, 192## ID[\"EPSG\", 9122]]]]" 193``` 194 195```{r} 196## PROJ4/GDAL2 CRS arguments: 197## PROJ4/GDAL2 +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 198``` 199 200We know that `+datum`, `+towgs84`, `+nadgrids` and `+init` are fragile, so we'll try one: 201 202```{r} 203(crs <- CRS("+proj=longlat +towgs84=0,0,0")) 204``` 205 206```{r} 207## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum WGS_1984 in CRS definition, 208## but +towgs84= values preserved 209 210## CRS arguments: 211## +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 212``` 213 214The comment here includes the `+towgs84` parameters (three of them), while the PROJ string gives seven. 215 216```{r} 217cat(strwrap(gsub(",", ", ", (comment(crs)))), sep="\n") 218``` 219 220```{r} 221## [1] "BOUNDCRS[SOURCECRS[GEOGCRS[\"unknown\", DATUM[\"World Geodetic System 1984\", 222## ELLIPSOID[\"WGS 84\", 6378137, 298.257223563, LENGTHUNIT[\"metre\", 1]], ID[\"EPSG\", 223## 6326]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\", 0.0174532925199433], 224## ID[\"EPSG\", 8901]], CS[ellipsoidal, 2], AXIS[\"longitude\", east, ORDER[1], 225## ANGLEUNIT[\"degree\", 0.0174532925199433, ID[\"EPSG\", 9122]]], AXIS[\"latitude\", 226## north, ORDER[2], ANGLEUNIT[\"degree\", 0.0174532925199433, ID[\"EPSG\", 9122]]]]], 227## TARGETCRS[GEOGCRS[\"WGS 84\", DATUM[\"World Geodetic System 1984\", ELLIPSOID[\"WGS 228## 84\", 6378137, 298.257223563, LENGTHUNIT[\"metre\", 1]]], PRIMEM[\"Greenwich\", 0, 229## ANGLEUNIT[\"degree\", 0.0174532925199433]], CS[ellipsoidal, 2], AXIS[\"latitude\", 230## north, ORDER[1], ANGLEUNIT[\"degree\", 0.0174532925199433]], AXIS[\"longitude\", east, 231## ORDER[2], ANGLEUNIT[\"degree\", 0.0174532925199433]], ID[\"EPSG\", 4326]]], 232## ABRIDGEDTRANSFORMATION[\"Transformation from unknown to WGS84\", METHOD[\"Geocentric 233## translations (geog2D domain)\", ID[\"EPSG\", 9603]], PARAMETER[\"X-axis translation\", 234## 0, ID[\"EPSG\", 8605]], PARAMETER[\"Y-axis translation\", 0, ID[\"EPSG\", 8606]], 235## PARAMETER[\"Z-axis translation\", 0, ID[\"EPSG\", 8607]]]]" 236``` 237 238```{r} 239## PROJ4/GDAL2 CRS arguments: 240## PROJ4/GDAL2 +proj=longlat +towgs84=0,0,0 +ellps=WGS84 241``` 242 243The `+init` value is still accepted, but not repeated in the output: 244 245```{r} 246(crs <- CRS("+init=epsg:4326")) 247``` 248 249```{r} 250## CRS arguments: +proj=longlat +datum=WGS84 +no_defs 251``` 252 253```{r} 254cat(strwrap(gsub(",", ", ", (comment(crs)))), sep="\n") 255``` 256 257```{r} 258## [1] "GEOGCRS[\"WGS 84\", DATUM[\"World Geodetic System 1984\", ELLIPSOID[\"WGS 84\", 259## 6378137, 298.257223563, LENGTHUNIT[\"metre\", 1]], ID[\"EPSG\", 6326]], 260## PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\", 0.0174532925199433], ID[\"EPSG\", 261## 8901]], CS[ellipsoidal, 2], AXIS[\"longitude\", east, ORDER[1], ANGLEUNIT[\"degree\", 262## 0.0174532925199433, ID[\"EPSG\", 9122]]], AXIS[\"latitude\", north, ORDER[2], 263## ANGLEUNIT[\"degree\", 0.0174532925199433, ID[\"EPSG\", 9122]]], 264## USAGE[SCOPE[\"unknown\"], AREA[\"World\"], BBOX[-90, -180, 90, 180]]]" 265``` 266 267```{r} 268## PROJ4/GDAL2 CRS arguments: 269## PROJ4/GDAL2 +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 270## PROJ4/GDAL2 +towgs84=0,0,0 271``` 272 273An early warning of difficulties with discarded `+datum` values came with the GB datum: 274 275```{r, eval=FALSE} 276(crs <- CRS("+init=epsg:27700")) 277``` 278 279```{r} 280## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): 281## Discarded datum OSGB_1936 in CRS definition 282 283## CRS arguments: 284## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 285## +y_0=-100000 +ellps=airy +units=m +no_defs 286``` 287 288Here, the comment has all the information that would be needed to carry out coordinate operations, but the PROJ string is defective for GDAL3/PROJ6, giving at best ballpark accuracy. 289 290```{r, eval=FALSE} 291cat(strwrap(gsub(",", ", ", (comment(crs)))), sep="\n") 292``` 293 294```{r} 295## [1] "PROJCRS[\"OSGB 1936 / British National Grid\", BASEGEOGCRS[\"OSGB 1936\", 296## DATUM[\"OSGB 1936\", ELLIPSOID[\"Airy 1830\", 6377563.396, 299.3249646, 297## LENGTHUNIT[\"metre\", 1]]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\", 298## 0.0174532925199433]], ID[\"EPSG\", 4277]], CONVERSION[\"British National Grid\", 299## METHOD[\"Transverse Mercator\", ID[\"EPSG\", 9807]], PARAMETER[\"Latitude of natural 300## origin\", 49, ANGLEUNIT[\"degree\", 0.0174532925199433], ID[\"EPSG\", 8801]], 301## PARAMETER[\"Longitude of natural origin\", -2, ANGLEUNIT[\"degree\", 302## 0.0174532925199433], ID[\"EPSG\", 8802]], PARAMETER[\"Scale factor at natural 303## origin\", 0.9996012717, SCALEUNIT[\"unity\", 1], ID[\"EPSG\", 8805]], 304## PARAMETER[\"False easting\", 400000, LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8806]], 305## PARAMETER[\"False northing\", -100000, LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8807]], 306## ID[\"EPSG\", 19916]], CS[Cartesian, 2], AXIS[\"(E)\", east, ORDER[1], 307## LENGTHUNIT[\"metre\", 1, ID[\"EPSG\", 9001]]], AXIS[\"(N)\", north, ORDER[2], 308## LENGTHUNIT[\"metre\", 1, ID[\"EPSG\", 9001]]], USAGE[SCOPE[\"unknown\"], AREA[\"UK - 309## Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"], BBOX[49.75, -9.2, 61.14, 2.88]]]" 310``` 311 312In PROJ4/GDAL2, the `+datum` and `+towgs84` values give much better transformation accuracy from the PROJ string. 313 314```{r} 315## PROJ4/GDAL2 CRS arguments: 316## PROJ4/GDAL2 +init=epsg:27700 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 317## PROJ4/GDAL2 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs 318## PROJ4/GDAL2 +ellps=airy 319## PROJ4/GDAL2 +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894 320``` 321 322From **sp** 1.3-3 (my fork), `CRS()` takes a third argument to the legacy two, `projargs=` and `doCheckCRSArgs=TRUE`: `SRS_string=NULL`. This can take any string that `rgdal::showSRID()` can handle. The warning follows use of `exportToProj4()`: 323 324```{r} 325run <- rgdal::new_proj_and_gdal() 326``` 327 328 329```{r, eval=run} 330(crs <- CRS(SRS_string = "EPSG:27700")) 331``` 332 333```{r} 334## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"): 335## Discarded datum OSGB_1936 in CRS definition 336 337## CRS arguments: 338## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 339## +y_0=-100000 +ellps=airy +units=m +no_defs 340``` 341 342```{r, eval=run} 343cat(strwrap(gsub(",", ", ", (comment(crs)))), sep="\n") 344``` 345 346```{r} 347## [1] "PROJCRS[\"OSGB 1936 / British National Grid\", BASEGEOGCRS[\"OSGB 1936\", 348## DATUM[\"OSGB 1936\", ELLIPSOID[\"Airy 1830\", 6377563.396, 299.3249646, 349## LENGTHUNIT[\"metre\", 1]]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\", 350## 0.0174532925199433]], ID[\"EPSG\", 4277]], CONVERSION[\"British National Grid\", 351## METHOD[\"Transverse Mercator\", ID[\"EPSG\", 9807]], PARAMETER[\"Latitude of natural 352## origin\", 49, ANGLEUNIT[\"degree\", 0.0174532925199433], ID[\"EPSG\", 8801]], 353## PARAMETER[\"Longitude of natural origin\", -2, ANGLEUNIT[\"degree\", 354## 0.0174532925199433], ID[\"EPSG\", 8802]], PARAMETER[\"Scale factor at natural 355## origin\", 0.9996012717, SCALEUNIT[\"unity\", 1], ID[\"EPSG\", 8805]], 356## PARAMETER[\"False easting\", 400000, LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8806]], 357## PARAMETER[\"False northing\", -100000, LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8807]]], 358## CS[Cartesian, 2], AXIS[\"(E)\", east, ORDER[1], LENGTHUNIT[\"metre\", 1]], 359## AXIS[\"(N)\", north, ORDER[2], LENGTHUNIT[\"metre\", 1]], USAGE[SCOPE[\"unknown\"], 360## AREA[\"UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"], BBOX[49.75, 361## -9.2, 61.14, 2.88]], ID[\"EPSG\", 27700]]" 362``` 363 364We can also use a more detailed PROJ string, but without any improvement of the output PROJ representation; the comment is OK: 365 366```{r, eval=run} 367(crs <- CRS("+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs")) 368``` 369 370```{r} 371## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): 372## Discarded datum OSGB_1936 in CRS definition 373 374## CRS arguments: 375## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 376## +y_0=-100000 +ellps=airy +units=m +no_defs 377``` 378 379 380```{r, eval=run} 381cat(strwrap(gsub(",", ", ", (comment(crs)))), sep="\n") 382``` 383 384```{r} 385## [1] "PROJCRS[\"unknown\", BASEGEOGCRS[\"unknown\", DATUM[\"OSGB 1936\", 386## ELLIPSOID[\"Airy 1830\", 6377563.396, 299.3249646, LENGTHUNIT[\"metre\", 1]], 387## ID[\"EPSG\", 6277]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\", 388## 0.0174532925199433], ID[\"EPSG\", 8901]]], CONVERSION[\"unknown\", METHOD[\"Transverse 389## Mercator\", ID[\"EPSG\", 9807]], PARAMETER[\"Latitude of natural origin\", 49, 390## ANGLEUNIT[\"degree\", 0.0174532925199433], ID[\"EPSG\", 8801]], PARAMETER[\"Longitude 391## of natural origin\", -2, ANGLEUNIT[\"degree\", 0.0174532925199433], ID[\"EPSG\", 392## 8802]], PARAMETER[\"Scale factor at natural origin\", 0.9996012717, 393## SCALEUNIT[\"unity\", 1], ID[\"EPSG\", 8805]], PARAMETER[\"False easting\", 400000, 394## LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8806]], PARAMETER[\"False northing\", -100000, 395## LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8807]]], CS[Cartesian, 2], AXIS[\"(E)\", east, 396## ORDER[1], LENGTHUNIT[\"metre\", 1, ID[\"EPSG\", 9001]]], AXIS[\"(N)\", north, 397## ORDER[2], LENGTHUNIT[\"metre\", 1, ID[\"EPSG\", 9001]]]]" 398``` 399 400`rgdal::showSRID()` is quite versatile, so we can display a multiline WKT string as well: 401 402```{r, eval=run} 403if (packageVersion("rgdal") >= "1.5.1") cat(rgdal::showSRID("+init=epsg:27700", format="WKT2_2018", multiline="YES", prefer_proj=FALSE), "\n") 404``` 405 406```{r} 407## PROJCRS["OSGB 1936 / British National Grid", 408## BASEGEOGCRS["OSGB 1936", 409## DATUM["OSGB 1936", 410## ELLIPSOID["Airy 1830",6377563.396,299.3249646, 411## LENGTHUNIT["metre",1]]], 412## PRIMEM["Greenwich",0, 413## ANGLEUNIT["degree",0.0174532925199433]], 414## ID["EPSG",4277]], 415## CONVERSION["British National Grid", 416## METHOD["Transverse Mercator", 417## ID["EPSG",9807]], 418## PARAMETER["Latitude of natural origin",49, 419## ANGLEUNIT["degree",0.0174532925199433], 420## ID["EPSG",8801]], 421## PARAMETER["Longitude of natural origin",-2, 422## ANGLEUNIT["degree",0.0174532925199433], 423## ID["EPSG",8802]], 424## PARAMETER["Scale factor at natural origin",0.9996012717, 425## SCALEUNIT["unity",1], 426## ID["EPSG",8805]], 427## PARAMETER["False easting",400000, 428## LENGTHUNIT["metre",1], 429## ID["EPSG",8806]], 430## PARAMETER["False northing",-100000, 431## LENGTHUNIT["metre",1], 432## ID["EPSG",8807]], 433## ID["EPSG",19916]], 434## CS[Cartesian,2], 435## AXIS["(E)",east, 436## ORDER[1], 437## LENGTHUNIT["metre",1, 438## ID["EPSG",9001]]], 439## AXIS["(N)",north, 440## ORDER[2], 441## LENGTHUNIT["metre",1, 442## ID["EPSG",9001]]], 443## USAGE[ 444## SCOPE["unknown"], 445## AREA["UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E"], 446## BBOX[49.75,-9.2,61.14,2.88]]] 447``` 448 449So far, `"CRS"` objects created on creation from **sp**, and from reading rasters and vectors in **rgdal**, are furnished with WKT comments. These are used to instantiate SRS objects when writing raster and vector objects. What remains is to convert the compiled `transform()` function and `spTransform()` methods to use the WKT comments if available instead of the PROJ string in the `"CRS"` object in each `"Spatial"` object. 450 451## Coordinate operations 452 453```{r} 454library(rgdal) 455``` 456 457 458```{r} 459## rgdal: version: 1.5-1, (SVN revision 870) 460## Geospatial Data Abstraction Library extensions to R successfully loaded 461## Loaded GDAL runtime: GDAL 3.0.2, released 2019/10/28 462## Path to GDAL shared files: 463## GDAL binary built with GEOS: TRUE 464## Loaded PROJ.4 runtime: Rel. 6.2.1, November 1st, 2019, [PJ_VERSION: 621] 465## Path to PROJ.4 shared files: (autodetected) 466## Linking to sp version: 1.3-3 467``` 468 469 470```{r} 471run <- new_proj_and_gdal() 472``` 473 474We can first show what happens when searching for instantiable coordinate operations. These are coordinate operations that can be created by look-up in the PROJ database given the information in the input description. Here the datum has been discarded. 475 476```{r, eval=run} 477(discarded_datum <- showSRID("EPSG:27700", "PROJ")) 478``` 479 480```{r} 481## Warning in showSRID("EPSG:27700", "PROJ"): Discarded datum OSGB_1936 in CRS 482## definition 483``` 484 485Consequently, when searching for coordinate operations to transform to geographical coordinates and the WGS84 datum, and using `importFromProj4()` to instantiate the source SRS, the only operation found only has the Airy ellipse information to guide its search. The `list_coordOps()` function takes two SRS descriptions, and returns coordinate operations found by look-up. We need to add the `type` tag to a PROJ string here. 486 487```{r, eval=run} 488(x <- list_coordOps(paste0(discarded_datum, " +type=crs"), "EPSG:4326")) 489``` 490 491```{r} 492## Candidate coordinate operations found: 1 493## Strict containment: FALSE 494## Visualization order: TRUE 495## Source: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 496## +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs 497## Target: EPSG:4326 498## Best instantiable operation has only ballpark accuracy 499## Description: Inverse of unknown + Ballpark geographic offset from unknown to 500## WGS 84 + axis order change (2D) 501## Definition: +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 502## +k=0.9996012717 +x_0=400000 +y_0=-100000 503## +ellps=airy +step +proj=unitconvert +xy_in=rad 504## +xy_out=deg 505``` 506 507The `best_instantiable_coordOp()` function retuns the best instantiable coordinate operation, in this case the only one, with a description attribute. This is the current best available in calling `spTransform()` with `"CRS"` objects with discarded `+datum=` tags. 508 509```{r, eval=run} 510best_instantiable_coordOp(x) 511``` 512 513```{r} 514## Warning in best_instantiable_coordOp(x): Best instantiable operation has only 515## ballpark accuracy 516## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 517## +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=unitconvert +xy_in=rad +xy_out=deg" 518## attr(,"description") 519## [1] "Inverse of unknown + Ballpark geographic offset from unknown to WGS 84 + 520## axis order change (2D)" 521``` 522 523If we add back the discarded `+datum=` tag with a valid value, we give the search process what it needs to find more coordinate operations. 524 525```{r, eval=run} 526list_coordOps(paste0(discarded_datum, " +datum=OSGB36 +type=crs"), "EPSG:4326") 527``` 528 529Now we have 7 operations, one the ballpark accuracy operation with unknown accuracy found above, 5 others that can be instantiated, of which the best has 2m accuracy, and an operation that cannot be instantiated because a publically available grid is missing. On download and installation of this grid, 1m accuracy could be achieved. 530 531```{r} 532## Candidate coordinate operations found: 7 533## Strict containment: FALSE 534## Visualization order: TRUE 535## Source: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 536## +y_0=-100000 +ellps=airy +units=m +no_defs 537## +datum=OSGB36 +type=crs 538## Target: EPSG:4326 539## Best instantiable operation has accuracy: 2 m 540## Description: Inverse of unknown + axis order change (2D) + OSGB 1936 to WGS 541## 84 (6) + axis order change (2D) 542## Definition: +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 543## +k=0.9996012717 +x_0=400000 +y_0=-100000 544## +ellps=airy +step +proj=push +v_3 +step +proj=cart 545## +ellps=airy +step +proj=helmert +x=446.448 546## +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842 547## +s=-20.489 +convention=position_vector +step +inv 548## +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step 549## +proj=unitconvert +xy_in=rad +xy_out=deg 550## Operation 6 is lacking 1 grid with accuracy 1 m 551## Missing grid: OSTN15_NTv2_OSGBtoETRS.gsb 552## URL: https://download.osgeo.org/proj/proj-datumgrid-europe-latest.zip 553``` 554 555If we pass the source WKT string to `list_coordOps()` we get 6 operations back, now without the ballpark accuracy operation; this would be the situation in which WKT comments if present were used to construct the source and target SRS: 556 557```{r, eval=run} 558wkt_source_datum <- showSRID("EPSG:27700", "WKT2") 559wkt_target_datum <- showSRID("EPSG:4326", "WKT2") 560(x <- list_coordOps(wkt_source_datum, wkt_target_datum)) 561``` 562 563```{r} 564## Candidate coordinate operations found: 6 565## Strict containment: FALSE 566## Visualization order: TRUE 567## Source: PROJCRS["OSGB 1936 / British National Grid", BASEGEOGCRS["OSGB 568## 1936", DATUM["OSGB 1936", ELLIPSOID["Airy 1830", 569## 6377563.396, 299.3249646, LENGTHUNIT["metre", 1]]], 570## PRIMEM["Greenwich", 0, ANGLEUNIT["degree", 571## 0.0174532925199433]], ID["EPSG", 4277]], 572## CONVERSION["British National Grid", METHOD["Transverse 573## Mercator", ID["EPSG", 9807]], PARAMETER["Latitude of 574## natural origin", 49, ANGLEUNIT["degree", 575## 0.0174532925199433], ID["EPSG", 8801]], 576## PARAMETER["Longitude of natural origin", -2, 577## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG", 578## 8802]], PARAMETER["Scale factor at natural origin", 579## 0.9996012717, SCALEUNIT["unity", 1], ID["EPSG", 8805]], 580## PARAMETER["False easting", 400000, LENGTHUNIT["metre", 581## 1], ID["EPSG", 8806]], PARAMETER["False northing", 582## -100000, LENGTHUNIT["metre", 1], ID["EPSG", 8807]]], 583## CS[Cartesian, 2], AXIS["(E)", east, ORDER[1], 584## LENGTHUNIT["metre", 1]], AXIS["(N)", north, ORDER[2], 585## LENGTHUNIT["metre", 1]], USAGE[SCOPE["unknown"], 586## AREA["UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W 587## to 3°33'E"], BBOX[49.75, -9.2, 61.14, 2.88]], 588## ID["EPSG", 27700]] 589## Target: GEOGCRS["WGS 84", DATUM["World Geodetic System 1984", 590## ELLIPSOID["WGS 84", 6378137, 298.257223563, 591## LENGTHUNIT["metre", 1]]], PRIMEM["Greenwich", 0, 592## ANGLEUNIT["degree", 0.0174532925199433]], 593## CS[ellipsoidal, 2], AXIS["geodetic latitude (Lat)", 594## north, ORDER[1], ANGLEUNIT["degree", 595## 0.0174532925199433]], AXIS["geodetic longitude (Lon)", 596## east, ORDER[2], ANGLEUNIT["degree", 597## 0.0174532925199433]], USAGE[SCOPE["unknown"], 598## AREA["World"], BBOX[-90, -180, 90, 180]], ID["EPSG", 599## 4326]] 600## Best instantiable operation has accuracy: 2 m 601## Description: Inverse of British National Grid + OSGB 1936 to WGS 84 (6) + 602## axis order change (2D) 603## Definition: +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 604## +k=0.9996012717 +x_0=400000 +y_0=-100000 605## +ellps=airy +step +proj=push +v_3 +step +proj=cart 606## +ellps=airy +step +proj=helmert +x=446.448 607## +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842 608## +s=-20.489 +convention=position_vector +step +inv 609## +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step 610## +proj=unitconvert +xy_in=rad +xy_out=deg 611## Operation 6 is lacking 1 grid with accuracy 1 m 612## Missing grid: OSTN15_NTv2_OSGBtoETRS.gsb 613## URL: https://download.osgeo.org/proj/proj-datumgrid-europe-latest.zip 614``` 615 616```{r, eval=run} 617best_instantiable_coordOp(x) 618``` 619 620```{r} 621## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 622## +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=push +v_3 +step +proj=cart 623## +ellps=airy +step +proj=helmert +x=446.448 +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 624## +rz=0.842 +s=-20.489 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 625## +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg" 626## attr(,"description") 627## [1] "Inverse of British National Grid + OSGB 1936 to WGS 84 (6) + axis order change (2D)" 628``` 629 630Turning to a different example, a Brazilian UTM25S Corrego Alegre 1970-72 datum, and looking for coordinate operations transform to UTM25S SIRGAS 2000, we see that a `+towgs84=` tag was preserved. 631 632```{r, eval=run} 633discarded_datum <- showSRID("EPSG:22525", "PROJ") 634``` 635 636```{r} 637## Warning in showSRID("EPSG:22525", "PROJ"): Discarded datum Corrego_Alegre_1970-72 in CRS definition, 638## but +towgs84= values preserved 639``` 640 641This meant that while only ballpark accuracy is reported, a Helmert transformation is carried out: 642 643```{r, eval=run} 644(x <- list_coordOps(paste0(discarded_datum, " +type=crs"), "EPSG:31985")) 645``` 646 647 648```{r} 649## Candidate coordinate operations found: 1 650## Strict containment: FALSE 651## Visualization order: TRUE 652## Source: +proj=utm +zone=25 +south +ellps=intl 653## +towgs84=-205.57,168.77,-4.12,0,0,0,0 +units=m +no_defs 654## +type=crs 655## Target: EPSG:31985 656## Best instantiable operation has only ballpark accuracy 657## Description: Inverse of UTM zone 25S + Transformation from unknown to WGS84 658## + Inverse of SIRGAS 2000 to WGS 84 (1) + UTM zone 659## 25S 660## Definition: +proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl 661## +step +proj=push +v_3 +step +proj=cart +ellps=intl 662## +step +proj=helmert +x=-205.57 +y=168.77 +z=-4.12 663## +rx=0 +ry=0 +rz=0 +s=0 +convention=position_vector 664## +step +inv +proj=cart +ellps=GRS80 +step +proj=pop 665## +v_3 +step +proj=utm +zone=25 +south +ellps=GRS80 666``` 667 668 669```{r, eval=run} 670best_instantiable_coordOp(x) 671``` 672 673```{r} 674## Warning in best_instantiable_coordOp(x): Best instantiable operation has only 675## ballpark accuracy 676## [1] "+proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl +step 677## +proj=push +v_3 +step +proj=cart +ellps=intl +step +proj=helmert +x=-205.57 678## +y=168.77 +z=-4.12 +rx=0 +ry=0 +rz=0 +s=0 +convention=position_vector +step +inv 679## +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=utm +zone=25 +south 680## +ellps=GRS80" 681## attr(,"description") 682## [1] "Inverse of UTM zone 25S + Transformation from unknown to WGS84 + Inverse of 683## SIRGAS 2000 to WGS 84 (1) + UTM zone 25S" 684``` 685 686If we move to input WKT strings when searching for coordinate operations, we get to the same result with unknown accuracy but using a Helmert transformation: 687 688```{r, eval=run} 689wkt_source_datum <- showSRID("EPSG:22525", "WKT2") 690wkt_target_datum <- showSRID("EPSG:31985", "WKT2") 691(x <- list_coordOps(wkt_source_datum, wkt_target_datum)) 692``` 693 694```{r} 695## Candidate coordinate operations found: 1 696## Strict containment: FALSE 697## Visualization order: TRUE 698## Source: BOUNDCRS[SOURCECRS[PROJCRS["Corrego Alegre 1970-72 / UTM zone 699## 25S", BASEGEOGCRS["Corrego Alegre 1970-72", 700## DATUM["Corrego Alegre 1970-72", 701## ELLIPSOID["International 1924", 6378388, 297, 702## LENGTHUNIT["metre", 1]]], PRIMEM["Greenwich", 0, 703## ANGLEUNIT["degree", 0.0174532925199433]], ID["EPSG", 704## 4225]], CONVERSION["UTM zone 25S", METHOD["Transverse 705## Mercator", ID["EPSG", 9807]], PARAMETER["Latitude of 706## natural origin", 0, ANGLEUNIT["degree", 707## 0.0174532925199433], ID["EPSG", 8801]], 708## PARAMETER["Longitude of natural origin", -33, 709## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG", 710## 8802]], PARAMETER["Scale factor at natural origin", 711## 0.9996, SCALEUNIT["unity", 1], ID["EPSG", 8805]], 712## PARAMETER["False easting", 500000, LENGTHUNIT["metre", 713## 1], ID["EPSG", 8806]], PARAMETER["False northing", 714## 10000000, LENGTHUNIT["metre", 1], ID["EPSG", 8807]]], 715## CS[Cartesian, 2], AXIS["(E)", east, ORDER[1], 716## LENGTHUNIT["metre", 1]], AXIS["(N)", north, ORDER[2], 717## LENGTHUNIT["metre", 1]], USAGE[SCOPE["unknown"], 718## AREA["Brazil - east of 36°W onshore"], BBOX[-10.1, -36, 719## -4.99, -34.74]], ID["EPSG", 22525]]], 720## TARGETCRS[GEOGCRS["WGS 84", DATUM["World Geodetic 721## System 1984", ELLIPSOID["WGS 84", 6378137, 722## 298.257223563, LENGTHUNIT["metre", 1]]], 723## PRIMEM["Greenwich", 0, ANGLEUNIT["degree", 724## 0.0174532925199433]], CS[ellipsoidal, 2], 725## AXIS["latitude", north, ORDER[1], ANGLEUNIT["degree", 726## 0.0174532925199433]], AXIS["longitude", east, ORDER[2], 727## ANGLEUNIT["degree", 0.0174532925199433]], ID["EPSG", 728## 4326]]], ABRIDGEDTRANSFORMATION["Corrego Alegre 1970-72 729## to WGS 84 (3)", VERSION["PBS-Bra 1983"], 730## METHOD["Geocentric translations (geog2D domain)", 731## ID["EPSG", 9603]], PARAMETER["X-axis translation", 732## -205.57, ID["EPSG", 8605]], PARAMETER["Y-axis 733## translation", 168.77, ID["EPSG", 8606]], 734## PARAMETER["Z-axis translation", -4.12, ID["EPSG", 735## 8607]], USAGE[SCOPE["Medium and small scale mapping."], 736## AREA["Brazil - Corrego Alegre 1970-1972"], BBOX[-33.78, 737## -58.16, -2.68, -34.74]], ID["EPSG", 6192], 738## REMARK["Formed by concatenation of tfms codes 6191 and 739## 1877. Used by Petrobras and ANP until February 2005 740## when replaced by Corrego Alegre 1970-72 to WGS 84 (4) 741## (tfm code 6194)."]]] 742## Target: 743## BOUNDCRS[SOURCECRS[PROJCRS["SIRGAS 2000 / UTM zone 25S", 744## BASEGEOGCRS["SIRGAS 2000", DATUM["Sistema de Referencia 745## Geocentrico para las AmericaS 2000", ELLIPSOID["GRS 746## 1980", 6378137, 298.257222101, LENGTHUNIT["metre", 747## 1]]], PRIMEM["Greenwich", 0, ANGLEUNIT["degree", 748## 0.0174532925199433]], ID["EPSG", 4674]], 749## CONVERSION["UTM zone 25S", METHOD["Transverse 750## Mercator", ID["EPSG", 9807]], PARAMETER["Latitude of 751## natural origin", 0, ANGLEUNIT["degree", 752## 0.0174532925199433], ID["EPSG", 8801]], 753## PARAMETER["Longitude of natural origin", -33, 754## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG", 755## 8802]], PARAMETER["Scale factor at natural origin", 756## 0.9996, SCALEUNIT["unity", 1], ID["EPSG", 8805]], 757## PARAMETER["False easting", 500000, LENGTHUNIT["metre", 758## 1], ID["EPSG", 8806]], PARAMETER["False northing", 759## 10000000, LENGTHUNIT["metre", 1], ID["EPSG", 8807]]], 760## CS[Cartesian, 2], AXIS["(E)", east, ORDER[1], 761## LENGTHUNIT["metre", 1]], AXIS["(N)", north, ORDER[2], 762## LENGTHUNIT["metre", 1]], USAGE[SCOPE["unknown"], 763## AREA["Brazil - 36°W to 30°W"], BBOX[-23.8, -36, 4.19, 764## -29.99]], ID["EPSG", 31985]]], TARGETCRS[GEOGCRS["WGS 765## 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 766## 84", 6378137, 298.257223563, LENGTHUNIT["metre", 1]]], 767## PRIMEM["Greenwich", 0, ANGLEUNIT["degree", 768## 0.0174532925199433]], CS[ellipsoidal, 2], 769## AXIS["latitude", north, ORDER[1], ANGLEUNIT["degree", 770## 0.0174532925199433]], AXIS["longitude", east, ORDER[2], 771## ANGLEUNIT["degree", 0.0174532925199433]], ID["EPSG", 772## 4326]]], ABRIDGEDTRANSFORMATION["SIRGAS 2000 to WGS 84 773## (1)", VERSION["OGP-C&S America"], METHOD["Geocentric 774## translations (geog2D domain)", ID["EPSG", 9603]], 775## PARAMETER["X-axis translation", 0, ID["EPSG", 8605]], 776## PARAMETER["Y-axis translation", 0, ID["EPSG", 8606]], 777## PARAMETER["Z-axis translation", 0, ID["EPSG", 8607]], 778## USAGE[SCOPE["Accuracy 1m."], AREA["Latin America - 779## SIRGAS 2000 by country"], BBOX[-59.87, -122.19, 32.72, 780## -25.28]], ID["EPSG", 15894]]] 781## Best instantiable operation has only ballpark accuracy 782## Description: Inverse of UTM zone 25S + Corrego Alegre 1970-72 to WGS 84 (3) 783## + Inverse of SIRGAS 2000 to WGS 84 (1) + UTM zone 784## 25S 785## Definition: +proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl 786## +step +proj=push +v_3 +step +proj=cart +ellps=intl 787## +step +proj=helmert +x=-205.57 +y=168.77 +z=-4.12 788## +step +inv +proj=cart +ellps=GRS80 +step +proj=pop 789## +v_3 +step +proj=utm +zone=25 +south +ellps=GRS80 790``` 791 792If we define the look-up terms directly, in this case and unlike that for the OSGB36 datum, we find more operations than when using WKT strings. Now we find one operation with known accuracy (5m), and that a further operation would be possible if a required grid had been present, this grid is not in the PROJ download archive, as no URL is given. 793 794```{r, eval=run} 795(x <- list_coordOps("EPSG:22525", "EPSG:31985")) 796``` 797 798```{r} 799## Candidate coordinate operations found: 2 800## Strict containment: FALSE 801## Visualization order: TRUE 802## Source: EPSG:22525 803## Target: EPSG:31985 804## Best instantiable operation has accuracy: 5 m 805## Description: Inverse of UTM zone 25S + Corrego Alegre 1970-72 to SIRGAS 2000 806## (2) + UTM zone 25S 807## Definition: +proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl 808## +step +proj=push +v_3 +step +proj=cart +ellps=intl 809## +step +proj=helmert +x=-206.05 +y=168.28 +z=-3.82 810## +step +inv +proj=cart +ellps=GRS80 +step +proj=pop 811## +v_3 +step +proj=utm +zone=25 +south +ellps=GRS80 812## Operation 2 is lacking 1 grid with accuracy 2 m 813## Missing grid: CA7072_003.gsb 814``` 815 816 817 818```{r, eval=run} 819best_instantiable_coordOp(x) 820``` 821 822 823```{r} 824## [1] "+proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl +step 825## +proj=push +v_3 +step +proj=cart +ellps=intl +step +proj=helmert +x=-206.05 826## +y=168.28 +z=-3.82 +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step 827## +proj=utm +zone=25 +south +ellps=GRS80" 828## attr(,"description") 829## [1] "Inverse of UTM zone 25S + Corrego Alegre 1970-72 to SIRGAS 2000 (2) + 830## UTM zone 25S" 831``` 832 833## Transformations 834 835The problem raised by the modernisation of PROJ can be encapsulated in this example of Scotland (once again the commented out runs are from a PROJ62/GDAL30 platform, the live runs from the platform building the vignette): 836 837```{r} 838scot_BNG <- readOGR(system.file("vectors", package="rgdal"), "scot_BNG") 839``` 840 841```{r} 842## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS = 843## dumpSRS): Discarded datum OSGB_1936 in CRS definition: +proj=tmerc +lat_0=49 844## +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs 845## OGR data source with driver: ESRI Shapefile 846## Source: "/home/rsb/lib/r_libs/rgdal/vectors", layer: "scot_BNG" 847## with 56 features 848## It has 13 fields 849## Warning in OGRSpatialRef(dsn = dsn, layer = layer, morphFromESRI = 850## morphFromESRI, : Discarded datum OSGB_1936 in CRS definition: +proj=tmerc 851## +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy 852## +units=m +no_defs 853## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded 854## datum Unknown_based_on_Airy_1830_ellipsoid in CRS definition 855``` 856 857On reading with PROJ6/GDAL3, we discard the `+datum` tag, so losing information and ending up with a bounding box with positional accuracy degraded to an unknown extent, using the pre-PROJ6 adaptation CRS representation: 858 859```{r} 860(load_status <- get_transform_wkt_comment()) 861``` 862 863```{r} 864## [1] TRUE 865``` 866 867(Again, commented output blocks are run on a PROJ6/GDAL3 platform). This demonstration was impacted by PROJ 6.3.0 fragility, with the `+ellps=` tag vanishing (with the `+units=` tag) apparently arbitrarily. A tentative diagnosis is that the datum node of the OGRSpatialReference object becomes unavailable, although why some change in PROJ leads to this aberation is unknown. 868 869```{r} 870set_transform_wkt_comment(FALSE) 871scot_LL <- spTransform(scot_BNG, CRS("+proj=longlat +datum=WGS84")) 872(b0 <- bbox(scot_LL)) 873``` 874 875```{r} 876## min max 877## x -8.621387 -0.753056 878## y 54.626555 60.843843 879``` 880 881We might fake the PROJ string by extracting the `"CRS"` object: 882 883```{r} 884(crs <- slot(scot_BNG, "proj4string")) 885``` 886 887```{r} 888## CRS arguments: 889## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 890## +y_0=-100000 +ellps=airy +units=m +no_defs 891``` 892 893 894```{r} 895slot(crs, "projargs") 896``` 897 898```{r} 899## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 900## +ellps=airy +units=m +no_defs" 901``` 902 903and updating the PROJ string in-place if we know a sensible incantation: 904 905```{r} 906slot(crs, "projargs") <- paste0(slot(crs, "projargs"), " +datum=OSGB36") 907slot(scot_BNG, "proj4string") <- crs 908scot_LL1 <- spTransform(scot_BNG, CRS("+proj=longlat +datum=WGS84")) 909(b1 <- bbox(scot_LL1)) 910``` 911 912```{r} 913## min max 914## x -8.622158 -0.755071 915## y 54.626633 60.843232 916``` 917 918The two bounding boxes differ a good deal: 919 920```{r} 921all.equal(b0, b1, scale=1) 922``` 923 924```{r} 925## [1] "Mean absolute difference: 0.0008689328" 926``` 927 928The SW corner of Scotland differs by about 50m, the NE corner by almost 130m: 929 930```{r} 931diag(spDists(t(b0), t(b1), longlat=TRUE))*1000 932``` 933 934```{r} 935## [1] 50.56708 128.99003 936``` 937 938 939Re-instating the use of WKT comments lets us use the transformation path which instantiates the coordinate operation from the WKT strings in the `"CRS"` object comments: 940 941```{r} 942set_transform_wkt_comment(load_status) 943``` 944 945```{r} 946scot_BNG <- readOGR(system.file("vectors", package="rgdal"), "scot_BNG") 947``` 948 949```{r} 950## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS = 951## dumpSRS): Discarded datum OSGB_1936 in CRS definition: +proj=tmerc +lat_0=49 952## +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs 953## OGR data source with driver: ESRI Shapefile 954## Source: "/home/rsb/lib/r_libs/rgdal/vectors", layer: "scot_BNG" 955## with 56 features 956## It has 13 fields 957## Warning in OGRSpatialRef(dsn = dsn, layer = layer, morphFromESRI = 958## morphFromESRI, : Discarded datum OSGB_1936 in CRS definition: +proj=tmerc 959## +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy 960## +units=m +no_defs 961## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded 962## datum Unknown_based_on_Airy_1830_ellipsoid in CRS definition 963``` 964 965In **rgdal** 1.5-2, the coordinate operation lookup is done once for the first matrix of coordinates, and passed then on to subsequent transformations: 966 967```{r} 968system.time(scot_LL2 <- spTransform(scot_BNG, CRS("+proj=longlat +datum=WGS84"))) 969``` 970 971```{r} 972## user system elapsed 973## 0.095 0.000 0.096 974``` 975 976 977Timings for **rgdal** 1.5-1, when look-up was repeated for each matrix of coordinates: 978 979```{r} 980## user system elapsed 981## 3.431 0.037 3.546 982``` 983 984```{r} 985(b2 <- bbox(scot_LL2)) 986``` 987 988```{r} 989## min max 990## x -8.622158 -0.7550709 991## y 54.626633 60.8432318 992``` 993 994```{r} 995all.equal(b1, b2, scale=1) 996``` 997 998```{r} 999## [1] "Mean absolute difference: 3.361822e-08" 1000``` 1001 1002The coordinate operation last used may be retrieved with (default empty): 1003 1004```{r} 1005get_last_coordOp() 1006``` 1007 1008```{r} 1009## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 1010## +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=push +v_3 +step +proj=cart 1011## +ellps=airy +step +proj=helmert +x=446.448 +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 1012## +rz=0.842 +s=-20.489 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 1013## +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg" 1014``` 1015 1016If a given coordinate operation needs to be repeated, a good deal of time can be saved, as searches for coordinate operations take place once for `"SpatialPoints"` objects, but for `"SpatialLines"` and `"SpatialPolygons"` objects once for each "`Line"` or `"Polygon"` object: 1017 1018```{r, eval=run} 1019system.time(scot_LL3 <- spTransform(scot_BNG, CRS("+proj=longlat +datum=WGS84"), 1020 coordOp=get_last_coordOp())) 1021``` 1022 1023```{r} 1024## user system elapsed 1025## 0.066 0.002 0.070 1026``` 1027 1028 1029```{r, eval=run} 1030all.equal(b2, bbox(scot_LL3), scale=1) 1031``` 1032 1033```{r} 1034## [1] TRUE 1035``` 1036 1037The possibility of speeding up frequently used coordinate operations shows how listing candidate coordinate operations from a direct search lets us use functions from the previous section: 1038 1039 1040```{r, eval=run} 1041wkt_source_datum <- comment(slot(scot_BNG, "proj4string")) 1042wkt_target_datum <- comment(CRS("+proj=longlat +datum=WGS84")) 1043x <- list_coordOps(wkt_source_datum, wkt_target_datum) 1044system.time(scot_LL4 <- spTransform(scot_BNG, CRS("+proj=longlat +datum=WGS84"), 1045 coordOp=best_instantiable_coordOp(x))) 1046``` 1047 1048```{r} 1049## user system elapsed 1050## 0.066 0.002 0.069 1051``` 1052 1053 1054```{r, eval=run} 1055all.equal(b2, bbox(scot_LL4), scale=1) 1056``` 1057 1058```{r} 1059## [1] TRUE 1060``` 1061 1062## Axis swapping 1063 1064Luigi Ranghetti and Lorenzo Busetto provided [input](http://rpubs.com/ranghetti/sptransform-issue) with regard to axis swapping in **rgdal**. The task is to provide round-trip coordinate identity for four CRS. Here, we use just **sp**/**rgdal** to create the input objects: 1065 1066 1067```{r} 1068library(rgdal) 1069``` 1070 1071```{r} 1072packageVersion("rgdal") 1073rgdal_extSoftVersion() 1074``` 1075 1076 1077```{r} 1078pt0_lonlat <- SpatialPoints(matrix(c(10,46), nrow=1), proj4string= CRS("+init=epsg:4326")) 1079pt0_laea89 <- spTransform(pt0_lonlat, CRS("+init=epsg:3035")) 1080pt0_wgs32n <- spTransform(pt0_lonlat, CRS("+init=epsg:32632")) 1081pt0_psmerc <- spTransform(pt0_lonlat, CRS("+init=epsg:3857")) 1082``` 1083 1084```{r} 1085coordinates(pt0_lonlat) 1086``` 1087 1088```{r} 1089coordinates(pt0_laea89) 1090``` 1091 1092 1093```{r} 1094coordinates(pt0_wgs32n) 1095``` 1096 1097 1098```{r} 1099coordinates(pt0_psmerc) 1100``` 1101 1102We put the input (source) objects into a list: 1103 1104```{r} 1105from <- list(pt0_lonlat=pt0_lonlat, pt0_psmerc=pt0_psmerc, pt0_wgs32n=pt0_wgs32n, pt0_laea89=pt0_laea89) 1106names(from) 1107``` 1108 1109```{r} 1110sapply(from, proj4string) 1111``` 1112 1113and the target EPSG codes into a vector: 1114 1115```{r} 1116to <- c("4326", "3857", "32632", "3035") 1117``` 1118 1119From SVN revision 903, a global **rgdal** option can be set to enforce visualization order in `spTransform()` methods (the default value when the package is loaded is `TRUE`) when no `coordOp=` argument is given, and the coordinate operation has to be found automatically: 1120 1121```{r} 1122get_enforce_xy() 1123set_enforce_xy(FALSE) 1124get_enforce_xy() 1125``` 1126 1127```{r} 1128run <- TRUE 1129if (packageVersion("sp") < "1.3.3") run <- FALSE 1130if (packageVersion("rgdal") < "1.5.3") run <- FALSE 1131if (run && !rgdal::new_proj_and_gdal()) run <- FALSE 1132``` 1133 1134 1135Setting it to false gives failing round-trips for two of the four CRS: 1136 1137```{r, warning=FALSE, message=FALSE, eval=run} 1138out_EPSG_non_viz <- matrix(as.logical(NA), nrow=4, ncol=4) 1139colnames(out_EPSG_non_viz) <- names(from) 1140rownames(out_EPSG_non_viz) <- to 1141for (j in seq(along=from)) { 1142 for (i in seq(along=to)) { 1143 pt1 <- spTransform(from[[j]], CRS(SRS_string = paste0("EPSG:", to[i]))) 1144 out_EPSG_non_viz[i, j] <- isTRUE(all.equal(coordinates(from[[i]]), 1145 coordinates(pt1))) 1146 } 1147} 1148out_EPSG_non_viz 1149``` 1150 1151Resetting to the default (`TRUE`) value gives round-trip success: 1152 1153```{r} 1154set_enforce_xy(TRUE) 1155get_enforce_xy() 1156``` 1157 1158```{r, warning=FALSE, message=FALSE, eval=run} 1159out_EPSG <- matrix(as.logical(NA), nrow=4, ncol=4) 1160colnames(out_EPSG) <- names(from) 1161rownames(out_EPSG) <- to 1162for (j in seq(along=from)) { 1163 for (i in seq(along=to)) { 1164 pt1 <- spTransform(from[[j]], CRS(SRS_string = paste0("EPSG:", to[i]))) 1165 out_EPSG[i, j] <- isTRUE(all.equal(coordinates(from[[i]]), 1166 coordinates(pt1))) 1167 } 1168} 1169out_EPSG 1170``` 1171 1172If we instantiate the target CRS simply using a PROJ string with `+init=` rather than using the `CRS()` `SRS_string=` argument, it appears that visualization order is chosen anyway, whatever the status of `enforce_xy`: 1173 1174```{r} 1175get_enforce_xy() 1176set_enforce_xy(FALSE) 1177get_enforce_xy() 1178``` 1179 1180```{r, warning=FALSE, message=FALSE} 1181out_init_non_viz <- matrix(as.logical(NA), nrow=4, ncol=4) 1182colnames(out_init_non_viz) <- names(from) 1183rownames(out_init_non_viz) <- to 1184for (j in seq(along=from)) { 1185 for (i in seq(along=to)) { 1186 pt1 <- spTransform(from[[j]], CRS(paste0("+init=epsg:", to[i]))) 1187 out_init_non_viz[i, j] <- isTRUE(all.equal(coordinates(from[[i]]), 1188 coordinates(pt1))) 1189 } 1190} 1191out_init_non_viz 1192``` 1193 1194```{r} 1195set_enforce_xy(TRUE) 1196get_enforce_xy() 1197``` 1198 1199 1200```{r, warning=FALSE, message=FALSE} 1201out_init <- matrix(as.logical(NA), nrow=4, ncol=4) 1202colnames(out_init) <- names(from) 1203rownames(out_init) <- to 1204for (j in seq(along=from)) { 1205 for (i in seq(along=to)) { 1206 pt1 <- spTransform(from[[j]], CRS(paste0("+init=epsg:", to[i]))) 1207 out_init[i, j] <- isTRUE(all.equal(coordinates(from[[i]]), 1208 coordinates(pt1))) 1209 } 1210} 1211out_init 1212``` 1213 1214Finally, by listing candidate coordinate operations first, and choosing the best one, use can be made of the `list_coordOps()` `visualization_order=` argument, with default value `TRUE`; here it is given explicitly. This gives round trip success: 1215 1216```{r, warning=FALSE, message=FALSE, eval=run} 1217out_EPSG <- matrix(as.logical(NA), nrow=4, ncol=4) 1218colnames(out_EPSG) <- names(from) 1219rownames(out_EPSG) <- to 1220for (j in seq(along=from)) { 1221 for (i in seq(along=to)) { 1222 coo1 <- list_coordOps(comment(slot(from[[j]], "proj4string")), 1223 comment(CRS(SRS_string = paste0("EPSG:", to[i]))), 1224 visualization_order=TRUE) 1225 pt1 <- spTransform(from[[j]], CRS(SRS_string = paste0("EPSG:", to[i])), 1226 coordOp=best_instantiable_coordOp(coo1)) 1227 out_EPSG[i, j] <- isTRUE(all.equal(coordinates(from[[i]]), 1228 coordinates(pt1))) 1229 } 1230} 1231out_EPSG 1232``` 1233 1234The setting of the `list_coordOps()` `visualization_order=` argument overrides the global option value (as would setting `enforce_xy=` in calls to `spTransform()`): 1235 1236```{r, warning=FALSE, message=FALSE, eval=run} 1237out_EPSG_non_viz1 <- matrix(as.logical(NA), nrow=4, ncol=4) 1238colnames(out_EPSG_non_viz1) <- names(from) 1239rownames(out_EPSG_non_viz1) <- to 1240for (j in seq(along=from)) { 1241 for (i in seq(along=to)) { 1242 coo1 <- list_coordOps(comment(slot(from[[j]], "proj4string")), 1243 comment(CRS(SRS_string = paste0("EPSG:", to[i]))), 1244 visualization_order=FALSE) 1245 pt1 <- spTransform(from[[j]], CRS(SRS_string = paste0("EPSG:", to[i])), 1246 coordOp=best_instantiable_coordOp(coo1)) 1247 out_EPSG_non_viz1[i, j] <- isTRUE(all.equal(coordinates(from[[i]]), 1248 coordinates(pt1))) 1249 } 1250} 1251out_EPSG_non_viz1 1252``` 1253 1254 1255```{r, warning=FALSE, message=FALSE, eval=run} 1256out_EPSG_non_viz2 <- matrix(as.logical(NA), nrow=4, ncol=4) 1257colnames(out_EPSG_non_viz2) <- names(from) 1258rownames(out_EPSG_non_viz2) <- to 1259for (j in seq(along=from)) { 1260 for (i in seq(along=to)) { 1261 pt1 <- spTransform(from[[j]], CRS(SRS_string = paste0("EPSG:", to[i])), 1262 enforce_xy=FALSE) 1263 out_EPSG_non_viz2[i, j] <- isTRUE(all.equal(coordinates(from[[i]]), 1264 coordinates(pt1))) 1265 } 1266} 1267out_EPSG_non_viz2 1268``` 1269 1270## Adapting `project()` for PROJ 6+ 1271 1272It was originally intended to retire `project()`, because of the need to focus on `spTransform()`. Because it turned out that more packages than anticipated use `project()`, it has been adapted for the new setting. The declared projection will be accepted as a PROJ string, a WKT2 string, or similar, and new PROJ functions are used first to extract the geographical CRS from the given projected CRS, and a transformation found either forward from geographical to projected or inverse from projected to geographical. It now uses `proj_trans()` internally to convert single coordinates, so permitting the handling of out-of-domain coordinates. Note that no adaptation has been made for `legacy=FALSE` because as yet no Windows 32-bit build of PROJ or GDAL have been tried - `legacy=FALSE` uses the pre-PROJ6 compiled function `transform()`. The `verbose=` argument shows the chosen transformation pipeline: 1273 1274```{r} 1275ll <- structure(c(12.1823368669203, 11.9149630062421, 12.3186076188739, 127612.6207597184845, 12.9955172054652, 12.6316117692658, 12.4680041846297, 127712.4366882666609, NA, NA, -5.78993051516384, -5.03798674888479, 1278-4.60623015708619, -4.43802336997614, -4.78110320396188, -4.99127125409291, 1279-5.24836150474498, -5.68430388755925, NA, NA), .Dim = c(10L, 12802L), .Dimnames = list(NULL, c("longitude", "latitude"))) 1281try(xy0 <- project(ll, "+proj=moll", legacy=TRUE, verbose=TRUE)) 1282if (!exists("xy0")) xy0 <- structure(c(1217100.8468177, 1191302.229156, 12831232143.28841193, 1262546.27733232, 1299648.82357849, 1263011.18154638, 12841246343.17808186, 1242654.33986052, NA, NA, -715428.207551599, 1285-622613.577983058, -569301.605757784, -548528.530156422, -590895.949857199, 1286-616845.926397351, -648585.161643274, -702393.1160979, NA, NA), 1287.Dim = c(10L, 2L), .Dimnames = list(NULL, c("longitude", "latitude"))) 1288try(ll0 <- project(xy0, "+proj=moll", inv=TRUE, legacy=TRUE, verbose=TRUE)) 1289if (exists("ll0")) all.equal(ll, ll0) 1290``` 1291 1292It is possible to use the `coordOp=` argument to pass through a known transformation pipeline: 1293 1294```{r} 1295try(xy1 <- project(ll, "+proj=moll", legacy=TRUE, coordOp=paste("+proj=pipeline +step", 1296"+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=moll +lon_0=0 +x_0=0 +y_0=0 ellps=WGS84"))) 1297try(ll1 <- project(xy1, "+proj=moll", inv=TRUE, legacy=TRUE, coordOp=paste("+proj=pipeline +step", 1298"+inv +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg"))) 1299if (exists("ll1")) all.equal(ll, ll1) 1300``` 1301 1302```{r} 1303WKT <- CRS("+proj=moll") 1304try(xy2 <- project(ll, comment(WKT), legacy=TRUE, verbose=TRUE)) 1305try(ll2 <- project(xy1, comment(WKT), inv=TRUE, legacy=TRUE, verbose=TRUE)) 1306if (exists("ll2")) all.equal(ll, ll2) 1307``` 1308 1309 1310