Introduction

Changes in the representation of coordinate reference systems (CRS), and of operations on coordinates, that have been occurring over decades must now be implemented in the way spatial objects are handled in R packages. Up to the 1990s, most spatial data simply used the coordinates given by the local mapping authority; for example, the Meuse bank data set used a planar representation in metres, which turned out to be EPSG:28992, “Amersfoort / RD New”. A major resource for finding out why CRS were specified as they were are Clifford J. Mugnier's columns in Photogrammetric Engineering & Remote Sensing, references to which are available in rgdal; the Netherlands were covered in the February 2003 column:

library(rgdal)
## Loading required package: sp
## rgdal: version: 1.5-12, (SVN revision 1018)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.1.1, released 2020/06/22
## Path to GDAL shared files: /usr/local/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.1.0, August 1st, 2020, [PJ_VERSION: 710]
## Path to PROJ shared files: /home/rsb/.local/share/proj:/usr/local/share/proj:/usr/local/share/proj
## PROJ CDN enabled:FALSE
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
data("GridsDatums")
GridsDatums[grep("Netherlands", GridsDatums$country),]
##                               country      month year ISO
## 50 Aruba and the Netherlands Antilles     (July) 2002 ABW
## 57     The Kingdom of the Netherlands (February) 2003 NLD

While most national mapping agencies defined their own standard geographical and projected CRS, supranational bodies, such as military alliances and colonial administrations often imposed some regularity to facilitate operations across national boundaries. This also led to the creation of the European Petroleum Survey Group (EPSG), because maritime jurisdiction was not orderly, and mattered when countries sharing the same coastal shelf tried to assign conflicting exploration concessions. Experts from oil companies accumulated vast experience, which fed through to the International Standards Organization (ISO, especially TC 211) and the Open Geospatial Consortium (OGC).

Defining the CRS became necessary when integrating other data with a different CRS, and for displaying on a web map background. Many legacy file formats, such as the ESRI Shapefile format, did not mandate the inclusion of the CRS of positional data. Most open source software then used PROJ.4 strings as a flexible representation, but as internationally accepted standards have been reached, in particular ISO 19111, and improved over time by iteration, it is really necessary to change to a modern text representation, known as WKT2 (2019). Now it looks as though almost all corporations and mapping agencies accommodate this representation, and it has been adopted by sp through rgdal, sf and other packages.

demo(meuse, ask=FALSE, package="sp", echo=FALSE)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum Amersfoort in CRS definition
library(mapview)
mapview(meuse, zcol="zinc")

plot of chunk unnamed-chunk-3

Coordinate reference system representation

The mapproj package provided coordinate reference system and projection support for the maps package. From mapproj/src/map.h, line 20, we can see that the eccentricity of the Earth is defined as 0.08227185422, corrresponding to the Clark 1866 ellipsoid (Iliffe and Lott 2008):

ellps <- projInfo("ellps")
(clrk66 <- ellps[ellps$name=="clrk66",])
##      name       major         ell description
## 16 clrk66 a=6378206.4 b=6356583.8 Clarke 1866
#eval(parse(text=clrk66$major))
#eval(parse(text=clrk66$ell))
a <- 6378206.4
b <- 6356583.8
print(sqrt((a^2-b^2)/a^2), digits=10)
## [1] 0.08227185422

With a very few exceptions, projections included in mapproj::mapproject() use the Clarke 1866 ellipsoid, with the remainder using a sphere with the Clarke 1866 major axis radius. The function returns coordinates for visualization in an unknown metric; no inverse projections are available.

Like many other free and open source software projects around 2000, R spatial development chose to use the best open source library and infrastructure then available, PROJ.4. Version 4.4 was published by Frank Warmerdam in 2000, based on Gerald Evenden's earlier work. This earlier work was a library for forward and inverse projection using a key-value string interface to describe the required projection (Evenden 1990). The key-value string is taken as +key=value, where =value could be omitted for some keys, and the definition of each projection is built up from space-separated key-value string, such as +proj=utm +zone=25 +south for Universal Transverse Mercator zone 25 in the southern hemisphere.

PROJ.4 2000-2018

Unlike mapproj, PROJ.4 had begun to introduce the +ellps= key in addition to projection parameters before PROJ.4.4, as some users would not treat Clark 1866 as their natural preference. The need for more care in geodetic specification became pressing from 2000, when civilian use of GPS became as accurate as military use; GPS was typically registered to a more modern geographical CRS than digitized printed maps had been, and most mapping agencies scrambled to update their products and services to “GPS coordinates”.

The ellps= key was followed by the +datum=, +nadgrids= and +towgs84= keys in successive releases to attempt to specify the geodetic model. The +init= key appeared to permit the look-up of sets of key-value strings in the packaged version of a given table with a known authority, typically +init=epsg:<code>, where <code> was the EPSG code of the coordinate reference system. Where +towgs84= was given, a three or seven parameter transformation to the WGS84 datum was provided as a comma-separated string, so that the coordinate reference system also included the inverse coordinate operation from projected coordinates to geographical coordinates defined by the WGS84 datum. This led to the need for a placeholder for geographical coordinates, set as +proj=longlat (or lonlat, or perhaps with reversed axis order latlong or latlon). Some of the issues involved were discussed in a 2010 blog post by Frank Warmerdam.

The PROJ.4 framework functioned well for projection before it was expected to handle datum tranformation too. Within the remit of single mapping agencies, some adaptation could still provide help, say using +nadgrids= in parts of North America (NAD, North American Datum), but where positional data from multiple agencies was being integrated, the framework was showing its age. For example, from about 2010 it was observed that the +datum= and +towgs84= keys sometimes provided contradictory values, leading functions in GDAL reading raster files to prefer +towgs84= values to +datum= values. For users for whom accuracy of better than about 150m was irrelevant, using coordinate reference systems correctly defined in terms of the underlying geographical coordinates was less important, but this is still about five Landsat cells.

While the representation of coordinate reference systems (sometimes supplemented by coordinate operations to transform the underlying geographical coordinates to WGS84) as PROJ.4 strings continued to work adequately, changes were occurring. Many file formats chose to use WKT (well-known text) string representations, starting from the 2007 edition of the ISO standard (ISO 2019). This placed the file reading and writing functions offered by GDAL under stress, especially the exportToProj4() function, as an increasing number of specification components really did not map adequately from the WKT string representation to the PROJ.4 string representation. Another change was that pivoting through a chosen hub when transforming coordinates (in PROJ.4 WGS84) meant that accuracy was lost transforming from the source to the hub, and more was lost from the hub to the target. Why not transform from source to target in one step if possible?

PRØJ

Signalling changes, PROJ.4 changed its name to PRØJ, and began burning through major version numbers. PROJ 5 (2018) introduced transformation pipelines (Knudsen and Evers 2017; Evers and Knudsen 2017), representing coordinate operations using a syntax similar to PROJ.4 strings, but showing the whole operation pipeline. PROJ 6 (2019) followed up by shifting from ad-hoc text files for holding coordinate reference system and coordinate operation metadata to an SQLite database. In an increasing number of cases, more accurate coordinate operations could be supported using open access transformation grids, and the grid files needed were now tabulated in the SQLite database. This database is distributed with PROJ, and kept in a directory on the PROJ search path, usually the only or final directory (getPROJ4libPath() returns the current search path in an attribute).

shpr <- strsplit(attr(getPROJ4libPath(), "search_path"), ifelse(.Platform$OS.type == "unix", ":", ";"))[[1]]
shpr
## [1] "/home/rsb/.local/share/proj" "/usr/local/share/proj"      
## [3] "/usr/local/share/proj"
library(RSQLite)
db <- dbConnect(SQLite(), dbname=file.path(shpr[length(shpr)], "proj.db"))
dbListTables(db)
##  [1] "alias_name"                               
##  [2] "area"                                     
##  [3] "authority_list"                           
##  [4] "authority_to_authority_preference"        
##  [5] "axis"                                     
##  [6] "celestial_body"                           
##  [7] "compound_crs"                             
##  [8] "concatenated_operation"                   
##  [9] "concatenated_operation_step"              
## [10] "conversion"                               
## [11] "conversion_method"                        
## [12] "conversion_param"                         
## [13] "conversion_table"                         
## [14] "coordinate_operation_method"              
## [15] "coordinate_operation_view"                
## [16] "coordinate_operation_with_conversion_view"
## [17] "coordinate_system"                        
## [18] "crs_view"                                 
## [19] "deprecation"                              
## [20] "ellipsoid"                                
## [21] "geodetic_crs"                             
## [22] "geodetic_datum"                           
## [23] "geoid_model"                              
## [24] "grid_alternatives"                        
## [25] "grid_packages"                            
## [26] "grid_transformation"                      
## [27] "helmert_transformation"                   
## [28] "helmert_transformation_table"             
## [29] "metadata"                                 
## [30] "object_view"                              
## [31] "other_transformation"                     
## [32] "prime_meridian"                           
## [33] "projected_crs"                            
## [34] "supersession"                             
## [35] "unit_of_measure"                          
## [36] "vertical_crs"                             
## [37] "vertical_datum"
dbReadTable(db, "metadata")
##            key                                                    value
## 1 EPSG.VERSION                                                  v9.8.12
## 2    EPSG.DATE                                               2020-06-19
## 3 ESRI.VERSION                                            ArcMap 10.8.1
## 4    ESRI.DATE                                               2020-05-24
## 5  IGNF.SOURCE https://geodesie.ign.fr/contenu/fichiers/IGNF.v3.1.0.xml
## 6 IGNF.VERSION                                                    3.1.0
## 7    IGNF.DATE                                               2019-05-24

PROJ 7 (2020) reconfigured the transformation grids, now using the Geodetic TIFF Grid (GTG) format, and created pathways for on-demand download (typically using a content download netwoek (CDN) over CURL) of chunks of such grids to a local user-writable cache held in another SQLite database. After little change from the late 1990s to early 2018, PROJ.4 has incremented its major version number three times in three years, and by 2021 (PROJ 8), the pre-existing application programming interface will be history. In addition, GDAL 3 (2019) has tightened its links with PROJ >= 6, and exportToProj4() now says:

“Use of this function is discouraged. Its behavior 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/api/ogrspatialref.html).

For these reasons, sf and sp are changing from PROJ.4 strings representing coordinate reference systems to WKT2:2019 strings, as described in this r-spatial blog. Most users who had been relying on file reading to set the coordinate reference systems of objects will not notice much difference, and legacy PROJ.4 strings can still be used to create new, authority-free definitions if need be.

Coordinate operations

The introduction of interactive mapping using mapview and tmap among other packages highlights the need to set the coordinate reference system (CRS) of objects correctly, so that zooming does not reveal embarrassing divergences. Using the location of the Broad Street pump disabled by Dr John Snow to stop a cholera epidemic in Soho, London, in 1854 (Brody et al. 2000), we can start to see what steps are being taken. The point location of the pump is given in projected coordinates, defined in the British National Grid. The workflow used by mapview::mapview() is to transform first to WGS84 (EPSG:4326) using sf::st_transform() if need be, before permitting leaflet to project to Web Mercator (EPSG:3857) internally.

b_pump <- readOGR(system.file("vectors/b_pump.gpkg", package="rgdal"))
## 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.999601 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
## OGR data source with driver: GPKG 
## Source: "/tmp/RtmpFom9DX/Rinst3b3a16978ffef/rgdal/vectors/b_pump.gpkg", layer: "b_pump"
## with 1 features
## It has 1 fields
## Integer64 fields read as strings:  cat

Reading the file loses the PROJ.4 +datum= key-value pair, but the WKT2:2019 string is complete:

proj4string(b_pump)
## Warning in proj4string(b_pump): CRS object has comment, which is lost in output
## [1] "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs"
if (packageVersion("sp") > "1.4.1") {
  WKT <- wkt(b_pump)
} else {
  WKT <- comment(slot(b_pump, "proj4string"))
}
cat(WKT, "\n")
## PROJCRS["Transverse_Mercator",
##     BASEGEOGCRS["GCS_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["unnamed",
##         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.999601,
##             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["easting",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["northing",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

Pipelines

PROJ.4 assumed that the from/source and to/target coordinate reference system definitions involved in coordinate operations each contained the specifications necessary to get from source to WGS84 and then on from WGS84 to target. PR$\phi$J drops this assumption, searching among many candidate coordinate operations for viable pipelines. The search is conducted using the tables given in the proj.db SQLite database, which is now backed by authorities, and regularly updated at each release to the current upstream state. The tables are searched to find lists of candidates.

cov <- dbReadTable(db, "coordinate_operation_view")
cov[grep("OSGB", cov$name), c(1, 3, 4, 9, 16)]
##                  table_name   code                           name
## 162     grid_transformation   5338        OSGB 1936 to ETRS89 (1)
## 163     grid_transformation   5339        OSGB 1936 to WGS 84 (7)
## 213     grid_transformation   7709        OSGB 1936 to ETRS89 (2)
## 214     grid_transformation   7710        OSGB 1936 to WGS 84 (9)
## 517     grid_transformation 108057  ETRS_1989_To_OSGB_1936_OSTN15
## 518     grid_transformation 108058   WGS_1984_To_OSGB_1936_OSTN15
## 519     grid_transformation 108059      OSGB_1936_To_Xrail84_NTv2
## 720  helmert_transformation   1195        OSGB 1936 to WGS 84 (1)
## 721  helmert_transformation   1196        OSGB 1936 to WGS 84 (2)
## 722  helmert_transformation   1197        OSGB 1936 to WGS 84 (3)
## 723  helmert_transformation   1198        OSGB 1936 to WGS 84 (4)
## 724  helmert_transformation   1199        OSGB 1936 to WGS 84 (5)
## 822  helmert_transformation   1314        OSGB 1936 to WGS 84 (6)
## 823  helmert_transformation   1315      OSGB 1936 to ED50 (UKOOA)
## 1377 helmert_transformation   5622        OSGB 1936 to WGS 84 (8)
## 2094 helmert_transformation 108089 OSGB_1936_To_WGS_1984_8_BAD_DX
## 2249 helmert_transformation 108336 OSGB_1936_To_WGS_1984_NGA_7PAR
##                                         method_name accuracy
## 162                                            NTv2     0.03
## 163                                            NTv2     1.00
## 213                                            NTv2     0.03
## 214                                            NTv2     1.00
## 517                                            NTv2     0.03
## 518                                            NTv2     1.00
## 519                                            NTv2     0.50
## 720         Geocentric translations (geog2D domain)    21.00
## 721         Geocentric translations (geog2D domain)    10.00
## 722         Geocentric translations (geog2D domain)    21.00
## 723         Geocentric translations (geog2D domain)    18.00
## 724         Geocentric translations (geog2D domain)    35.00
## 822  Position Vector transformation (geog2D domain)     2.00
## 823  Position Vector transformation (geog2D domain)     2.00
## 1377        Geocentric translations (geog2D domain)     3.00
## 2094        Geocentric translations (geog2D domain)     5.00
## 2249      Coordinate Frame rotation (geog2D domain)    21.00

The same search can be conducted directly without using RSQLite to query the database tables, searching by source and target CRS, and in the near future also by area of interest. If we search using only the degraded PROJ.4 string, we only find a ballpark accuracy coordinate operation, yielding a pipeline with two steps, inverse projection to geographical coordinates in radians, and conversion from radians to degrees. Note that rgdal::spTransform() and its wrapper sp::spTransform() use PROJ for coordinate operations:

list_coordOps(paste0(proj4string(b_pump), " +type=crs"), "EPSG:4326")
## Warning in proj4string(b_pump): CRS object has comment, which is lost in output
## Candidate coordinate operations found:  1 
## Strict containment:  FALSE 
## Visualization order:  TRUE 
## Source: +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +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.999601 +x_0=400000 +y_0=-100000 +ellps=airy
##              +step +proj=unitconvert +xy_in=rad +xy_out=deg

The description component “+ axis order change (2D)” refers to the EPSG:4326 definition, which specifies Northings/Latitude as the first axis, and Eastings/Longitude as the second axis; in sp/rgdal workflows, it is assumed that GIS/visualization order with Eastings/Longitude as the first 2D axis and Northings/Latitude as the second axis is preferred. Because the input data are in GIS/visualization order already, the steps to swap axes to standards conformity and then back to GIS/visualization order cancel each other out.

Setting the internal control option set_transform_wkt_comment() to FALSE, we use only the degraded PROJ.4 string when transforming. spTransform() undertakes the same search, chooses the best instantiable coordinate operation on its first pass, then uses that pipeline on all objects. The pipeline specification of the coordinate operation may be retrieved using get_last_coordOp():

set_transform_wkt_comment(FALSE)
isballpark <- spTransform(b_pump, CRS(SRS_string="EPSG:4326"))
get_last_coordOp()
## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +x_0=400000 +y_0=-100000 +ellps=airy +step +proj=unitconvert +xy_in=rad +xy_out=deg"

The coordinate returned is unfortunately in Ingestre Place, not Broad Street.

print(coordinates(isballpark), digits=10)
##          coords.x1  coords.x2
## [1,] -0.1351070464 51.5127926

Let us repeat the search using the WKT2 string; here we see that providing a well-specified CRS representation allows us to choose 2m accuracy for the coordinate operation. Further, we can also see that, had we had access to a named grid, we could have achieved 1m accuracy:

list_coordOps(WKT, "EPSG:4326")
## Candidate coordinate operations found:  8 
## Strict containment:  FALSE 
## Visualization order:  TRUE 
## Source: PROJCRS["Transverse_Mercator", BASEGEOGCRS["GCS_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["unnamed", 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.999601, 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["easting", east, ORDER[1],
##         LENGTHUNIT["metre", 1, ID["EPSG", 9001]]],
##         AXIS["northing", north, ORDER[2], LENGTHUNIT["metre",
##         1, ID["EPSG", 9001]]]]
## Target: EPSG:4326 
## Best instantiable operation has accuracy: 2 m
## Description: Inverse of unnamed + 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.999601 +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 7 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

The Helmert transformation has parameters retrieved from the PROJ SQLite database (code 1314):

helm <- dbReadTable(db, "helmert_transformation_table")
helm[helm$code == "1314",c(1:3, 15:17, 20:22, 25)]
##     auth_name code                    name      tx       ty     tz   rx    ry
## 239      EPSG 1314 OSGB 1936 to WGS 84 (6) 446.448 -125.157 542.06 0.15 0.247
##        rz scale_difference
## 239 0.842          -20.489
dbDisconnect(db)

Using the WKT2 CRS representation, we can achieve 2m accuracy (or as the table says in other fields: “Oil exploration. Accuracy better than 4m and generally better than 2m”:

set_transform_wkt_comment(TRUE)
is2m <- spTransform(b_pump, CRS(SRS_string="EPSG:4326"))
get_last_coordOp()
## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +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"

The output point is close to the Broad Street pump:

print(coordinates(is2m), digits=10)
##          coords.x1   coords.x2
## [1,] -0.1367127374 51.51330218

It is over 100m West-North-West of the Ingestre Place position:

c(spDists(isballpark, is2m)*1000)
## [1] 125.0577
## Checking rgeos availability: TRUE
c(maptools::gzAzimuth(coordinates(isballpark), coordinates(is2m)))
## coords.x1 
## -62.98022

This was about as good as one could get prior to PROJ 7 without downloading the missing grid file manually, and installing the downloaded file in a directory that would usually not be user-writable. The whole set of grids can be downloaded and installed manually for workgroups needing to be sure that the same grids are available to all users, as has been the case in the past as well.

The rgdal::project() uses the underlying geographical coordinate reference system, and does not transform, so using the degraded PROJ.4 string and WKT2 give the same output, and because the input is projected, we take the inverse:

(a <- project(coordinates(b_pump), proj4string(b_pump), inv=TRUE, verbose=TRUE))
## Warning in proj4string(b_pump): CRS object has comment, which is lost in output
## proj=pipeline step inv proj=tmerc lat_0=49 lon_0=-2 k=0.999601
## x_0=400000 y_0=-100000 ellps=airy step proj=unitconvert xy_in=rad
## xy_out=deg
##      coords.x1 coords.x2
## [1,] -0.135107  51.51279
(b <- project(coordinates(b_pump), WKT, inv=TRUE))
##      coords.x1 coords.x2
## [1,] -0.135107  51.51279

The projected points only inverse project the projected coordinates using the specified projection

all.equal(a, b)
## [1] TRUE
c(spDists(coordinates(isballpark), a)*1000)
## [1] 0

Transformation grid files

PROJ 7 introduced on-demand downloading of (chunks of) transformation grids from a content delivery network to a user-writable directory on the PROJ search path (usually the first path component). The status of the downloaded grids is stored in another SQLite database, cache.db. Let is unlink this file, and check that rgdal is running with on-demand downloading disabled:

run
## [1] TRUE
shpr
## [1] "/home/rsb/.local/share/proj" "/usr/local/share/proj"      
## [3] "/usr/local/share/proj"
unlink(file.path(shpr[1], "cache.db"))
rgdal:::is_proj_network_enabled()
## [1] FALSE

We can (for now, this will be handled later in a function) enable on-demand download setting an environment variable; we see that with this setting, network download is enabled:

Sys.setenv("PROJ_NETWORK"="ON")
rgdal:::is_proj_network_enabled()
## [1] FALSE

When we then search for candidate coordinate operations, we see that the operation using an absent grid now sees that download is enabled, and proposes the 1m accuracy candidate, because the required grid can be downloaded:

list_coordOps(WKT, "EPSG:4326")
## Candidate coordinate operations found:  8 
## Strict containment:  FALSE 
## Visualization order:  TRUE 
## Source: PROJCRS["Transverse_Mercator", BASEGEOGCRS["GCS_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["unnamed", 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.999601, 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["easting", east, ORDER[1],
##         LENGTHUNIT["metre", 1, ID["EPSG", 9001]]],
##         AXIS["northing", north, ORDER[2], LENGTHUNIT["metre",
##         1, ID["EPSG", 9001]]]]
## Target: EPSG:4326 
## Best instantiable operation has accuracy: 2 m
## Description: Inverse of unnamed + 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.999601 +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 7 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

On making the transformation, we may see that the coordinate operation takes longer than expected, because on first pass the grid is downloaded from the network:

system.time(is1m <- spTransform(b_pump, CRS(SRS_string="EPSG:4326")))
##    user  system elapsed 
##    0.03    0.00    0.03

The coordinate operation used now specifies the grid in the pipeline:

get_last_coordOp()
## [1] "+proj=pipeline +step +inv +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.999601 +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"

The coordinate values differ little from the 2m accuracy Helmert pipeline:

print(coordinates(is1m), digits=10)
##          coords.x1   coords.x2
## [1,] -0.1367127374 51.51330218

as we can see, the 1m accuracy point is 1.7m from the 2m accuracy point, just to the West:

c(spDists(is2m, is1m)*1000)
## [1] 0
c(maptools::gzAzimuth(coordinates(is1m), coordinates(is2m)))
## coords.x1 
##        NA

If we look in the SQLite database of downloaded grids, we see that the grid components that were downloaded. Here we have not yet used the area of interest to limit the number of chunks involved:

library(RSQLite)
db <- dbConnect(SQLite(), dbname=file.path(shpr[1], "cache.db"))
(tbls <- dbListTables(db))
## character(0)
if ("chunks" %in% tbls) dbReadTable(db, "chunks")
dbDisconnect(db)

Finally, we disable grid download to return to the status existing when rgdal was attached.

Sys.setenv("PROJ_NETWORK"="OFF")
rgdal:::is_proj_network_enabled()
## [1] FALSE

The outcome positions are shown here; at zoom 18 we can see that the 1m accuracy green point matches the Open Street Map location of the pump very well:

library(mapview)
mapview(is2m, map.type="OpenStreetMap", legend=FALSE) + mapview(is1m, col.regions="green", legend=FALSE) + mapview(isballpark, col.regions="red", legend=FALSE)

plot of chunk unnamed-chunk-41

References

Brody, H., M. R. Rip, P. Vinten-Johansen, N. Paneth, and S. Rachman. 2000. “Map-Making and Myth-Making in Broad Street: The London Cholera Epidemic, 1854.” Lancet 356: 64–68.

Evenden, Gerald I. 1990. Cartographic Projection Procedures for the UNIX Environment — a User’s Manual. (http://download.osgeo.org/proj/OF90-284.pdf).

Evers, Kristian, and Thomas Knudsen. 2017. Transformation Pipelines for Proj.4. (https://www.fig.net/resources/proceedings/fig_proceedings/fig2017/papers/iss6b/ISS6B_evers_knudsen_9156.pdf).

Iliffe, Jonathan, and Roger Lott. 2008. Datums and Map Projections: For Remote Sensing, GIS and Surveying. Boca Raton: CRC.

ISO. 2019. ISO 19111:2019 Geographic Information – Referencing by Coordinates. (https://www.iso.org/standard/74039.html).

Knudsen, Thomas, and Kristian Evers. 2017. Transformation Pipelines for Proj.4. (https://meetingorganizer.copernicus.org/EGU2017/EGU2017-8050.pdf).

Warmerdam, Frank. 2010. Slaying the Datum Shift Dragon. (http://fwarmerdam.blogspot.com/2010/03/in-last-few-weeks-i-believe-i-have-made.html).