I am trying to convert coordinates in EPSG:4283 to EPSG:4326 in sf in R, on 2 computers with different dependency versions.
Computer A has GDAL 3.8.2, GEOS 3.11.2 and PROJ 9.31, Windows 10
Computer B has GDAL 3.5.2, GEOS 3.9.3 and PROJ 8.2.1, Windows 8
The default st_transform output for computer A does not alter the input geometry:
> library(sf)
> in_data <- st_as_sf(data.frame(x=151.262387590000003, y=-33.912053239999998766), coords=c('x', 'y'), crs='EPSG:4283')
> st_transform(in_data, 'EPSG:4326')$geometry[[1]]
POINT (151.262387590000003 -33.912053239999998766)
Whereas on computer B there is a difference after 5dp (which is significant for this task):
> library(sf)
> in_data <- st_as_sf(data.frame(x=151.262387590000003, y=-33.912053239999998766), coords=c('x', 'y'), crs='EPSG:4283')
> st_transform(in_data, 'EPSG:4326')$geometry[[1]]
POINT (151.26239300639633711 -33.912040561252297266)
To try and figure out what was causing this, I first printed out the available pipelines for both computers with sf_proj_pipelines('EPSG:4283', 'EPSG:4326') which gave the same 2 operations for both (albeit in different orders, with computer A having the no operation pipeline first):
[1] "+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=0.06155 +y=-0.01087 +z=-0.04019 +rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 +s=-0.009994 +convention=coordinate_frame +step +inv +proj=cart +ellps=WGS84 +step +proj=pop +v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg"
[2] "proj=noop +ellps=GRS80"
This gives me the options, but not what was specifically used by st_transform (worth noting that both of those options had the same "accuracy" value on both computers). Using rgdal's list_coordOps I tried to find the 'best' candidate pipeline. The output for computer A was:
Candidate coordinate operations found: 3
Strict containment: FALSE
Visualization order: TRUE
Source: EPSG:4283
Target: EPSG:4326
Best instantiable operation has accuracy: 3 m
Description: axis order change (2D) + GDA94 to WGS 84 (1) + axis order change (2D)
Definition: +proj=noop +ellps=GRS80
Operation 3 is lacking 1 grid with accuracy 3 m
Missing grid: au_icsm_GDA94_GDA2020_conformal_and_distortion.tif
URL: https://cdn.proj.org/au_icsm_GDA94_GDA2020_conformal_and_distortion.tif
And computer B:
Candidate coordinate operations found: 3
Strict containment: FALSE
Visualization order: TRUE
Source: EPSG:4283
Target: EPSG:4326
Best instantiable operation has accuracy: 3 m
Description: axis order change (2D) + GDA94 to WGS 84 (2) + axis order change (2D)
Definition: +proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=push +v_3 +step +proj=cart +ellps=GRS80 +step +proj=helmert +x=0.06155 +y=-0.01087
+z=-0.04019 +rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 +s=-0.009994 +convention=coordinate_frame +step +inv +proj=cart +ellps=WGS84 +step +proj=pop
+v_3 +step +proj=unitconvert +xy_in=rad +xy_out=deg
Operation 3 is lacking 1 grid with accuracy 3 m
Missing grid: au_icsm_GDA94_GDA2020_conformal_and_distortion.tif
URL: https://cdn.proj.org/au_icsm_GDA94_GDA2020_conformal_and_distortion.tif
When I try and enter the Definition of the best candidate from computer B on an st_transform on computer A:
Warning message:
In st_transform.sfc(st_geometry(x), crs, ...) :
pipeline not found in PROJ-suggested candidate transformations
And even when entering it on computer B:
hgridshift: could not find required grid(s).
pipeline: Pipeline: Bad step definition: proj=hgridshift (File not found or invalid)
hgridshift: could not find required grid(s).
pipeline: Pipeline: Bad step definition: proj=hgridshift (File not found or invalid)
Warning message:
In st_transform.sfc(st_geometry(x), crs, ...) :
pipeline not found in PROJ-suggested candidate transformations
Both of which result in empty geometries.
In order to get output on computer B to match that of computer A, I can use a no operation projstring:
+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=noop +ellps=GRS80
But I can't figure out how to get computer A to match the output of computer B. I suspect computer B is missing the grid required (I've tried sf_proj_network(TRUE) to get grids or other proj info online but this doesn't change anything) to make the shift it puts forward as the best candidate, so it falls back on a different transformation that I cannot find the definition for. But I'm not sure.
What pipeline can I use, or other method, to get results on computer A to match that of computer B, without changing the underlying PROJ or GDAL versions?
