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