R spatial
follows GDAL and PROJ development 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.
This 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.
For 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, https://github.com/r-spatial/sf/issues/1328 and many others.
We’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.
Because 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.
‘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 and Knudsen 2017).
Following the introduction of geodetic modules and pipelines in PROJ 5 (Knudsen and Evers 2017; Evers and Knudsen 2017), PROJ 6 moves further. Changes in the legacy PROJ representation and WGS84 transformation hub have been coordinated through the GDAL barn raising 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.
There are very useful postings on the PROJ mailing list from Martin Desruisseaux, first proposing clarifications and a follow-up including a summary:
- “Early binding” ≈ hub transformation technique.
- “Late binding” ≈ hub transformation technique NOT used, replaced by a more complex technique consisting in searching parameters in the EPSG database after the transformation context (source, target, epoch, area of interest) is known.
- The problem of hub transformation technique is independent of WGS84. It is caused by the fact that transformations to/from the hub are approximate. Any other hub we could invent in replacement of WGS84 will have the same problem, unless we can invent a hub for which transformations are exact (I think that if such hub existed, we would have already heard about it).
The solution proposed by ISO 19111 (in my understanding) is:
- Forget about hub (WGS84 or other), unless the simplicity of early-binding is considered more important than accuracy.
- Associating a CRS to a coordinate set (geometry or raster) is no longer sufficient. A {CRS, epoch} tuple must be associated. ISO 19111 calls this tuple “Coordinate metadata”. From a programmatic API point of view, this means that getCoordinateReferenceSystem() method in Geometry objects (for instance) needs to be replaced by a getCoordinateMetadata() method.
For users of North American spatial data, this ESRI news item gives a broad-brush picture of some of the motivations and oncoming changes.
The 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:
Warning 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.
(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).
This 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.
My 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.
## [1] '2.1.4'The "CRS" class definition is unchanged going forward to
maintain backward compatibility.
## Class "CRS" [package "sp"]
##
## Slots:
##
## Name: projargs
## Class: character## GDAL GDAL_with_GEOS PROJ sp EPSG
## "3.8.4" "TRUE" "9.4.0" "2.1-4" "v11.004"## [1] '1.6.7'The 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.
It 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.
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +ellps=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["unknown",
## DATUM["Unknown based on WGS 84 ellipsoid",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1],
## ID["EPSG",7030]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"):
## Discarded datum Unknown_based_on_WGS84_ellipsoid in CRS definition
## CRS arguments: +proj=longlat +ellps=WGS84 +no_defsIn 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).
## GEOGCRS["unknown", DATUM["Unknown based on WGS 84 ellipsoid",
## ELLIPSOID["WGS 84", 6378137, 298.257223563, LENGTHUNIT["metre", 1],
## ID["EPSG", 7030]]], PRIMEM["Greenwich", 0, ANGLEUNIT["degree",
## 0.0174532925199433], ID["EPSG", 8901]], CS[ellipsoidal, 2],
## AXIS["longitude", east, ORDER[1], ANGLEUNIT["degree",
## 0.0174532925199433, ID["EPSG", 9122]]], AXIS["latitude", north,
## ORDER[2], ANGLEUNIT["degree", 0.0174532925199433, ID["EPSG", 9122]]]]## [1] "GEOGCRS[\"unknown\", DATUM[\"Unknown based on WGS84 ellipsoid\",
## ELLIPSOID[\"WGS 84\", 6378137, 298.257223563, LENGTHUNIT[\"metre\", 1],
## ID[\"EPSG\", 7030]]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\",
## 0.0174532925199433], ID[\"EPSG\", 8901]], CS[ellipsoidal, 2], AXIS[\"longitude\",
## east, ORDER[1], ANGLEUNIT[\"degree\", 0.0174532925199433, ID[\"EPSG\", 9122]]],
## AXIS[\"latitude\", north, ORDER[2], ANGLEUNIT[\"degree\", 0.0174532925199433,
## ID[\"EPSG\", 9122]]]]"For 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.
A 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.
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]## GEOGCRS["unknown", DATUM["World Geodetic System 1984", ELLIPSOID["WGS
## 84", 6378137, 298.257223563, LENGTHUNIT["metre", 1]], ID["EPSG",
## 6326]], PRIMEM["Greenwich", 0, ANGLEUNIT["degree", 0.0174532925199433],
## ID["EPSG", 8901]], CS[ellipsoidal, 2], AXIS["longitude", east,
## ORDER[1], ANGLEUNIT["degree", 0.0174532925199433, ID["EPSG", 9122]]],
## AXIS["latitude", north, ORDER[2], ANGLEUNIT["degree",
## 0.0174532925199433, ID["EPSG", 9122]]]]## [1] "GEOGCRS[\"unknown\", DATUM[\"World Geodetic System 1984\",
## ELLIPSOID[\"WGS 84\", 6378137, 298.257223563, LENGTHUNIT[\"metre\", 1]],
## ID[\"EPSG\", 6326]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\",
## 0.0174532925199433], ID[\"EPSG\", 8901]], CS[ellipsoidal, 2], AXIS[\"longitude\",
## east, ORDER[1], ANGLEUNIT[\"degree\", 0.0174532925199433, ID[\"EPSG\", 9122]]],
## AXIS[\"latitude\", north, ORDER[2], ANGLEUNIT[\"degree\", 0.0174532925199433,
## ID[\"EPSG\", 9122]]]]"We know that +datum, +towgs84,
+nadgrids and +init are fragile, so we’ll try
one:
## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
## WKT2 2019 representation:
## BOUNDCRS[
## SOURCECRS[
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Geocentric translations (geog2D domain)",
## ID["EPSG",9603]],
## PARAMETER["X-axis translation",0,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",0,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",0,
## ID["EPSG",8607]]]]## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum WGS_1984 in CRS definition,
## but +towgs84= values preserved
## CRS arguments:
## +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defsThe comment here includes the +towgs84 parameters (three
of them), while the PROJ string gives seven.
## BOUNDCRS[ SOURCECRS[ GEOGCRS["unknown", DATUM["World Geodetic System
## 1984", ELLIPSOID["WGS 84", 6378137, 298.257223563, LENGTHUNIT["metre",
## 1]], ID["EPSG", 6326]], PRIMEM["Greenwich", 0, ANGLEUNIT["degree",
## 0.0174532925199433], ID["EPSG", 8901]], CS[ellipsoidal, 2],
## AXIS["longitude", east, ORDER[1], ANGLEUNIT["degree",
## 0.0174532925199433, ID["EPSG", 9122]]], AXIS["latitude", north,
## ORDER[2], ANGLEUNIT["degree", 0.0174532925199433, ID["EPSG", 9122]]]]],
## TARGETCRS[ GEOGCRS["WGS 84", DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84", 6378137, 298.257223563, LENGTHUNIT["metre", 1]]],
## PRIMEM["Greenwich", 0, ANGLEUNIT["degree", 0.0174532925199433]],
## CS[ellipsoidal, 2], AXIS["latitude", north, ORDER[1],
## ANGLEUNIT["degree", 0.0174532925199433]], AXIS["longitude", east,
## ORDER[2], ANGLEUNIT["degree", 0.0174532925199433]], ID["EPSG", 4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Geocentric translations (geog2D domain)", ID["EPSG", 9603]],
## PARAMETER["X-axis translation", 0, ID["EPSG", 8605]], PARAMETER["Y-axis
## translation", 0, ID["EPSG", 8606]], PARAMETER["Z-axis translation", 0,
## ID["EPSG", 8607]]]]## [1] "BOUNDCRS[SOURCECRS[GEOGCRS[\"unknown\", DATUM[\"World Geodetic System 1984\",
## ELLIPSOID[\"WGS 84\", 6378137, 298.257223563, LENGTHUNIT[\"metre\", 1]], ID[\"EPSG\",
## 6326]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\", 0.0174532925199433],
## ID[\"EPSG\", 8901]], CS[ellipsoidal, 2], AXIS[\"longitude\", east, ORDER[1],
## ANGLEUNIT[\"degree\", 0.0174532925199433, ID[\"EPSG\", 9122]]], AXIS[\"latitude\",
## north, ORDER[2], ANGLEUNIT[\"degree\", 0.0174532925199433, ID[\"EPSG\", 9122]]]]],
## TARGETCRS[GEOGCRS[\"WGS 84\", DATUM[\"World Geodetic System 1984\", ELLIPSOID[\"WGS
## 84\", 6378137, 298.257223563, LENGTHUNIT[\"metre\", 1]]], PRIMEM[\"Greenwich\", 0,
## ANGLEUNIT[\"degree\", 0.0174532925199433]], CS[ellipsoidal, 2], AXIS[\"latitude\",
## north, ORDER[1], ANGLEUNIT[\"degree\", 0.0174532925199433]], AXIS[\"longitude\", east,
## ORDER[2], ANGLEUNIT[\"degree\", 0.0174532925199433]], ID[\"EPSG\", 4326]]],
## ABRIDGEDTRANSFORMATION[\"Transformation from unknown to WGS84\", METHOD[\"Geocentric
## translations (geog2D domain)\", ID[\"EPSG\", 9603]], PARAMETER[\"X-axis translation\",
## 0, ID[\"EPSG\", 8605]], PARAMETER[\"Y-axis translation\", 0, ID[\"EPSG\", 8606]],
## PARAMETER[\"Z-axis translation\", 0, ID[\"EPSG\", 8607]]]]"The +init value is still accepted, but not repeated in
the output:
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)",
## ID["EPSG",1166]],
## MEMBER["World Geodetic System 1984 (G730)",
## ID["EPSG",1152]],
## MEMBER["World Geodetic System 1984 (G873)",
## ID["EPSG",1153]],
## MEMBER["World Geodetic System 1984 (G1150)",
## ID["EPSG",1154]],
## MEMBER["World Geodetic System 1984 (G1674)",
## ID["EPSG",1155]],
## MEMBER["World Geodetic System 1984 (G1762)",
## ID["EPSG",1156]],
## MEMBER["World Geodetic System 1984 (G2139)",
## ID["EPSG",1309]],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1],
## ID["EPSG",7030]],
## ENSEMBLEACCURACY[2.0],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## USAGE[
## SCOPE["unknown"],
## AREA["World."],
## BBOX[-90,-180,90,180]]]## GEOGCRS["WGS 84", ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)", ID["EPSG", 1166]],
## MEMBER["World Geodetic System 1984 (G730)", ID["EPSG", 1152]],
## MEMBER["World Geodetic System 1984 (G873)", ID["EPSG", 1153]],
## MEMBER["World Geodetic System 1984 (G1150)", ID["EPSG", 1154]],
## MEMBER["World Geodetic System 1984 (G1674)", ID["EPSG", 1155]],
## MEMBER["World Geodetic System 1984 (G1762)", ID["EPSG", 1156]],
## MEMBER["World Geodetic System 1984 (G2139)", ID["EPSG", 1309]],
## ELLIPSOID["WGS 84", 6378137, 298.257223563, LENGTHUNIT["metre", 1],
## ID["EPSG", 7030]], ENSEMBLEACCURACY[2.0], ID["EPSG", 6326]],
## PRIMEM["Greenwich", 0, ANGLEUNIT["degree", 0.0174532925199433],
## ID["EPSG", 8901]], CS[ellipsoidal, 2], AXIS["longitude", east,
## ORDER[1], ANGLEUNIT["degree", 0.0174532925199433, ID["EPSG", 9122]]],
## AXIS["latitude", north, ORDER[2], ANGLEUNIT["degree",
## 0.0174532925199433, ID["EPSG", 9122]]], USAGE[ SCOPE["unknown"],
## AREA["World."], BBOX[-90, -180, 90, 180]]]## [1] "GEOGCRS[\"WGS 84\", DATUM[\"World Geodetic System 1984\", ELLIPSOID[\"WGS 84\",
## 6378137, 298.257223563, LENGTHUNIT[\"metre\", 1]], ID[\"EPSG\", 6326]],
## PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\", 0.0174532925199433], ID[\"EPSG\",
## 8901]], CS[ellipsoidal, 2], AXIS[\"longitude\", east, ORDER[1], ANGLEUNIT[\"degree\",
## 0.0174532925199433, ID[\"EPSG\", 9122]]], AXIS[\"latitude\", north, ORDER[2],
## ANGLEUNIT[\"degree\", 0.0174532925199433, ID[\"EPSG\", 9122]]],
## USAGE[SCOPE[\"unknown\"], AREA[\"World\"], BBOX[-90, -180, 90, 180]]]"## PROJ4/GDAL2 CRS arguments:
## PROJ4/GDAL2 +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
## PROJ4/GDAL2 +towgs84=0,0,0An early warning of difficulties with discarded +datum
values came with the GB datum:
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"):
## Discarded datum OSGB_1936 in CRS definition
## CRS arguments:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defsHere, 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.
## [1] "PROJCRS[\"OSGB 1936 / British National Grid\", BASEGEOGCRS[\"OSGB 1936\",
## DATUM[\"OSGB 1936\", ELLIPSOID[\"Airy 1830\", 6377563.396, 299.3249646,
## LENGTHUNIT[\"metre\", 1]]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\",
## 0.0174532925199433]], ID[\"EPSG\", 4277]], CONVERSION[\"British National Grid\",
## METHOD[\"Transverse Mercator\", ID[\"EPSG\", 9807]], PARAMETER[\"Latitude of natural
## origin\", 49, ANGLEUNIT[\"degree\", 0.0174532925199433], ID[\"EPSG\", 8801]],
## PARAMETER[\"Longitude of natural origin\", -2, ANGLEUNIT[\"degree\",
## 0.0174532925199433], ID[\"EPSG\", 8802]], PARAMETER[\"Scale factor at natural
## origin\", 0.9996012717, SCALEUNIT[\"unity\", 1], ID[\"EPSG\", 8805]],
## PARAMETER[\"False easting\", 400000, LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8806]],
## PARAMETER[\"False northing\", -100000, LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8807]],
## ID[\"EPSG\", 19916]], CS[Cartesian, 2], AXIS[\"(E)\", east, ORDER[1],
## LENGTHUNIT[\"metre\", 1, ID[\"EPSG\", 9001]]], AXIS[\"(N)\", north, ORDER[2],
## LENGTHUNIT[\"metre\", 1, ID[\"EPSG\", 9001]]], USAGE[SCOPE[\"unknown\"], AREA[\"UK -
## 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]]]"In PROJ4/GDAL2, the +datum and +towgs84
values give much better transformation accuracy from the PROJ
string.
## PROJ4/GDAL2 CRS arguments:
## PROJ4/GDAL2 +init=epsg:27700 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717
## PROJ4/GDAL2 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs
## PROJ4/GDAL2 +ellps=airy
## PROJ4/GDAL2 +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894From 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():
## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["OSGB36 / British National Grid",
## BASEGEOGCRS["OSGB36",
## DATUM["Ordnance Survey of Great Britain 1936",
## ELLIPSOID["Airy 1830",6377563.396,299.3249646,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4277]],
## CONVERSION["British National Grid",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",49,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-2,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996012717,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",400000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",-100000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],
## BBOX[49.75,-9.01,61.01,2.01]],
## ID["EPSG",27700]]## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO"):
## Discarded datum OSGB_1936 in CRS definition
## CRS arguments:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs## PROJCRS["OSGB36 / British National Grid", BASEGEOGCRS["OSGB36",
## DATUM["Ordnance Survey of Great Britain 1936", ELLIPSOID["Airy 1830",
## 6377563.396, 299.3249646, LENGTHUNIT["metre", 1]]], PRIMEM["Greenwich",
## 0, ANGLEUNIT["degree", 0.0174532925199433]], ID["EPSG", 4277]],
## CONVERSION["British National Grid", METHOD["Transverse Mercator",
## ID["EPSG", 9807]], PARAMETER["Latitude of natural origin", 49,
## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG", 8801]],
## PARAMETER["Longitude of natural origin", -2, ANGLEUNIT["degree",
## 0.0174532925199433], ID["EPSG", 8802]], PARAMETER["Scale factor at
## natural origin", 0.9996012717, SCALEUNIT["unity", 1], ID["EPSG",
## 8805]], PARAMETER["False easting", 400000, LENGTHUNIT["metre", 1],
## ID["EPSG", 8806]], PARAMETER["False northing", -100000,
## LENGTHUNIT["metre", 1], ID["EPSG", 8807]]], CS[Cartesian, 2],
## AXIS["(E)", east, ORDER[1], LENGTHUNIT["metre", 1]], AXIS["(N)", north,
## ORDER[2], LENGTHUNIT["metre", 1]], USAGE[ SCOPE["Engineering survey,
## topographic mapping."], AREA["United Kingdom (UK) - offshore to
## boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great
## Britain (England, Wales and Scotland). Isle of Man onshore."],
## BBOX[49.75, -9.01, 61.01, 2.01]], ID["EPSG", 27700]]## [1] "PROJCRS[\"OSGB 1936 / British National Grid\", BASEGEOGCRS[\"OSGB 1936\",
## DATUM[\"OSGB 1936\", ELLIPSOID[\"Airy 1830\", 6377563.396, 299.3249646,
## LENGTHUNIT[\"metre\", 1]]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\",
## 0.0174532925199433]], ID[\"EPSG\", 4277]], CONVERSION[\"British National Grid\",
## METHOD[\"Transverse Mercator\", ID[\"EPSG\", 9807]], PARAMETER[\"Latitude of natural
## origin\", 49, ANGLEUNIT[\"degree\", 0.0174532925199433], ID[\"EPSG\", 8801]],
## PARAMETER[\"Longitude of natural origin\", -2, ANGLEUNIT[\"degree\",
## 0.0174532925199433], ID[\"EPSG\", 8802]], PARAMETER[\"Scale factor at natural
## origin\", 0.9996012717, SCALEUNIT[\"unity\", 1], ID[\"EPSG\", 8805]],
## PARAMETER[\"False easting\", 400000, LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8806]],
## PARAMETER[\"False northing\", -100000, LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8807]]],
## CS[Cartesian, 2], AXIS[\"(E)\", east, ORDER[1], LENGTHUNIT[\"metre\", 1]],
## AXIS[\"(N)\", north, ORDER[2], LENGTHUNIT[\"metre\", 1]], USAGE[SCOPE[\"unknown\"],
## AREA[\"UK - 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]], ID[\"EPSG\", 27700]]"We can also use a more detailed PROJ string, but without any improvement of the output PROJ representation; the comment is OK:
(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"))## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["OSGB 1936",
## ELLIPSOID["Airy 1830",6377563.396,299.3249646,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6277]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["unknown",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",49,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-2,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996012717,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",400000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",-100000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"):
## Discarded datum OSGB_1936 in CRS definition
## CRS arguments:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs## PROJCRS["unknown", BASEGEOGCRS["unknown", DATUM["OSGB 1936",
## ELLIPSOID["Airy 1830", 6377563.396, 299.3249646, LENGTHUNIT["metre",
## 1]], ID["EPSG", 6277]], PRIMEM["Greenwich", 0, ANGLEUNIT["degree",
## 0.0174532925199433], ID["EPSG", 8901]]], CONVERSION["unknown",
## METHOD["Transverse Mercator", ID["EPSG", 9807]], PARAMETER["Latitude of
## natural origin", 49, ANGLEUNIT["degree", 0.0174532925199433],
## ID["EPSG", 8801]], PARAMETER["Longitude of natural origin", -2,
## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG", 8802]],
## PARAMETER["Scale factor at natural origin", 0.9996012717,
## SCALEUNIT["unity", 1], ID["EPSG", 8805]], PARAMETER["False easting",
## 400000, LENGTHUNIT["metre", 1], ID["EPSG", 8806]], PARAMETER["False
## northing", -100000, LENGTHUNIT["metre", 1], ID["EPSG", 8807]]],
## CS[Cartesian, 2], AXIS["(E)", east, ORDER[1], LENGTHUNIT["metre", 1,
## ID["EPSG", 9001]]], AXIS["(N)", north, ORDER[2], LENGTHUNIT["metre", 1,
## ID["EPSG", 9001]]]]## [1] "PROJCRS[\"unknown\", BASEGEOGCRS[\"unknown\", DATUM[\"OSGB 1936\",
## ELLIPSOID[\"Airy 1830\", 6377563.396, 299.3249646, LENGTHUNIT[\"metre\", 1]],
## ID[\"EPSG\", 6277]], PRIMEM[\"Greenwich\", 0, ANGLEUNIT[\"degree\",
## 0.0174532925199433], ID[\"EPSG\", 8901]]], CONVERSION[\"unknown\", METHOD[\"Transverse
## Mercator\", ID[\"EPSG\", 9807]], PARAMETER[\"Latitude of natural origin\", 49,
## ANGLEUNIT[\"degree\", 0.0174532925199433], ID[\"EPSG\", 8801]], PARAMETER[\"Longitude
## of natural origin\", -2, ANGLEUNIT[\"degree\", 0.0174532925199433], ID[\"EPSG\",
## 8802]], PARAMETER[\"Scale factor at natural origin\", 0.9996012717,
## SCALEUNIT[\"unity\", 1], ID[\"EPSG\", 8805]], PARAMETER[\"False easting\", 400000,
## LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8806]], PARAMETER[\"False northing\", -100000,
## LENGTHUNIT[\"metre\", 1], ID[\"EPSG\", 8807]]], CS[Cartesian, 2], AXIS[\"(E)\", east,
## ORDER[1], LENGTHUNIT[\"metre\", 1, ID[\"EPSG\", 9001]]], AXIS[\"(N)\", north,
## ORDER[2], LENGTHUNIT[\"metre\", 1, ID[\"EPSG\", 9001]]]]"rgdal::showSRID() is quite versatile, so we can display
a multiline WKT string as well:
if (packageVersion("rgdal") >= "1.5.1") cat(rgdal::showSRID("+init=epsg:27700", format="WKT2_2018", multiline="YES", prefer_proj=FALSE), "\n")## PROJCRS["OSGB36 / British National Grid",
## BASEGEOGCRS["OSGB36",
## DATUM["Ordnance Survey of Great Britain 1936",
## ELLIPSOID["Airy 1830",6377563.396,299.3249646,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4277]],
## CONVERSION["British National Grid",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",49,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-2,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996012717,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",400000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",-100000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",19916]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## USAGE[
## SCOPE["unknown"],
## AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],
## BBOX[49.75,-9.01,61.01,2.01]]]## PROJCRS["OSGB 1936 / British National Grid",
## BASEGEOGCRS["OSGB 1936",
## DATUM["OSGB 1936",
## ELLIPSOID["Airy 1830",6377563.396,299.3249646,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4277]],
## CONVERSION["British National Grid",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",49,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-2,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996012717,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",400000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",-100000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",19916]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## USAGE[
## SCOPE["unknown"],
## AREA["UK - 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]]]So 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.
## rgdal: version: 1.5-1, (SVN revision 870)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.2, released 2019/10/28
## Path to GDAL shared files:
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 6.2.1, November 1st, 2019, [PJ_VERSION: 621]
## Path to PROJ.4 shared files: (autodetected)
## Linking to sp version: 1.3-3We 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.
## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs"Consequently, 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.
## Candidate coordinate operations found: 1
## Strict containment: FALSE
## Visualization order: TRUE
## Source: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs
## +type=crs
## Target: EPSG:4326
## Best instantiable operation has only ballpark accuracy
## Description: Inverse of unknown + Ballpark geographic offset from unknown to
## WGS 84 + axis order change (2D)
## Definition: +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2
## +k=0.9996012717 +x_0=400000 +y_0=-100000
## +ellps=airy +step +proj=unitconvert +xy_in=rad
## +xy_out=deg## Candidate coordinate operations found: 1
## Strict containment: FALSE
## Visualization order: TRUE
## Source: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs
## Target: EPSG:4326
## Best instantiable operation has only ballpark accuracy
## Description: Inverse of unknown + Ballpark geographic offset from unknown to
## WGS 84 + axis order change (2D)
## Definition: +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2
## +k=0.9996012717 +x_0=400000 +y_0=-100000
## +ellps=airy +step +proj=unitconvert +xy_in=rad
## +xy_out=degThe 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.
## Warning in best_instantiable_coordOp(x): Best instantiable operation has only
## ballpark accuracy## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=unitconvert +xy_in=rad +xy_out=deg"
## attr(,"description")
## [1] "Inverse of unknown + Ballpark geographic offset from unknown to WGS 84 + axis order change (2D)"## Warning in best_instantiable_coordOp(x): Best instantiable operation has only
## ballpark accuracy
## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717
## +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=unitconvert +xy_in=rad +xy_out=deg"
## attr(,"description")
## [1] "Inverse of unknown + Ballpark geographic offset from unknown to WGS 84 +
## axis order change (2D)"If we add back the discarded +datum= tag with a valid
value, we give the search process what it needs to find more coordinate
operations.
## Candidate coordinate operations found: 9
## Strict containment: FALSE
## Visualization order: TRUE
## Source: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs +type=crs
## +datum=OSGB36 +type=crs
## Target: EPSG:4326
## Best instantiable operation has accuracy: 2 m
## Description: Inverse of unknown + axis order change (2D) + OSGB36 to WGS 84
## (6) + axis order change (2D)
## Definition: +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2
## +k=0.9996012717 +x_0=400000 +y_0=-100000
## +ellps=airy +step +proj=push +v_3 +step +proj=cart
## +ellps=airy +step +proj=helmert +x=446.448
## +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842
## +s=-20.489 +convention=position_vector +step +inv
## +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step
## +proj=unitconvert +xy_in=rad +xy_out=deg
## Operation 8 is lacking 1 grid with accuracy 1 m
## Missing grid: uk_os_OSTN15_NTv2_OSGBtoETRS.tif
## URL: https://cdn.proj.org/uk_os_OSTN15_NTv2_OSGBtoETRS.tifNow 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.
## Candidate coordinate operations found: 7
## Strict containment: FALSE
## Visualization order: TRUE
## Source: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs
## +datum=OSGB36 +type=crs
## Target: EPSG:4326
## Best instantiable operation has accuracy: 2 m
## Description: Inverse of unknown + axis order change (2D) + OSGB 1936 to WGS
## 84 (6) + axis order change (2D)
## Definition: +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2
## +k=0.9996012717 +x_0=400000 +y_0=-100000
## +ellps=airy +step +proj=push +v_3 +step +proj=cart
## +ellps=airy +step +proj=helmert +x=446.448
## +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842
## +s=-20.489 +convention=position_vector +step +inv
## +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step
## +proj=unitconvert +xy_in=rad +xy_out=deg
## Operation 6 is lacking 1 grid with accuracy 1 m
## Missing grid: OSTN15_NTv2_OSGBtoETRS.gsb
## URL: https://download.osgeo.org/proj/proj-datumgrid-europe-latest.zipIf 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:
wkt_source_datum <- showSRID("EPSG:27700", "WKT2")
wkt_target_datum <- showSRID("EPSG:4326", "WKT2")
(x <- list_coordOps(wkt_source_datum, wkt_target_datum))## Candidate coordinate operations found: 9
## Strict containment: FALSE
## Visualization order: TRUE
## Source: PROJCRS["OSGB36 / British National Grid", BASEGEOGCRS["OSGB36",
## DATUM["Ordnance Survey of Great Britain 1936",
## ELLIPSOID["Airy 1830", 6377563.396, 299.3249646,
## LENGTHUNIT["metre", 1]]], PRIMEM["Greenwich", 0,
## ANGLEUNIT["degree", 0.0174532925199433]], ID["EPSG",
## 4277]], CONVERSION["British National Grid",
## METHOD["Transverse Mercator", ID["EPSG", 9807]],
## PARAMETER["Latitude of natural origin", 49,
## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG",
## 8801]], PARAMETER["Longitude of natural origin", -2,
## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG",
## 8802]], PARAMETER["Scale factor at natural origin",
## 0.9996012717, SCALEUNIT["unity", 1], ID["EPSG", 8805]],
## PARAMETER["False easting", 400000, LENGTHUNIT["metre",
## 1], ID["EPSG", 8806]], PARAMETER["False northing",
## -100000, LENGTHUNIT["metre", 1], ID["EPSG", 8807]]],
## CS[Cartesian, 2], AXIS["(E)", east, ORDER[1],
## LENGTHUNIT["metre", 1]], AXIS["(N)", north, ORDER[2],
## LENGTHUNIT["metre", 1]], USAGE[ SCOPE["Engineering
## survey, topographic mapping."], AREA["United Kingdom
## (UK) - offshore to boundary of UKCS within 49°45'N to
## 61°N and 9°W to 2°E; onshore Great Britain (England,
## Wales and Scotland). Isle of Man onshore."],
## BBOX[49.75, -9.01, 61.01, 2.01]], ID["EPSG", 27700]]
## Target:
## GEOGCRS["WGS 84 (with axis order normalized for
## visualization)", ENSEMBLE["World Geodetic System 1984
## ensemble", MEMBER["World Geodetic System 1984
## (Transit)", ID["EPSG", 1166]], MEMBER["World Geodetic
## System 1984 (G730)", ID["EPSG", 1152]], MEMBER["World
## Geodetic System 1984 (G873)", ID["EPSG", 1153]],
## MEMBER["World Geodetic System 1984 (G1150)", ID["EPSG",
## 1154]], MEMBER["World Geodetic System 1984 (G1674)",
## ID["EPSG", 1155]], MEMBER["World Geodetic System 1984
## (G1762)", ID["EPSG", 1156]], MEMBER["World Geodetic
## System 1984 (G2139)", ID["EPSG", 1309]], ELLIPSOID["WGS
## 84", 6378137, 298.257223563, LENGTHUNIT["metre", 1],
## ID["EPSG", 7030]], ENSEMBLEACCURACY[2.0], ID["EPSG",
## 6326]], PRIMEM["Greenwich", 0, ANGLEUNIT["degree",
## 0.0174532925199433], ID["EPSG", 8901]], CS[ellipsoidal,
## 2], AXIS["geodetic longitude (Lon)", east, ORDER[1],
## ANGLEUNIT["degree", 0.0174532925199433, ID["EPSG",
## 9122]]], AXIS["geodetic latitude (Lat)", north,
## ORDER[2], ANGLEUNIT["degree", 0.0174532925199433,
## ID["EPSG", 9122]]], USAGE[ SCOPE["Horizontal component
## of 3D system."], AREA["World."], BBOX[-90, -180, 90,
## 180]], REMARK["Axis order reversed compared to
## EPSG:4326"]]
## Best instantiable operation has accuracy: 2 m
## Description: Inverse of British National Grid + OSGB36 to WGS 84 (6) + axis
## order change (2D)
## Definition: +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2
## +k=0.9996012717 +x_0=400000 +y_0=-100000
## +ellps=airy +step +proj=push +v_3 +step +proj=cart
## +ellps=airy +step +proj=helmert +x=446.448
## +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842
## +s=-20.489 +convention=position_vector +step +inv
## +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step
## +proj=unitconvert +xy_in=rad +xy_out=deg
## Operation 8 is lacking 1 grid with accuracy 1 m
## Missing grid: uk_os_OSTN15_NTv2_OSGBtoETRS.tif
## URL: https://cdn.proj.org/uk_os_OSTN15_NTv2_OSGBtoETRS.tif## Candidate coordinate operations found: 6
## Strict containment: FALSE
## Visualization order: TRUE
## Source: PROJCRS["OSGB 1936 / British National Grid", BASEGEOGCRS["OSGB
## 1936", DATUM["OSGB 1936", ELLIPSOID["Airy 1830",
## 6377563.396, 299.3249646, LENGTHUNIT["metre", 1]]],
## PRIMEM["Greenwich", 0, ANGLEUNIT["degree",
## 0.0174532925199433]], ID["EPSG", 4277]],
## CONVERSION["British National Grid", METHOD["Transverse
## Mercator", ID["EPSG", 9807]], PARAMETER["Latitude of
## natural origin", 49, ANGLEUNIT["degree",
## 0.0174532925199433], ID["EPSG", 8801]],
## PARAMETER["Longitude of natural origin", -2,
## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG",
## 8802]], PARAMETER["Scale factor at natural origin",
## 0.9996012717, SCALEUNIT["unity", 1], ID["EPSG", 8805]],
## PARAMETER["False easting", 400000, LENGTHUNIT["metre",
## 1], ID["EPSG", 8806]], PARAMETER["False northing",
## -100000, LENGTHUNIT["metre", 1], ID["EPSG", 8807]]],
## CS[Cartesian, 2], AXIS["(E)", east, ORDER[1],
## LENGTHUNIT["metre", 1]], AXIS["(N)", north, ORDER[2],
## LENGTHUNIT["metre", 1]], USAGE[SCOPE["unknown"],
## AREA["UK - 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]],
## ID["EPSG", 27700]]
## Target: GEOGCRS["WGS 84", DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84", 6378137, 298.257223563,
## LENGTHUNIT["metre", 1]]], PRIMEM["Greenwich", 0,
## ANGLEUNIT["degree", 0.0174532925199433]],
## CS[ellipsoidal, 2], AXIS["geodetic latitude (Lat)",
## north, ORDER[1], ANGLEUNIT["degree",
## 0.0174532925199433]], AXIS["geodetic longitude (Lon)",
## east, ORDER[2], ANGLEUNIT["degree",
## 0.0174532925199433]], USAGE[SCOPE["unknown"],
## AREA["World"], BBOX[-90, -180, 90, 180]], ID["EPSG",
## 4326]]
## Best instantiable operation has accuracy: 2 m
## Description: Inverse of British National Grid + OSGB 1936 to WGS 84 (6) +
## axis order change (2D)
## Definition: +proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2
## +k=0.9996012717 +x_0=400000 +y_0=-100000
## +ellps=airy +step +proj=push +v_3 +step +proj=cart
## +ellps=airy +step +proj=helmert +x=446.448
## +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842
## +s=-20.489 +convention=position_vector +step +inv
## +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step
## +proj=unitconvert +xy_in=rad +xy_out=deg
## Operation 6 is lacking 1 grid with accuracy 1 m
## Missing grid: OSTN15_NTv2_OSGBtoETRS.gsb
## URL: https://download.osgeo.org/proj/proj-datumgrid-europe-latest.zip## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=push +v_3 +step +proj=cart +ellps=airy +step +proj=helmert +x=446.448 +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842 +s=-20.489 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg"
## attr(,"description")
## [1] "Inverse of British National Grid + OSGB36 to WGS 84 (6) + axis order change (2D)"## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717
## +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=push +v_3 +step +proj=cart
## +ellps=airy +step +proj=helmert +x=446.448 +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247
## +rz=0.842 +s=-20.489 +convention=position_vector +step +inv +proj=cart +ellps=WGS84
## +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg"
## attr(,"description")
## [1] "Inverse of British National Grid + OSGB 1936 to WGS 84 (6) + axis order change (2D)"Turning 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.
## Warning in showSRID("EPSG:22525", "PROJ"): Discarded datum Corrego_Alegre_1970-72 in CRS definition,
## but +towgs84= values preservedThis meant that while only ballpark accuracy is reported, a Helmert transformation is carried out:
## Candidate coordinate operations found: 1
## Strict containment: FALSE
## Visualization order: TRUE
## Source: +proj=utm +zone=25 +south +ellps=intl +units=m +no_defs
## +type=crs +type=crs
## Target: EPSG:31985
## Best instantiable operation has only ballpark accuracy
## Description: Inverse of UTM zone 25S + Ballpark geographic offset from
## unknown to SIRGAS 2000 + UTM zone 25S
## Definition: +proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl
## +step +proj=utm +zone=25 +south +ellps=GRS80## Candidate coordinate operations found: 1
## Strict containment: FALSE
## Visualization order: TRUE
## Source: +proj=utm +zone=25 +south +ellps=intl
## +towgs84=-205.57,168.77,-4.12,0,0,0,0 +units=m +no_defs
## +type=crs
## Target: EPSG:31985
## Best instantiable operation has only ballpark accuracy
## Description: Inverse of UTM zone 25S + Transformation from unknown to WGS84
## + Inverse of SIRGAS 2000 to WGS 84 (1) + UTM zone
## 25S
## Definition: +proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl
## +step +proj=push +v_3 +step +proj=cart +ellps=intl
## +step +proj=helmert +x=-205.57 +y=168.77 +z=-4.12
## +rx=0 +ry=0 +rz=0 +s=0 +convention=position_vector
## +step +inv +proj=cart +ellps=GRS80 +step +proj=pop
## +v_3 +step +proj=utm +zone=25 +south +ellps=GRS80## Warning in best_instantiable_coordOp(x): Best instantiable operation has only
## ballpark accuracy## [1] "+proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl +step +proj=utm +zone=25 +south +ellps=GRS80"
## attr(,"description")
## [1] "Inverse of UTM zone 25S + Ballpark geographic offset from unknown to SIRGAS 2000 + UTM zone 25S"## Warning in best_instantiable_coordOp(x): Best instantiable operation has only
## ballpark accuracy
## [1] "+proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl +step
## +proj=push +v_3 +step +proj=cart +ellps=intl +step +proj=helmert +x=-205.57
## +y=168.77 +z=-4.12 +rx=0 +ry=0 +rz=0 +s=0 +convention=position_vector +step +inv
## +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=utm +zone=25 +south
## +ellps=GRS80"
## attr(,"description")
## [1] "Inverse of UTM zone 25S + Transformation from unknown to WGS84 + Inverse of
## SIRGAS 2000 to WGS 84 (1) + UTM zone 25S"If 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:
wkt_source_datum <- showSRID("EPSG:22525", "WKT2")
wkt_target_datum <- showSRID("EPSG:31985", "WKT2")
(x <- list_coordOps(wkt_source_datum, wkt_target_datum))## Candidate coordinate operations found: 2
## Strict containment: FALSE
## Visualization order: TRUE
## Source: PROJCRS["Corrego Alegre 1970-72 / UTM zone 25S",
## BASEGEOGCRS["Corrego Alegre 1970-72", DATUM["Corrego
## Alegre 1970-72", ELLIPSOID["International 1924",
## 6378388, 297, LENGTHUNIT["metre", 1]]],
## PRIMEM["Greenwich", 0, ANGLEUNIT["degree",
## 0.0174532925199433]], ID["EPSG", 4225]],
## CONVERSION["UTM zone 25S", METHOD["Transverse
## Mercator", ID["EPSG", 9807]], PARAMETER["Latitude of
## natural origin", 0, ANGLEUNIT["degree",
## 0.0174532925199433], ID["EPSG", 8801]],
## PARAMETER["Longitude of natural origin", -33,
## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG",
## 8802]], PARAMETER["Scale factor at natural origin",
## 0.9996, SCALEUNIT["unity", 1], ID["EPSG", 8805]],
## PARAMETER["False easting", 500000, LENGTHUNIT["metre",
## 1], ID["EPSG", 8806]], PARAMETER["False northing",
## 10000000, LENGTHUNIT["metre", 1], ID["EPSG", 8807]]],
## CS[Cartesian, 2], AXIS["(E)", east, ORDER[1],
## LENGTHUNIT["metre", 1]], AXIS["(N)", north, ORDER[2],
## LENGTHUNIT["metre", 1]], USAGE[ SCOPE["Engineering
## survey, topographic mapping."], AREA["Brazil - onshore
## east of 36°W ."], BBOX[-10.1, -36, -4.99, -34.74]],
## ID["EPSG", 22525]]
## Target:
## PROJCRS["SIRGAS 2000 / UTM zone 25S", BASEGEOGCRS["SIRGAS
## 2000", DATUM["Sistema de Referencia Geocentrico para
## las AmericaS 2000", ELLIPSOID["GRS 1980", 6378137,
## 298.257222101, LENGTHUNIT["metre", 1]]],
## PRIMEM["Greenwich", 0, ANGLEUNIT["degree",
## 0.0174532925199433]], ID["EPSG", 4674]],
## CONVERSION["UTM zone 25S", METHOD["Transverse
## Mercator", ID["EPSG", 9807]], PARAMETER["Latitude of
## natural origin", 0, ANGLEUNIT["degree",
## 0.0174532925199433], ID["EPSG", 8801]],
## PARAMETER["Longitude of natural origin", -33,
## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG",
## 8802]], PARAMETER["Scale factor at natural origin",
## 0.9996, SCALEUNIT["unity", 1], ID["EPSG", 8805]],
## PARAMETER["False easting", 500000, LENGTHUNIT["metre",
## 1], ID["EPSG", 8806]], PARAMETER["False northing",
## 10000000, LENGTHUNIT["metre", 1], ID["EPSG", 8807]]],
## CS[Cartesian, 2], AXIS["(E)", east, ORDER[1],
## LENGTHUNIT["metre", 1]], AXIS["(N)", north, ORDER[2],
## LENGTHUNIT["metre", 1]], USAGE[ SCOPE["Engineering
## survey, topographic mapping."], AREA["Brazil - between
## 36°W and 30°W, northern and southern hemispheres,
## onshore and offshore."], BBOX[-23.8, -36, 4.19,
## -29.99]], ID["EPSG", 31985]]
## Best instantiable operation has accuracy: 5 m
## Description: Inverse of UTM zone 25S + Corrego Alegre 1970-72 to SIRGAS 2000
## (2) + UTM zone 25S
## Definition: +proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl
## +step +proj=push +v_3 +step +proj=cart +ellps=intl
## +step +proj=helmert +x=-206.05 +y=168.28 +z=-3.82
## +step +inv +proj=cart +ellps=GRS80 +step +proj=pop
## +v_3 +step +proj=utm +zone=25 +south +ellps=GRS80
## Operation 2 is lacking 1 grid with accuracy 2 m
## Missing grid: br_ibge_CA7072_003.tif
## URL: https://cdn.proj.org/br_ibge_CA7072_003.tif## Candidate coordinate operations found: 1
## Strict containment: FALSE
## Visualization order: TRUE
## Source: BOUNDCRS[SOURCECRS[PROJCRS["Corrego Alegre 1970-72 / UTM zone
## 25S", BASEGEOGCRS["Corrego Alegre 1970-72",
## DATUM["Corrego Alegre 1970-72",
## ELLIPSOID["International 1924", 6378388, 297,
## LENGTHUNIT["metre", 1]]], PRIMEM["Greenwich", 0,
## ANGLEUNIT["degree", 0.0174532925199433]], ID["EPSG",
## 4225]], CONVERSION["UTM zone 25S", METHOD["Transverse
## Mercator", ID["EPSG", 9807]], PARAMETER["Latitude of
## natural origin", 0, ANGLEUNIT["degree",
## 0.0174532925199433], ID["EPSG", 8801]],
## PARAMETER["Longitude of natural origin", -33,
## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG",
## 8802]], PARAMETER["Scale factor at natural origin",
## 0.9996, SCALEUNIT["unity", 1], ID["EPSG", 8805]],
## PARAMETER["False easting", 500000, LENGTHUNIT["metre",
## 1], ID["EPSG", 8806]], PARAMETER["False northing",
## 10000000, LENGTHUNIT["metre", 1], ID["EPSG", 8807]]],
## CS[Cartesian, 2], AXIS["(E)", east, ORDER[1],
## LENGTHUNIT["metre", 1]], AXIS["(N)", north, ORDER[2],
## LENGTHUNIT["metre", 1]], USAGE[SCOPE["unknown"],
## AREA["Brazil - east of 36°W onshore"], BBOX[-10.1, -36,
## -4.99, -34.74]], ID["EPSG", 22525]]],
## TARGETCRS[GEOGCRS["WGS 84", DATUM["World Geodetic
## System 1984", ELLIPSOID["WGS 84", 6378137,
## 298.257223563, LENGTHUNIT["metre", 1]]],
## PRIMEM["Greenwich", 0, ANGLEUNIT["degree",
## 0.0174532925199433]], CS[ellipsoidal, 2],
## AXIS["latitude", north, ORDER[1], ANGLEUNIT["degree",
## 0.0174532925199433]], AXIS["longitude", east, ORDER[2],
## ANGLEUNIT["degree", 0.0174532925199433]], ID["EPSG",
## 4326]]], ABRIDGEDTRANSFORMATION["Corrego Alegre 1970-72
## to WGS 84 (3)", VERSION["PBS-Bra 1983"],
## METHOD["Geocentric translations (geog2D domain)",
## ID["EPSG", 9603]], PARAMETER["X-axis translation",
## -205.57, ID["EPSG", 8605]], PARAMETER["Y-axis
## translation", 168.77, ID["EPSG", 8606]],
## PARAMETER["Z-axis translation", -4.12, ID["EPSG",
## 8607]], USAGE[SCOPE["Medium and small scale mapping."],
## AREA["Brazil - Corrego Alegre 1970-1972"], BBOX[-33.78,
## -58.16, -2.68, -34.74]], ID["EPSG", 6192],
## REMARK["Formed by concatenation of tfms codes 6191 and
## 1877. Used by Petrobras and ANP until February 2005
## when replaced by Corrego Alegre 1970-72 to WGS 84 (4)
## (tfm code 6194)."]]]
## Target:
## BOUNDCRS[SOURCECRS[PROJCRS["SIRGAS 2000 / UTM zone 25S",
## BASEGEOGCRS["SIRGAS 2000", DATUM["Sistema de Referencia
## Geocentrico para las AmericaS 2000", ELLIPSOID["GRS
## 1980", 6378137, 298.257222101, LENGTHUNIT["metre",
## 1]]], PRIMEM["Greenwich", 0, ANGLEUNIT["degree",
## 0.0174532925199433]], ID["EPSG", 4674]],
## CONVERSION["UTM zone 25S", METHOD["Transverse
## Mercator", ID["EPSG", 9807]], PARAMETER["Latitude of
## natural origin", 0, ANGLEUNIT["degree",
## 0.0174532925199433], ID["EPSG", 8801]],
## PARAMETER["Longitude of natural origin", -33,
## ANGLEUNIT["degree", 0.0174532925199433], ID["EPSG",
## 8802]], PARAMETER["Scale factor at natural origin",
## 0.9996, SCALEUNIT["unity", 1], ID["EPSG", 8805]],
## PARAMETER["False easting", 500000, LENGTHUNIT["metre",
## 1], ID["EPSG", 8806]], PARAMETER["False northing",
## 10000000, LENGTHUNIT["metre", 1], ID["EPSG", 8807]]],
## CS[Cartesian, 2], AXIS["(E)", east, ORDER[1],
## LENGTHUNIT["metre", 1]], AXIS["(N)", north, ORDER[2],
## LENGTHUNIT["metre", 1]], USAGE[SCOPE["unknown"],
## AREA["Brazil - 36°W to 30°W"], BBOX[-23.8, -36, 4.19,
## -29.99]], ID["EPSG", 31985]]], TARGETCRS[GEOGCRS["WGS
## 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS
## 84", 6378137, 298.257223563, LENGTHUNIT["metre", 1]]],
## PRIMEM["Greenwich", 0, ANGLEUNIT["degree",
## 0.0174532925199433]], CS[ellipsoidal, 2],
## AXIS["latitude", north, ORDER[1], ANGLEUNIT["degree",
## 0.0174532925199433]], AXIS["longitude", east, ORDER[2],
## ANGLEUNIT["degree", 0.0174532925199433]], ID["EPSG",
## 4326]]], ABRIDGEDTRANSFORMATION["SIRGAS 2000 to WGS 84
## (1)", VERSION["OGP-C&S America"], METHOD["Geocentric
## translations (geog2D domain)", ID["EPSG", 9603]],
## PARAMETER["X-axis translation", 0, ID["EPSG", 8605]],
## PARAMETER["Y-axis translation", 0, ID["EPSG", 8606]],
## PARAMETER["Z-axis translation", 0, ID["EPSG", 8607]],
## USAGE[SCOPE["Accuracy 1m."], AREA["Latin America -
## SIRGAS 2000 by country"], BBOX[-59.87, -122.19, 32.72,
## -25.28]], ID["EPSG", 15894]]]
## Best instantiable operation has only ballpark accuracy
## Description: Inverse of UTM zone 25S + Corrego Alegre 1970-72 to WGS 84 (3)
## + Inverse of SIRGAS 2000 to WGS 84 (1) + UTM zone
## 25S
## Definition: +proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl
## +step +proj=push +v_3 +step +proj=cart +ellps=intl
## +step +proj=helmert +x=-205.57 +y=168.77 +z=-4.12
## +step +inv +proj=cart +ellps=GRS80 +step +proj=pop
## +v_3 +step +proj=utm +zone=25 +south +ellps=GRS80If 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.
## Candidate coordinate operations found: 2
## Strict containment: FALSE
## Visualization order: TRUE
## Source: EPSG:22525
## Target: EPSG:31985
## Best instantiable operation has accuracy: 5 m
## Description: Inverse of UTM zone 25S + Corrego Alegre 1970-72 to SIRGAS 2000
## (2) + UTM zone 25S
## Definition: +proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl
## +step +proj=push +v_3 +step +proj=cart +ellps=intl
## +step +proj=helmert +x=-206.05 +y=168.28 +z=-3.82
## +step +inv +proj=cart +ellps=GRS80 +step +proj=pop
## +v_3 +step +proj=utm +zone=25 +south +ellps=GRS80
## Operation 2 is lacking 1 grid with accuracy 2 m
## Missing grid: br_ibge_CA7072_003.tif
## URL: https://cdn.proj.org/br_ibge_CA7072_003.tif## Candidate coordinate operations found: 2
## Strict containment: FALSE
## Visualization order: TRUE
## Source: EPSG:22525
## Target: EPSG:31985
## Best instantiable operation has accuracy: 5 m
## Description: Inverse of UTM zone 25S + Corrego Alegre 1970-72 to SIRGAS 2000
## (2) + UTM zone 25S
## Definition: +proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl
## +step +proj=push +v_3 +step +proj=cart +ellps=intl
## +step +proj=helmert +x=-206.05 +y=168.28 +z=-3.82
## +step +inv +proj=cart +ellps=GRS80 +step +proj=pop
## +v_3 +step +proj=utm +zone=25 +south +ellps=GRS80
## Operation 2 is lacking 1 grid with accuracy 2 m
## Missing grid: CA7072_003.gsb## [1] "+proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl +step +proj=push +v_3 +step +proj=cart +ellps=intl +step +proj=helmert +x=-206.05 +y=168.28 +z=-3.82 +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=utm +zone=25 +south +ellps=GRS80"
## attr(,"description")
## [1] "Inverse of UTM zone 25S + Corrego Alegre 1970-72 to SIRGAS 2000 (2) + UTM zone 25S"## [1] "+proj=pipeline +step +inv +proj=utm +zone=25 +south +ellps=intl +step
## +proj=push +v_3 +step +proj=cart +ellps=intl +step +proj=helmert +x=-206.05
## +y=168.28 +z=-3.82 +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step
## +proj=utm +zone=25 +south +ellps=GRS80"
## attr(,"description")
## [1] "Inverse of UTM zone 25S + Corrego Alegre 1970-72 to SIRGAS 2000 (2) +
## UTM zone 25S"The 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):
## Warning in readOGR(system.file("vectors", package = "rgdal"), "scot_BNG"): OGR
## support is provided by the sf and terra packages among others## Warning in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv =
## use_iconv, : OGR support is provided by the sf and terra packages among others## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : OGR support is provided by the sf and terra packages among others## Warning in ogrListLayers(dsn): OGR support is provided by the sf and terra
## packages among others## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others## OGR data source with driver: ESRI Shapefile
## Source: "/tmp/RtmpMk1Iph/Rinstf89c9e4092/rgdal/vectors", layer: "scot_BNG"
## with 56 features
## It has 13 fields## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS): Discarded datum OSGB_1936 in CRS definition: +proj=tmerc +lat_0=49
## +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
## OGR data source with driver: ESRI Shapefile
## Source: "/home/rsb/lib/r_libs/rgdal/vectors", layer: "scot_BNG"
## with 56 features
## It has 13 fields
## Warning in OGRSpatialRef(dsn = dsn, layer = layer, morphFromESRI =
## morphFromESRI, : Discarded datum OSGB_1936 in CRS definition: +proj=tmerc
## +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy
## +units=m +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum Unknown_based_on_Airy_1830_ellipsoid in CRS definitionOn 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:
## [1] TRUE(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 aberration is unknown.
set_transform_wkt_comment(FALSE)
scot_LL <- spTransform(scot_BNG, CRS("+proj=longlat +datum=WGS84"))## Warning in project(t(bbox(obj))[, 1:2], tg, inv = TRUE, use_aoi = FALSE): PROJ
## support is provided by the sf and terra packages among others## min max
## x -8.621387 -0.753056
## y 54.626555 60.843843We might fake the PROJ string by extracting the "CRS"
object:
## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["OSGB36 / British National Grid",
## BASEGEOGCRS["OSGB36",
## DATUM["Ordnance Survey of Great Britain 1936",
## ELLIPSOID["Airy 1830",6377563.396,299.3249646,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4277]],
## CONVERSION["British National Grid",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",49,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-2,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996012717,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",400000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",-100000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],
## BBOX[49.75,-9.01,61.01,2.01]],
## ID["EPSG",27700]]## CRS arguments:
## +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000
## +y_0=-100000 +ellps=airy +units=m +no_defs## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000
## +ellps=airy +units=m +no_defs"and updating the PROJ string in-place if we know a sensible incantation:
slot(crs, "projargs") <- paste0(slot(crs, "projargs"), " +datum=OSGB36")
slot(scot_BNG, "proj4string") <- crs
scot_LL1 <- spTransform(scot_BNG, CRS("+proj=longlat +datum=WGS84"))## Warning in project(t(bbox(obj))[, 1:2], tg, inv = TRUE, use_aoi = FALSE): PROJ
## support is provided by the sf and terra packages among others## min max
## x -8.622158 -0.7550709
## y 54.626633 60.8432318The two bounding boxes differ a good deal:
## [1] "Mean absolute difference: 0.0008689015"The SW corner of Scotland differs by about 50m, the NE corner by almost 130m:
## [1] 50.56342 128.98675Re-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:
## Warning in readOGR(system.file("vectors", package = "rgdal"), "scot_BNG"): OGR
## support is provided by the sf and terra packages among others## Warning in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv =
## use_iconv, : OGR support is provided by the sf and terra packages among others## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : OGR support is provided by the sf and terra packages among others## Warning in ogrListLayers(dsn): OGR support is provided by the sf and terra
## packages among others## Warning in ogrFIDs(dsn = dsn, layer = layer): OGR support is provided by the sf
## and terra packages among others## OGR data source with driver: ESRI Shapefile
## Source: "/tmp/RtmpMk1Iph/Rinstf89c9e4092/rgdal/vectors", layer: "scot_BNG"
## with 56 features
## It has 13 fields## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS): Discarded datum OSGB_1936 in CRS definition: +proj=tmerc +lat_0=49
## +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
## OGR data source with driver: ESRI Shapefile
## Source: "/home/rsb/lib/r_libs/rgdal/vectors", layer: "scot_BNG"
## with 56 features
## It has 13 fields
## Warning in OGRSpatialRef(dsn = dsn, layer = layer, morphFromESRI =
## morphFromESRI, : Discarded datum OSGB_1936 in CRS definition: +proj=tmerc
## +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy
## +units=m +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum Unknown_based_on_Airy_1830_ellipsoid in CRS definitionIn rgdal 1.5-2, the coordinate operation lookup is done once for the first matrix of coordinates, and passed then on to subsequent transformations:
## Warning in project(t(bbox(obj))[, 1:2], tg, inv = TRUE, use_aoi = FALSE): PROJ
## support is provided by the sf and terra packages among others## user system elapsed
## 0.071 0.000 0.071Timings for rgdal 1.5-1, when look-up was repeated for each matrix of coordinates:
## min max
## x -8.622158 -0.7550709
## y 54.626633 60.8432318## [1] TRUEThe coordinate operation last used may be retrieved with (default empty):
## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=push +v_3 +step +proj=cart +ellps=airy +step +proj=helmert +x=446.448 +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247 +rz=0.842 +s=-20.489 +convention=position_vector +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg"## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717
## +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=push +v_3 +step +proj=cart
## +ellps=airy +step +proj=helmert +x=446.448 +y=-125.157 +z=542.06 +rx=0.15 +ry=0.247
## +rz=0.842 +s=-20.489 +convention=position_vector +step +inv +proj=cart +ellps=WGS84
## +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg"If 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:
system.time(scot_LL3 <- spTransform(scot_BNG, CRS("+proj=longlat +datum=WGS84"),
coordOp=get_last_coordOp()))## Warning in project(t(bbox(obj))[, 1:2], tg, inv = TRUE, use_aoi = FALSE): PROJ
## support is provided by the sf and terra packages among others## user system elapsed
## 0.048 0.000 0.049## [1] TRUEThe 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:
wkt_source_datum <- comment(slot(scot_BNG, "proj4string"))
wkt_target_datum <- comment(CRS("+proj=longlat +datum=WGS84"))
x <- list_coordOps(wkt_source_datum, wkt_target_datum)
system.time(scot_LL4 <- spTransform(scot_BNG, CRS("+proj=longlat +datum=WGS84"),
coordOp=best_instantiable_coordOp(x)))## Warning in project(t(bbox(obj))[, 1:2], tg, inv = TRUE, use_aoi = FALSE): PROJ
## support is provided by the sf and terra packages among others## user system elapsed
## 0.043 0.004 0.047## [1] TRUELuigi Ranghetti and Lorenzo Busetto provided input 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:
## [1] '1.6.7'## GDAL GDAL_with_GEOS PROJ sp EPSG
## "3.8.4" "TRUE" "9.4.0" "2.1-4" "v11.004"pt0_lonlat <- SpatialPoints(matrix(c(10,46), nrow=1), proj4string= CRS("+init=epsg:4326"))
pt0_laea89 <- spTransform(pt0_lonlat, CRS("+init=epsg:3035"))
pt0_wgs32n <- spTransform(pt0_lonlat, CRS("+init=epsg:32632"))
pt0_psmerc <- spTransform(pt0_lonlat, CRS("+init=epsg:3857"))## coords.x1 coords.x2
## [1,] 10 46## coords.x1 coords.x2
## [1,] 4321000 2543009## coords.x1 coords.x2
## [1,] 577432.2 5094534## coords.x1 coords.x2
## [1,] 1113195 5780349We put the input (source) objects into a list:
from <- list(pt0_lonlat=pt0_lonlat, pt0_psmerc=pt0_psmerc, pt0_wgs32n=pt0_wgs32n, pt0_laea89=pt0_laea89)
names(from)## [1] "pt0_lonlat" "pt0_psmerc" "pt0_wgs32n" "pt0_laea89"## pt0_lonlat
## "+proj=longlat +datum=WGS84 +no_defs"
## pt0_psmerc
## "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"
## pt0_wgs32n
## "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs"
## pt0_laea89
## "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"and the target EPSG codes into a vector:
From 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:
## [1] TRUE## [1] FALSErun <- TRUE
if (packageVersion("sp") < "1.3.3") run <- FALSE
if (packageVersion("rgdal") < "1.5.3") run <- FALSE
if (run && !rgdal::new_proj_and_gdal()) run <- FALSESetting it to false gives failing round-trips for two of the four CRS:
out_EPSG_non_viz <- matrix(as.logical(NA), nrow=4, ncol=4)
colnames(out_EPSG_non_viz) <- names(from)
rownames(out_EPSG_non_viz) <- to
for (j in seq(along=from)) {
for (i in seq(along=to)) {
pt1 <- spTransform(from[[j]], CRS(SRS_string = paste0("EPSG:", to[i])))
out_EPSG_non_viz[i, j] <- isTRUE(all.equal(coordinates(from[[i]]),
coordinates(pt1)))
}
}
out_EPSG_non_viz## pt0_lonlat pt0_psmerc pt0_wgs32n pt0_laea89
## 4326 FALSE FALSE FALSE FALSE
## 3857 TRUE TRUE TRUE TRUE
## 32632 TRUE TRUE TRUE TRUE
## 3035 FALSE FALSE FALSE FALSEResetting to the default (TRUE) value gives round-trip
success:
## [1] TRUEout_EPSG <- matrix(as.logical(NA), nrow=4, ncol=4)
colnames(out_EPSG) <- names(from)
rownames(out_EPSG) <- to
for (j in seq(along=from)) {
for (i in seq(along=to)) {
pt1 <- spTransform(from[[j]], CRS(SRS_string = paste0("EPSG:", to[i])))
out_EPSG[i, j] <- isTRUE(all.equal(coordinates(from[[i]]),
coordinates(pt1)))
}
}
out_EPSG## pt0_lonlat pt0_psmerc pt0_wgs32n pt0_laea89
## 4326 TRUE TRUE TRUE TRUE
## 3857 TRUE TRUE TRUE TRUE
## 32632 TRUE TRUE TRUE TRUE
## 3035 TRUE TRUE TRUE TRUEIf 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:
## [1] TRUE## [1] FALSEout_init_non_viz <- matrix(as.logical(NA), nrow=4, ncol=4)
colnames(out_init_non_viz) <- names(from)
rownames(out_init_non_viz) <- to
for (j in seq(along=from)) {
for (i in seq(along=to)) {
pt1 <- spTransform(from[[j]], CRS(paste0("+init=epsg:", to[i])))
out_init_non_viz[i, j] <- isTRUE(all.equal(coordinates(from[[i]]),
coordinates(pt1)))
}
}
out_init_non_viz## pt0_lonlat pt0_psmerc pt0_wgs32n pt0_laea89
## 4326 TRUE TRUE TRUE TRUE
## 3857 TRUE TRUE TRUE TRUE
## 32632 TRUE TRUE TRUE TRUE
## 3035 TRUE TRUE TRUE TRUE## [1] TRUEout_init <- matrix(as.logical(NA), nrow=4, ncol=4)
colnames(out_init) <- names(from)
rownames(out_init) <- to
for (j in seq(along=from)) {
for (i in seq(along=to)) {
pt1 <- spTransform(from[[j]], CRS(paste0("+init=epsg:", to[i])))
out_init[i, j] <- isTRUE(all.equal(coordinates(from[[i]]),
coordinates(pt1)))
}
}
out_init## pt0_lonlat pt0_psmerc pt0_wgs32n pt0_laea89
## 4326 TRUE TRUE TRUE TRUE
## 3857 TRUE TRUE TRUE TRUE
## 32632 TRUE TRUE TRUE TRUE
## 3035 TRUE TRUE TRUE TRUEFinally, 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:
out_EPSG <- matrix(as.logical(NA), nrow=4, ncol=4)
colnames(out_EPSG) <- names(from)
rownames(out_EPSG) <- to
for (j in seq(along=from)) {
for (i in seq(along=to)) {
coo1 <- list_coordOps(comment(slot(from[[j]], "proj4string")),
comment(CRS(SRS_string = paste0("EPSG:", to[i]))),
visualization_order=TRUE)
pt1 <- spTransform(from[[j]], CRS(SRS_string = paste0("EPSG:", to[i])),
coordOp=best_instantiable_coordOp(coo1))
out_EPSG[i, j] <- isTRUE(all.equal(coordinates(from[[i]]),
coordinates(pt1)))
}
}
out_EPSG## pt0_lonlat pt0_psmerc pt0_wgs32n pt0_laea89
## 4326 TRUE TRUE TRUE TRUE
## 3857 TRUE TRUE TRUE TRUE
## 32632 TRUE TRUE TRUE TRUE
## 3035 TRUE TRUE TRUE TRUEThe setting of the list_coordOps()
visualization_order= argument overrides the global option
value (as would setting enforce_xy= in calls to
spTransform()):
out_EPSG_non_viz1 <- matrix(as.logical(NA), nrow=4, ncol=4)
colnames(out_EPSG_non_viz1) <- names(from)
rownames(out_EPSG_non_viz1) <- to
for (j in seq(along=from)) {
for (i in seq(along=to)) {
coo1 <- list_coordOps(comment(slot(from[[j]], "proj4string")),
comment(CRS(SRS_string = paste0("EPSG:", to[i]))),
visualization_order=FALSE)
pt1 <- spTransform(from[[j]], CRS(SRS_string = paste0("EPSG:", to[i])),
coordOp=best_instantiable_coordOp(coo1))
out_EPSG_non_viz1[i, j] <- isTRUE(all.equal(coordinates(from[[i]]),
coordinates(pt1)))
}
}
out_EPSG_non_viz1## pt0_lonlat pt0_psmerc pt0_wgs32n pt0_laea89
## 4326 FALSE FALSE FALSE FALSE
## 3857 TRUE TRUE TRUE TRUE
## 32632 TRUE TRUE TRUE TRUE
## 3035 FALSE FALSE FALSE FALSEout_EPSG_non_viz2 <- matrix(as.logical(NA), nrow=4, ncol=4)
colnames(out_EPSG_non_viz2) <- names(from)
rownames(out_EPSG_non_viz2) <- to
for (j in seq(along=from)) {
for (i in seq(along=to)) {
pt1 <- spTransform(from[[j]], CRS(SRS_string = paste0("EPSG:", to[i])),
enforce_xy=FALSE)
out_EPSG_non_viz2[i, j] <- isTRUE(all.equal(coordinates(from[[i]]),
coordinates(pt1)))
}
}
out_EPSG_non_viz2## pt0_lonlat pt0_psmerc pt0_wgs32n pt0_laea89
## 4326 FALSE FALSE FALSE FALSE
## 3857 TRUE TRUE TRUE TRUE
## 32632 TRUE TRUE TRUE TRUE
## 3035 FALSE FALSE FALSE FALSEproject() for PROJ 6+It 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:
ll <- structure(c(12.1823368669203, 11.9149630062421, 12.3186076188739,
12.6207597184845, 12.9955172054652, 12.6316117692658, 12.4680041846297,
12.4366882666609, NA, NA, -5.78993051516384, -5.03798674888479,
-4.60623015708619, -4.43802336997614, -4.78110320396188, -4.99127125409291,
-5.24836150474498, -5.68430388755925, NA, NA), .Dim = c(10L,
2L), .Dimnames = list(NULL, c("longitude", "latitude")))
try(xy0 <- project(ll, "+proj=moll", legacy=TRUE, verbose=TRUE))## Warning in project(ll, "+proj=moll", legacy = TRUE, verbose = TRUE): PROJ
## support is provided by the sf and terra packages among others## +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step
## +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84if (!exists("xy0")) xy0 <- structure(c(1217100.8468177, 1191302.229156,
1232143.28841193, 1262546.27733232, 1299648.82357849, 1263011.18154638,
1246343.17808186, 1242654.33986052, NA, NA, -715428.207551599,
-622613.577983058, -569301.605757784, -548528.530156422, -590895.949857199,
-616845.926397351, -648585.161643274, -702393.1160979, NA, NA),
.Dim = c(10L, 2L), .Dimnames = list(NULL, c("longitude", "latitude")))
try(ll0 <- project(xy0, "+proj=moll", inv=TRUE, legacy=TRUE, verbose=TRUE))## Warning in project(xy0, "+proj=moll", inv = TRUE, legacy = TRUE, verbose =
## TRUE): PROJ support is provided by the sf and terra packages among others## Warning in project(bb, proj, inv = inv, use_aoi = FALSE): PROJ support is
## provided by the sf and terra packages among others## +proj=pipeline +step +inv +proj=moll +lon_0=0 +x_0=0 +y_0=0
## +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg## [1] TRUEIt is possible to use the coordOp= argument to pass
through a known transformation pipeline:
try(xy1 <- project(ll, "+proj=moll", legacy=TRUE, coordOp=paste("+proj=pipeline +step",
"+proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=moll +lon_0=0 +x_0=0 +y_0=0 ellps=WGS84")))## Warning in project(ll, "+proj=moll", legacy = TRUE, coordOp =
## paste("+proj=pipeline +step", : PROJ support is provided by the sf and terra
## packages among otherstry(ll1 <- project(xy1, "+proj=moll", inv=TRUE, legacy=TRUE, coordOp=paste("+proj=pipeline +step",
"+inv +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg")))## Warning in project(xy1, "+proj=moll", inv = TRUE, legacy = TRUE, coordOp =
## paste("+proj=pipeline +step", : PROJ support is provided by the sf and terra
## packages among others## [1] TRUE## Warning in project(ll, comment(WKT), legacy = TRUE, verbose = TRUE): PROJ
## support is provided by the sf and terra packages among others## +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step
## +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84## Warning in project(xy1, comment(WKT), inv = TRUE, legacy = TRUE, verbose =
## TRUE): PROJ support is provided by the sf and terra packages among others## Warning in project(bb, proj, inv = inv, use_aoi = FALSE): PROJ support is
## provided by the sf and terra packages among others## +proj=pipeline +step +inv +proj=moll +lon_0=0 +x_0=0 +y_0=0
## +ellps=WGS84 +step +proj=unitconvert +xy_in=rad +xy_out=deg## [1] TRUE