CDN grid st_transform() problems with Windows and macOS

[rsbivand]

See the thread starting https://stat.ethz.ch/pipermail/r-sig-geo/2021-October/028796.html, originally: https://stackoverflow.com/questions/69507752/how-to-apply-a-grid-transformation-in-sf-with-proj7

Uploaded an R file

cdn_211010_sf_rgdal.R.zip

# at the very beginning of the R session, define the user-writable
# directory
td <- tempfile()
dir.create(td)
Sys.setenv("PROJ_USER_WRITABLE_DIRECTORY"=td)   
library(rgdal)
rgdal_extSoftVersion()
# check the PROJ search paths and check that the user-writable directory
# is empty and that the CDN is not enabled
(pths <- get_proj_search_paths())
list.files(pths[1])
is_proj_CDN_enabled()
# define the objects and check that the search paths are still OK
library(sf)
get_proj_search_paths()
pt_27572 <- data.frame(X=55824.4970,Y=2394454.2120)
pt_2154 <- data.frame(X=107242.8310,Y=6832277.1820)
ptSf_27572 <- st_as_sf(pt_27572, coords = c("X", "Y"), crs=27572)
ptSf_2154 <- st_as_sf(pt_2154, coords = c("X", "Y"), crs=2154)
# enable the CDN and that the user-writable directory is still empty
enable_proj_CDN()
is_proj_CDN_enabled()
list.files(pths[1])
# list candidate coordinate operations
(o <- list_coordOps("EPSG:27572", "EPSG:2154"))
# transform using the first returned coordinate operation pipeline
ptSf_2154_grid <- st_as_sf(spTransform(as(ptSf_27572, "Spatial"),
  CRS("EPSG:2154"), coordOp=o[1, "definition"]))
# print the coordinate operation used
get_last_coordOp()
# check that the user-writable directory is populated
list.files(pths[1])
# check distance to defined point in target crs
st_distance(ptSf_2154_grid, ptSf_2154)

three Rout files from running this script on Linux, win10 and macOS

cdn_211010_sf_rgdal_Rout.zip

and a CSV file from running variants. It turned out that setting the pipeline was essential also on Linux, for unknown reasons, but that setting the pipeline didn't work for sf because st_transform was only looking at the bundled PROJ metadata directory.

cdn_211010_table.csv

[edzer]

I tried this, which made me believe that on Ubuntu, only setting the pipeline gets the CDN grid:

library(sf)
# Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
sf_proj_search_paths()
# [1] "/home/edzer/.local/share/proj" "/usr/share/proj"              
sf_proj_network()
# [1] FALSE

pt_27572 <- data.frame(X=55824.4970,Y=2394454.2120) # set coordinates of one point in EPSG:27572

pt_2154 <- data.frame(X=107242.8310,Y=6832277.1820) # known accurate

ptSf_27572 <- st_as_sf(pt_27572, coords = c("X", "Y"), crs = 27572)

ptSf_2154 <- st_as_sf(pt_2154, coords = c("X", "Y"), crs = 2154)

sf_proj_network(TRUE) # allow search for online datum grids in the
# [1] "https://cdn.proj.org"

# [1] "https://cdn.proj.org"

(p = sf_proj_pipelines("EPSG:27572", "EPSG:2154"))
# Candidate coordinate operations found:  3 
# Strict containment:     FALSE 
# Axis order auth compl:  FALSE 
# Source:  EPSG:27572 
# Target:  EPSG:2154 
# Best instantiable operation has accuracy: 1 m
# Description: Inverse of Lambert zone II + NTF (Paris) to NTF (1) + Inverse
#              of RGF93 to NTF (2) + Lambert-93
# Definition:  +proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8
#              +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000
#              +ellps=clrk80ign +pm=paris +step +proj=hgridshift
#              +grids=fr_ign_ntf_r93.tif +step +proj=lcc
#              +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44
#              +x_0=700000 +y_0=6600000 +ellps=GRS80
# OP: (fr_ign_gr3df97a.tif) seem to be found for accurate transformation from EPSG:27572 to EPSG:2154
# here we get a slightly different grid name

ptSf_2154_grid <- st_transform(ptSf_27572, crs=2154) # apply transformation

st_distance(ptSf_2154_grid, ptSf_2154) # incorrect (ungridded)
# Units: [m]
#         [,1]
# [1,] 73.7975
# transformation, the distance should be zero. 3.777 m is the known error
# for the ungridded transformation.

# apply transformation using pipeline:
(ptSf_pi = st_transform(ptSf_27572, "EPSG:2154", pipeline = p[1,]$definition))
# Simple feature collection with 1 feature and 0 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 107242.8 ymin: 6832277 xmax: 107242.8 ymax: 6832277
# Projected CRS: RGF93 / Lambert-93
#                   geometry
# 1 POINT (107242.8 6832277)

st_distance(ptSf_pi, ptSf_2154)
# Units: [m]
#              [,1]
# [1,] 0.0004808082

[florisvdh]

I've been struggling with similar things some time ago, using sf 1.0-2 on Linux Mint 20 (Ubuntu 20.04) with PROJ 7.2.1 and GDAL 3.2.1 (the latter is now updated to 3.3.2 on ubuntugis-unstable PPA). The sorting of the sf_proj_pipelines dataframe depended on the availability of grids in place, which makes sense, but having preset sf_proj_network(TRUE) did not trigger grid downloads when the grids were missing (tried several times in new R sessions). Only after manual grid downloading, sf_proj_pipelines() saw the grids, but the resulting dataframe was not always re-sorted accordingly if sf_proj_pipelines() had already been called before the download.

Behaviour seemed a bit unstable/different between iterations. I didn't yet find the time to delve deeper and try to isolate and reproduce problems, but since these observations seem related to this issue, it seemed better to inform you already.

[rsbivand]

Extracted example:

library(sf)
(sp <- sf_proj_search_paths())
sf_proj_network()
pt_27572 <- data.frame(X=55824.4970,Y=2394454.2120)
pt_2154 <- data.frame(X=107242.8310,Y=6832277.1820)
ptSf_27572 <- st_as_sf(pt_27572, coords = c("X", "Y"), crs = 27572)
ptSf_2154 <- st_as_sf(pt_2154, coords = c("X", "Y"), crs = 2154)
sf_proj_network(TRUE)
(p = sf_proj_pipelines("EPSG:27572", "EPSG:2154"))
ptSf_2154_grid <- st_transform(ptSf_27572, crs=2154)
st_distance(ptSf_2154_grid, ptSf_2154) 
(ptSf_pi = st_transform(ptSf_27572, "EPSG:2154", pipeline = p[1,]$definition))
st_distance(ptSf_pi, ptSf_2154)

Yes, this works on Linuxen, but still has the problem of the least accurate pipeline being chosen in st_transform() when no pipeline is given. This problem (for these crs) is the same for rgdal::spTransform(). Without setting a temporary user-writable directory, when CDN is OFF, any grids in the default user-writable directory should be used, but do not seem to be accessed:

library(RSQLite)
db <- dbConnect(SQLite(), dbname=file.path(sp[1], "cache.db"))
dbListTables(db)
dbReadTable(db, "chunks")
dbDisconnect(db)

Maybe fr_ign_ntf_r93.tif is used rather than fr_ign_gr3df97a.tif because it is smaller, less download?

I've also established that the 74m distance is omitting any datum transform, but:

ptSf_2154_via_CRS84 <- st_transform(st_transform(ptSf_27572, crs="OGC:CRS84"), crs=2154)
st_distance(ptSf_2154_via_CRS84, ptSf_2154)
# Units: [m]
#          [,1]
# [1,] 3.777444

where 27572 -> CRS84 is a 3-coefficient Helmert transformation. So st_transform() via GDAL seems to fall back to using the deprecated WGS84 hub when grids are not working on win10 and macOS.

[rsbivand]

win10:

> library(sf)
Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
> (sp <- sf_proj_search_paths())
[1] "C:/Users/RB/Documents/R/win-library/4.1/sf/proj"
> sf_proj_network()
[1] FALSE
> pt_27572 <- data.frame(X=55824.4970,Y=2394454.2120)
> pt_2154 <- data.frame(X=107242.8310,Y=6832277.1820)
> ptSf_27572 <- st_as_sf(pt_27572, coords = c("X", "Y"), crs = 27572)
> ptSf_2154 <- st_as_sf(pt_2154, coords = c("X", "Y"), crs = 2154)
> sf_proj_network(TRUE)
[1] "https://cdn.proj.org"
> (p = sf_proj_pipelines("EPSG:27572", "EPSG:2154"))
Candidate coordinate operations found:  3
Strict containment:     FALSE
Axis order auth compl:  FALSE
Source:  EPSG:27572
Target:  EPSG:2154
Best instantiable operation has accuracy: 1 m
Description: Inverse of Lambert zone II + NTF (Paris) to NTF (1) + NTF to
             RGF93 (1) + Lambert-93
Definition:  +proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8
             +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000
             +ellps=clrk80ign +pm=paris +step +proj=push +v_3
             +step +proj=cart +ellps=clrk80ign +step
             +proj=xyzgridshift +grids=fr_ign_gr3df97a.tif
             +grid_ref=output_crs +ellps=GRS80 +step +inv
             +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step
             +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44
             +x_0=700000 +y_0=6600000 +ellps=GRS80
> ptSf_2154_grid <- st_transform(ptSf_27572, crs=2154)
> st_distance(ptSf_2154_grid, ptSf_2154)
Units: [m]
         [,1]
[1,] 3.777346
> (ptSf_pi = st_transform(ptSf_27572, "EPSG:2154", pipeline = p[1,]$definition$
Error in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
  OGRCreateCoordinateTransformation(): transformation not available
In addition: Warning messages:
1: In CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
  GDAL Error 1: PROJ: proj_create: Error -38 (failed to load datum shift file): Pipeline: Bad step definition: proj=xyzgridshift (failed to load datum shift file)
2: In CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
  GDAL Error 6: Cannot instantiate pipeline +proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris +step +proj=push +v_3 +step +proj=cart +ellps=clrk80ign +step +proj=xyzgridshift +grids=fr_ign_gr3df97a.tif +grid_ref=output_crs +ellps=GRS80 +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80
> st_distance(ptSf_pi, ptSf_2154)
Error in st_crs(x) : object 'ptSf_pi' not found
>

[rsbivand]

macOS:

> library(sf)
Linking to GEOS 3.9.1, GDAL 3.2.3, PROJ 7.2.1
> (sp <- sf_proj_search_paths())
[1] "/Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/library/sf/proj"
> sf_proj_network()
[1] FALSE
> pt_27572 <- data.frame(X=55824.4970,Y=2394454.2120)
> pt_2154 <- data.frame(X=107242.8310,Y=6832277.1820)
> ptSf_27572 <- st_as_sf(pt_27572, coords = c("X", "Y"), crs = 27572)
> ptSf_2154 <- st_as_sf(pt_2154, coords = c("X", "Y"), crs = 2154)
> sf_proj_network(TRUE)
[1] "https://cdn.proj.org"
> (p = sf_proj_pipelines("EPSG:27572", "EPSG:2154"))
Candidate coordinate operations found:  3 
Strict containment:     FALSE 
Axis order auth compl:  FALSE 
Source:  EPSG:27572 
Target:  EPSG:2154 
Best instantiable operation has accuracy: 1 m
Description: Inverse of Lambert zone II + NTF (Paris) to NTF (1) + NTF to
             RGF93 (1) + Lambert-93
Definition:  +proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8
             +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000
             +ellps=clrk80ign +pm=paris +step +proj=push +v_3
             +step +proj=cart +ellps=clrk80ign +step
             +proj=xyzgridshift +grids=fr_ign_gr3df97a.tif
             +grid_ref=output_crs +ellps=GRS80 +step +inv
             +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step
             +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44
             +x_0=700000 +y_0=6600000 +ellps=GRS80
> ptSf_2154_grid <- st_transform(ptSf_27572, crs=2154)
> st_distance(ptSf_2154_grid, ptSf_2154) 
Units: [m]
         [,1]
[1,] 3.777346
> (ptSf_pi = st_transform(ptSf_27572, "EPSG:2154", pipeline = p[1,]$definition)) 
Error in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  : 
  OGRCreateCoordinateTransformation() returned NULL: PROJ available?
In addition: Warning messages:
1: In CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
  GDAL Error 1: PROJ: proj_create: Error -38 (failed to load datum shift file): Pipeline: Bad step definition: proj=xyzgridshift (failed to load datum shift file)
2: In CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
  GDAL Error 6: Cannot instantiate pipeline +proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris +step +proj=push +v_3 +step +proj=cart +ellps=clrk80ign +step +proj=xyzgridshift +grids=fr_ign_gr3df97a.tif +grid_ref=output_crs +ellps=GRS80 +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80
> st_distance(ptSf_pi, ptSf_2154)
Error in st_crs(x) : object 'ptSf_pi' not found

[florisvdh]

In Linux I get the same results as you (sf 1.0.3, linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1) with just one grid downloaded and available (and yes, auto-download worked!):

> dbReadTable(db, "chunks")
  id                                      url offset data_id data_size
1  1 https://cdn.proj.org/fr_ign_gr3df97a.tif      0       1     16384

It might be the grid naming mismatch that causes the ballpark transformation to be used - indeed it is (running below code after running your script):

> (ptSf_pi = st_transform(ptSf_27572, "EPSG:2154", pipeline = p[3,]$definition))
Simple feature collection with 1 feature and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 107316.6 ymin: 6832280 xmax: 107316.6 ymax: 6832280
Projected CRS: RGF93 / Lambert-93
                  geometry
1 POINT (107316.6 6832280)
> st_distance(ptSf_pi, ptSf_2154)
Units: [m]
        [,1]
[1,] 73.7975
> ptSf_pi
Simple feature collection with 1 feature and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 107316.6 ymin: 6832280 xmax: 107316.6 ymax: 6832280
Projected CRS: RGF93 / Lambert-93
                  geometry
1 POINT (107316.6 6832280)

This may be a bug in PROJ 7.2.1, since with its cs2cs tool one gets the same ballpark transformation:

$ echo 55824.4970 2394454.2120 | cs2cs EPSG:27572 EPSG:2154
107316.59   6832279.59 0.00

Moreover PROJ does not complain about a missing grid either, for the first pipeline, so it should have used the first pipeline it seems:

$ projinfo -s EPSG:27572 -t EPSG:2154 --summary
Candidate operations found: 3
unknown id, Inverse of Lambert zone II + NTF (Paris) to NTF (1) + Inverse of RGF93 to NTF (2) + Lambert-93, 1 m, France - onshore - mainland and Corsica.
unknown id, Inverse of Lambert zone II + NTF (Paris) to NTF (1) + NTF to RGF93 (1) + Lambert-93, 1 m, France - onshore - mainland and Corsica., at least one grid missing
unknown id, Inverse of Lambert zone II + Inverse of Transformation from RGF93 to NTF (Paris) altered to use prime meridian of RGF93 + Ballpark geographic offset from NTF (Paris) altered to use prime meridian of RGF93 to RGF93 + Lambert-93, unknown accuracy, World, has ballpark transformation

While writing a tutorial, it was observed that default st_transform() will trigger location dependent choice of pipeline, i.e. choosing a ballpark transformation for coordinates that fall outside the CRS's area of usage (geographic bounding box), and choosing the more accurate pipeline only for coordinates that fall within. (This makes sense but I think it is not documented to that detail in https://keen-swartz-3146c4.netlify.app/sf.html#pipelines.) However I checked the coordinates of the currently investigated point and such thing doesn't seem to be the case here - so excluding this as a possible cause.

Ideally this is re-tested with latest PROJ and if it's still there, seems like a bug to be solved at PROJ?

[rsbivand]

No, I have:

> sf_extSoftVersion()
          GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
 "3.10.0beta3"        "3.3.2"        "8.1.1"         "true"         "true" 
          PROJ 
       "8.1.1" 

And pipeline listing is with strict containment FALSE, so various odd things are going on.

[florisvdh]

The second pipeline is the one with fr_ign_gr3df97a.tif, which is said to be 'not found on the system':

$ projinfo -s EPSG:27572 -t EPSG:2154 -o PROJ
Candidate operations found: 3
-------------------------------------
Operation No. 1:

unknown id, Inverse of Lambert zone II + NTF (Paris) to NTF (1) + Inverse of RGF93 to NTF (2) + Lambert-93, 1 m, France - onshore - mainland and Corsica.

PROJ string:
+proj=pipeline
  +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742
        +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris
  +step +proj=hgridshift +grids=fr_ign_ntf_r93.tif
  +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000
        +y_0=6600000 +ellps=GRS80

-------------------------------------
Operation No. 2:

unknown id, Inverse of Lambert zone II + NTF (Paris) to NTF (1) + NTF to RGF93 (1) + Lambert-93, 1 m, France - onshore - mainland and Corsica., at least one grid missing

PROJ string:
+proj=pipeline
  +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742
        +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris
  +step +proj=push +v_3
  +step +proj=cart +ellps=clrk80ign
  +step +proj=xyzgridshift +grids=fr_ign_gr3df97a.tif +grid_ref=output_crs
        +ellps=GRS80
  +step +inv +proj=cart +ellps=GRS80
  +step +proj=pop +v_3
  +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000
        +y_0=6600000 +ellps=GRS80

Grid fr_ign_gr3df97a.tif needed but not found on the system. Can be obtained at https://cdn.proj.org/fr_ign_gr3df97a.tif

-------------------------------------
Operation No. 3:

unknown id, Inverse of Lambert zone II + Inverse of Transformation from RGF93 to NTF (Paris) altered to use prime meridian of RGF93 + Ballpark geographic offset from NTF (Paris) altered to use prime meridian of RGF93 to RGF93 + Lambert-93, unknown accuracy, World, has ballpark transformation

PROJ string:
+proj=pipeline
  +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742
        +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris
  +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000
        +y_0=6600000 +ellps=GRS80

While:

$ projinfo --searchpaths
/home/floris/.local/share/proj
/usr/share/proj

So PROJ seems to say it has r_ign_ntf_r93.tif (operation 1) - where is it? - and says it doesn't have fr_ign_gr3df97a.tif (operation 2) - but it's in the cache.db.

sf_proj_pipelines() has triggered download of fr_ign_gr3df97a.tif but not fr_ign_ntf_r93.tif. PROJ seems to think the first pipeline is available, the second is not and it performs the third (ballpark).

After deleting cache.db and in a new R session, with default sf_proj_network() as FALSE:

> library(sf)
Linking to GEOS 3.9.1, GDAL 3.3.2, PROJ 7.2.1
> (sp <- sf_proj_search_paths())
[1] "/home/floris/.local/share/proj" "/usr/share/proj"               
> sf_proj_network()
[1] FALSE
> pt_27572 <- data.frame(X=55824.4970,Y=2394454.2120)
> pt_2154 <- data.frame(X=107242.8310,Y=6832277.1820)
> ptSf_27572 <- st_as_sf(pt_27572, coords = c("X", "Y"), crs = 27572)
> ptSf_2154 <- st_as_sf(pt_2154, coords = c("X", "Y"), crs = 2154)
> (p = sf_proj_pipelines("EPSG:27572", "EPSG:2154"))
Candidate coordinate operations found:  3 
Strict containment:     FALSE 
Axis order auth compl:  FALSE 
Source:  EPSG:27572 
Target:  EPSG:2154 
Best instantiable operation has accuracy: 1 m
Description: Inverse of Lambert zone II + NTF (Paris) to NTF (1) + Inverse of RGF93 to NTF
             (2) + Lambert-93
Definition:  +proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0
             +k_0=0.99987742 +x_0=600000 +y_0=2200000 +ellps=clrk80ign
             +pm=paris +step +proj=hgridshift +grids=fr_ign_ntf_r93.tif +step
             +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000
             +y_0=6600000 +ellps=GRS80
Operation 2 is lacking 1 grid with accuracy 1 m
Missing grid: fr_ign_gr3df97a.tif 
URL: https://cdn.proj.org/fr_ign_gr3df97a.tif 
> tibble::as_tibble(p)
# A tibble: 3 × 9
  id         description definition has_inverse accuracy axis_order grid_count instantiable containment
  <chr>      <chr>       <chr>      <lgl>          <dbl> <lgl>           <int> <lgl>        <lgl>      
1 "pipeline" Inverse of… +proj=pip… TRUE               1 FALSE               1 TRUE         FALSE      
2 ""         Inverse of… +          FALSE              1 FALSE               1 FALSE        FALSE      
3 "pipeline" Inverse of… +proj=pip… TRUE              NA FALSE               0 TRUE         FALSE 

Now the second pipeline is not instantiable (odd/empty values appear for id and definition columns), but first and third remain the same.

And the second one has a reference to fr_ign_ntf_r93.tif:

> p[2,]
Candidate coordinate operations found:  1 
Strict containment:     NA 
Axis order auth compl:  NA 
Source:  EPSG:27572 
Target:  EPSG:2154 
Best instantiable operation has only ballpark accuracy 
Description: NA 
Definition:  NA 
Operation 1 is lacking 1 grid with accuracy 1 m
Missing grid: fr_ign_ntf_r93.tif 
Name: /usr/share/proj/ntf_r93.gsb 
URL: https://cdn.proj.org/fr_ign_ntf_r93.tif

Compare both above chunks: just depending on the way of outputting (p or p[2,]) another missing grid is referred for operation 2.

Probably some things are being switched here (by sf?), maybe due to equal accuracy between the first two operations? But at least PROJ seems to do some things wrong (selecting ballpark transformation - see cs2cs example before - and apparently wrong assumptions about missing and non-missing grids). Maybe a mistake in proj.db.

[rsbivand]

Same here: the local cache holds (part of) fr_ign_gr3df97a.tif, but it is not found; note also that cs2cs seems to use a grid when --no-ballpark is stipulated:

$ projinfo -s EPSG:27572 -t EPSG:2154 -o PROJ
Candidate operations found: 3
-------------------------------------
Operation No. 1:

unknown id, Inverse of Lambert zone II + NTF (Paris) to NTF (1) + Inverse of RGF93 to NTF (2) + Lambert-93, 1 m, France - onshore - mainland and Corsica.

PROJ string:
+proj=pipeline
  +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742
        +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris
  +step +proj=hgridshift +grids=fr_ign_ntf_r93.tif
  +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000
        +y_0=6600000 +ellps=GRS80

-------------------------------------
Operation No. 2:

unknown id, Inverse of Lambert zone II + NTF (Paris) to NTF (1) + NTF to RGF93 (1) + Lambert-93, 1 m, France - onshore - mainland and Corsica., at least one grid missing

PROJ string:
+proj=pipeline
  +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742
        +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris
  +step +proj=push +v_3
  +step +proj=cart +ellps=clrk80ign
  +step +proj=xyzgridshift +grids=fr_ign_gr3df97a.tif +grid_ref=output_crs
        +ellps=GRS80
  +step +inv +proj=cart +ellps=GRS80
  +step +proj=pop +v_3
  +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000
        +y_0=6600000 +ellps=GRS80

Grid fr_ign_gr3df97a.tif needed but not found on the system. Can be obtained at https://cdn.proj.org/fr_ign_gr3df97a.tif

-------------------------------------
Operation No. 3:

unknown id, Inverse of Lambert zone II + Inverse of Transformation from RGF93 to NTF (Paris) altered to use prime meridian of RGF93 + Ballpark geographic offset from NTF (Paris) altered to use prime meridian of RGF93 to RGF93 + Lambert-93, unknown accuracy, World, has ballpark transformation

PROJ string:
+proj=pipeline
  +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742
        +x_0=600000 +y_0=2200000 +ellps=clrk80ign +pm=paris
  +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000
        +y_0=6600000 +ellps=GRS80
$ projinfo --searchpaths
/home/rsb/.local/share/proj
/usr/local/share/proj
/usr/local/share/proj
$ ls /home/rsb/.local/share/proj
cache.db
$ sqlite3 /home/rsb/.local/share/proj/cache.db
SQLite version 3.34.1 2021-01-20 14:10:07
Enter ".help" for usage hints.
sqlite> select * from chunks;
1|https://cdn.proj.org/fr_ign_ntf_r93.tif|0|1|16384
2|https://cdn.proj.org/fr_ign_gr3df97a.tif|0|2|16384
3|https://cdn.proj.org/fr_ign_gr3df97a.tif|32768|3|16384
4|https://cdn.proj.org/fr_ign_gr3df97a.tif|65536|4|16384
sqlite> 
$ echo 55824.4970 2394454.2120 | cs2cs "EPSG:27572" "EPSG:2154"
107316.59   6832279.59 0.00
$ echo 55824.4970 2394454.2120 | cs2cs --no-ballpark "EPSG:27572" "EPSG:2154"
107242.83   6832277.18 0.00

[rsbivand]

And after st_write(ptSf_27572, "ptSf_27572.gpkg"); ogr2ogr chooses the wrong pipeline too: ptSf_27572.gpkg.zip

$ ogr2ogr -t_srs EPSG:2154 ptSf_2154.gpkg ptSf_27572.gpkg
$ ogrinfo -al ptSf_2154.gpkg
INFO: Open of `ptSf_2154.gpkg'
      using driver `GPKG' successful.

Layer name: ptSf_27572
Geometry: Point
Feature Count: 1
Extent: (107316.589068, 6832279.594280) - (107316.589068, 6832279.594280)
Layer SRS WKT:
PROJCRS["RGF93 / Lambert-93",
    BASEGEOGCRS["RGF93",
        DATUM["Reseau Geodesique Francais 1993",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4171]],
    CONVERSION["Lambert-93",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",46.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",3,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",49,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",44,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",700000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",6600000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["France - onshore and offshore, mainland and Corsica."],
        BBOX[41.15,-9.86,51.56,10.38]],
    ID["EPSG",2154]]
Data axis to CRS axis mapping: 1,2
FID Column = fid
Geometry Column = geom
OGRFeature(ptSf_27572):1
  POINT (107316.589067701 6832279.59427963)

[rsbivand]

@rouault - any ideas? Is the 27572 to 2154 transformation particularly difficult to resolve in 7.2.1 and 8.1.1? The previous comment is just PROJ and GDAL command line applications, but falls into the same holes as found when running PROJ and GDAL from R packages. Why might the ballpark accuracy be chosen before 1m accuracy when cache.db contains the required grids? Is a PROJ issue justified?

[rouault]

I'm a bit confused by all this, but a few things to have in mind:

• grids chunks cached in share/proj/cache.db are only taken into account if PROJ_NETWORK=ON. A grid there is ignored otherwise.
• the order of operations in projinfo output might depend on which grids are available in /usr/share/proj (or other searched directories), if PROJ_NETWORK=ON is set and the value of the --grid-check switch (known-available has the same effect as enabling PROJ_NETWORK=ON or having all the grids in /usr/share/proj)
• with PROJ master on Linux, echo 55824.4970 2394454.2120 | PROJ_LIB=../../PROJ-data/fr_ign:data PROJ_DEBUG=2 src/cs2cs EPSG:27572 EPSG:2154 outputs 107242.83 6832277.18 using the "NTF to RGF93 v1 (1)" transformation with the fr_ign_gr3df97a.tif grid . Same if setting PROJ_NETWORK=ON

[rsbivand]

Thanks @rouault ! With my built from source 8.1.1, I have:

$ echo 55824.4970 2394454.2120 | PROJ_NETWORK=ON PROJ_DEBUG=2 cs2cs "EPSG:27572" "EPSG:2154"
pj_open_lib(proj.db): call fopen(/usr/local/share/proj/proj.db) - succeeded
pj_open_lib(proj.ini): call fopen(/usr/local/share/proj/proj.ini) - succeeded
pj_open_lib(fr_ign_gr3df97a.tif): call fopen(fr_ign_gr3df97a.tif) - failed
Using https://cdn.proj.org/fr_ign_gr3df97a.tif
pj_open_lib(fr_ign_gr3df97a.tif): call fopen(fr_ign_gr3df97a.tif) - failed
Using https://cdn.proj.org/fr_ign_gr3df97a.tif
Using coordinate operation Inverse of Lambert zone II + Inverse of Transformation from RGF93 to NTF (Paris) altered to use prime meridian of RGF93 + Ballpark geographic offset from NTF (Paris) altered to use prime meridian of RGF93 to RGF93 + Lambert-93
107316.59   6832279.59 0.00

The confusion is real, I'm afraid. In the issue there are two problems, or even three - the problem of cached grids not being used, the problem of an apparent fallback to a WGS84 hub when using GDAL for transformation, and deployment problems when we build Windows and macOS static packages, meaning that we haven't defined the user-writable directory properly.

[rouault]

I can actually reproduce with 8.1.1 the issue with local grids not being used with that particular transformation, and confirm it works with PROJ master due to the fix of https://github.com/OSGeo/PROJ/pull/2884 that affects transformations when a CRS has a non-Greenwhich prime meridian. If you use a more eatern input coordinate, like 255824.4970 2394454.2120, you don't run into that bug.

[rsbivand]

With a source build of the PROJ master:

$ echo 55824.4970 2394454.2120 | PROJ_LIB=data PROJ_NETWORK=ON PROJ_DEBUG=2 bin/cs2cs EPSG:27572 EPSG:2154
pj_open_lib(proj.db): call fopen(data/proj.db) - succeeded
pj_open_lib(proj.ini): call fopen(data/proj.ini) - succeeded
pj_open_lib(fr_ign_gr3df97a.tif): call fopen(data/fr_ign_gr3df97a.tif) - failed
Using https://cdn.proj.org/fr_ign_gr3df97a.tif
pj_open_lib(fr_ign_gr3df97a.tif): call fopen(data/fr_ign_gr3df97a.tif) - failed
Using https://cdn.proj.org/fr_ign_gr3df97a.tif
Using coordinate operation Inverse of Lambert zone II + NTF (Paris) to NTF (1) + NTF to RGF93 v1 (1) + Lambert-93
107242.83   6832277.18 0.00

[rsbivand]

@edzer Further example (Fedora 34, PROJ 8.1.1, GDAL 3.3.2):

library(sf)
b_pump <- st_read(system.file("gpkg/b_pump.gpkg", package="sf"))
sf_proj_network(TRUE)
sf_proj_search_paths()
bp <- sf_proj_pipelines("EPSG:27700", "EPSG:32630")
bp$accuracy
o0 <- st_transform(b_pump, "EPSG:32630")
o1 <- st_transform(b_pump, "EPSG:32630", pipeline=bp[1, "definition"])
o7 <- st_transform(b_pump, "EPSG:32630", pipeline=bp[7, "definition"])
# o8 <- st_transform(b_pump, "EPSG:32630", pipeline=bp[8, "definition"]) # ERROR
o9 <- st_transform(b_pump, "EPSG:32630", pipeline=bp[9, "definition"])
st_distance(rbind(o0, o1, o7, o9))
print(st_coordinates(st_geometry(rbind(o0, o1, o7, o9))), digits=12)

Does the CDN work for you? No cache.db is created in the default location, where is it? For comparison:

library(rgdal)
b_pump <- readOGR(system.file("gpkg/b_pump.gpkg", package="sf"))
enable_proj_CDN()
bp <- list_coordOps("EPSG:27700", "EPSG:32630")
bp$accuracy
o0 <- spTransform(b_pump, CRS("EPSG:32630"))
get_last_coordOp()
o1 <- spTransform(b_pump, CRS("EPSG:32630"), coordOp=bp[1, "definition"])
o7 <- spTransform(b_pump, CRS("EPSG:32630"), coordOp=bp[7, "definition"])
o8 <- spTransform(b_pump, CRS("EPSG:32630"), coordOp=bp[8, "definition"])
o9 <- spTransform(b_pump, CRS("EPSG:32630"), coordOp=bp[9, "definition"])
spDists(rbind(o0, o1, o7, o8, o9))
print(coordinates(rbind(o0, o1, o7, o8, o9)), digits=12)

Here, proj.db created in default location, difference 0.08411507 from different +k= values (0.999601 in o0, 0.9996012717 in o8), both use uk_os_OSTN15_NTv2_OSGBtoETRS.tif, uk_os_OSTN15_NTv2_OSGBtoETRS.tif.

[rsbivand]

Is this: https://gdal.org/api/ogr_srs_api.html#_CPPv420OSRGetPROJAuxDbPathsv relevant?

[rouault]

No cache.db is created in the default location, where is it?

https://proj.org/usage/network.html?highlight=network#caching + https://proj.org/resource_files.html#user-writable-directory might help

Is this: https://gdal.org/api/ogr_srs_api.html#_CPPv420OSRGetPROJAuxDbPathsv relevant?

probably not. This is about storing extra CRS and coordinate operations in auxiliary databases. Not about grid caching

[rsbivand]

Thanks, I was monitoring /home/rsb/.local/share/proj after removing the existing cache.db. For the rgdal example, which only uses PROJ and does not touch OSR for transformation, a new cache.db was created in /home/rsb/.local/share/proj. However, sf uses OSR not PROJ for transformation (apart from direct projinfo-style pipeline lookup); nothing appeared in /home/rsb/.local/share/proj. For my Fedora 34 system, no XDG_DATA_HOME is set.

Trying OSRGetPROJAuxDbPaths(), nothing there:

> library(sf)
Linking to GEOS 3.10.0beta3, GDAL 3.3.2, PROJ 8.1.1
> sf:::CPL_get_proj_search_paths(character(0))
[1] "/home/rsb/.local/share/proj" "/usr/local/share/proj"      
[3] "/usr/local/share/proj"      
> sf:::CPL_get_proj_auxdb_paths()
character(0)

[rouault]

However, sf uses OSR not PROJ for transformation (apart from direct projinfo-style pipeline lookup)

PROJ networking must be enabled through the PROJ_NETWORK=ON environment variable (and not GDAL config option) to be enabled in OSR usage.

[rsbivand]

sf uses PROJ to enable the network: https://github.com/r-spatial/sf/blob/ea94efe83cf2ef575bf6e9d59e04a07301d06b28/src/proj.cpp#L193-L207 The equivalent rgdal code is in https://r-forge.r-project.org/scm/viewvc.php/pkg/src/proj_network7.cpp?view=markup&root=rgdal, and returns the output of proj_context_set_enable_network(PJ_DEFAULT_CTX, TRUE);.

[rsbivand]

And adding PROJ_NETWORK=ON brings sf and rgdal into alignment:

Sys.setenv("PROJ_NETWORK"="ON")
library(sf)
b_pump <- st_read(system.file("gpkg/b_pump.gpkg", package="sf"))
sf_proj_network()
sf_proj_search_paths()
bp <- sf_proj_pipelines("EPSG:27700", "EPSG:32630")
bp$accuracy
o0 <- st_transform(b_pump, "EPSG:32630")
o1 <- st_transform(b_pump, "EPSG:32630", pipeline=bp[1, "definition"])
o7 <- st_transform(b_pump, "EPSG:32630", pipeline=bp[7, "definition"])
o8 <- st_transform(b_pump, "EPSG:32630", pipeline=bp[8, "definition"])
o9 <- st_transform(b_pump, "EPSG:32630", pipeline=bp[9, "definition"])
st_distance(rbind(o0, o1, o7, o8, o9))
print(st_coordinates(st_geometry(rbind(o0, o1, o7, o8, o9))), digits=12)

0 is the default transformation, 1 and 7 are Helmert transformations, 8 is a grid transformation and 9 a ballpark transformation. With PROJ_NETWORK=ON projinfo --spatial-test intersects -s EPSG:27700 -t EPSG:32630 -o PROJ, the corresponding order is: 2, 8, 1, 9. The difference between default and the best from projinfo is +k=0.9996012717 in the specified pipeline and +k=0.999601 in the default (I think).

The same works on macOS, PROJ 7.2.1. Without setting PROJ_NETWORK=ON but with sf_proj_network(TRUE), o8 fails, and no cache.db appears in ${HOME}/Library/Application Support/proj; with PROJ_NETWORK=ON, it appears (same as Linux).

And the same works with winlibs R 4.1.1; same symptoms of no cache.db when not setting PROJ_NETWORK=ON but with sf_proj_network(TRUE).

Larger problems with UCRT3 R-4.2.0-devel; cache.db is created but not populated, so has no chunks. This is both for the CRAN Windows UCRT binary and for a local source install. The same problem exists for rgdal, so maybe the static linking of curl or sqlite3 or tiff isn't OK

[rsbivand]

A simple check function for the +pm-grid https://github.com/OSGeo/PROJ/pull/2884 bug in PROJfor sf:

sf_check_pm_grid_bug = function(x, to_crs) {
    p = sf_proj_pipelines(st_crs(x)$wkt, to_crs)
    res = rep(FALSE, nrow(p))
    pv = package_version(sf:::CPL_proj_version())
    if (pv >= "7.2.1" && pv <= "8.1.1") {
        res = grepl("\\+pm=", p$definition) & !grepl("greenwich", p$definition) & grepl("grid", p$definition)
    }
    res
}

In the context of the example at https://github.com/r-spatial/sf/issues/1815#issuecomment-940776012:

p = sf_proj_pipelines("EPSG:27572", "EPSG:2154")
o = sf_check_pm_grid_bug(ptSf_27572, "EPSG:2154")
p$accuracy[!o]

shows the residual accuracy from found pipelines, suggesting use of a manual hub avoiding the bug.

For rgdal:

rgdal_check_pm_grid_bug <- function(x, to_crs) {
    p <- list_coordOps(wkt(x), to_crs)
    res <- rep(FALSE, nrow(p))
    pv <- attr(getPROJ4VersionInfo(), "short")
    if (pv >= 721 && pv <= 811) {
        res <- grepl("\\+pm=", p$definition) & !grepl("greenwich", p$definition) & grepl("grid", p$definition)
    }
    res
}

[rsbivand]

Windows UCRT3 problem with PROJ network; share/proj/* and bin/cs2cs.exe copied from MXE cross-compilation to a separate directory, and run from the rtools42 ucrt64 console:

$ echo 529393.5 181020.6 | PROJ_LIB=./proj PROJ_NETWORK=ON PROJ_DEBUG=2 ./cs2cs.exe EPSG:27700 EPSG:32630pj_open_lib(proj.db): call fopen(./proj\proj.db) - succeeded
pj_open_lib(proj.ini): call fopen(./proj\proj.ini) - succeeded
Using coordinate operation Inverse of British National Grid + OSGB 1936 to WGS 84 (9) + UTM zone 30N
pj_open_lib(uk_os_OSTN15_NTv2_OSGBtoETRS.tif): call fopen(./proj\uk_os_OSTN15_NTv2_OSGBtoETRS.tif) - failed
pj_open_lib(OSTN15_NTv2_OSGBtoETRS.gsb): call fopen(./proj\OSTN15_NTv2_OSGBtoETRS.gsb) - failed
Cannot open https://cdn.proj.org/uk_os_OSTN15_NTv2_OSGBtoETRS.tif: SSL certificate problem: unable to get local issuer certificate
*       * inf

So the current MXE PROJ build is missing something as suspected.

Also as run from /C/rtools42/x86_64-w64-mingw32.static.posix with no need for copying:

 UCRT64 /C/rtools42/x86_64-w64-mingw32.static.posix
$ echo 529393.5 181020.6 | PROJ_LIB=share/proj PROJ_NETWORK=ON PROJ_DEBUG=2 bin/cs2cs.exe EPSG:27700 EPSG:32630
pj_open_lib(proj.db): call fopen(share/proj\proj.db) - succeeded
pj_open_lib(proj.ini): call fopen(share/proj\proj.ini) - succeeded
Using coordinate operation Inverse of British National Grid + OSGB 1936 to WGS 84 (9) + UTM zone 30N
pj_open_lib(uk_os_OSTN15_NTv2_OSGBtoETRS.tif): call fopen(share/proj\uk_os_OSTN15_NTv2_OSGBtoETRS.tif) - failed
pj_open_lib(OSTN15_NTv2_OSGBtoETRS.gsb): call fopen(share/proj\OSTN15_NTv2_OSGBtoETRS.gsb) - failed
Cannot open https://cdn.proj.org/uk_os_OSTN15_NTv2_OSGBtoETRS.tif: SSL certificate problem: unable to get local issuer certificate
*       * inf

[edzer]

So the current MXE PROJ build is missing something as suspected.

PROJ, or libcurl?

As of the sf issue, I guess that sf_proj_network(TRUE) should then always call Sys.setenv("PROJ_NETWORK"="ON") to inform GDAL, right?

[rouault]

PROJ, or libcurl?

At least some years ago, it was quite tricky to have curl to find the CA certificates on Windows installation. They are a number of manual options that were added in GDAL curl interface to control that, but we didn't push that to PROJ curl interface. Perhaps they are now ways to build curl on Windows with a hardcoded path to where the CA certificates are.

[rsbivand]

The PROJ, GEOS and GDAL static binaries are cross-compiled using MXE on Linux, and copied to Windows. So if the path was hardcoded during cross-compilation, it will not correspond on the target system, most likely. I'm in touch with the R-core member handling migration to UCRT using MXE hopefully from R 4.2.0 due April 2022, so I hope we'll be prepared for the UCRT roadbump before we see problems among users. I see that OSGeo4W typically grabs all the MS CRT versions (that ever existed?), UCRT simplified a lot (effectively moving to UTF-8 like everyone else) https://developer.r-project.org/Blog/public/2020/05/02/utf-8-support-on-windows/index.html.

[rsbivand]

@rouault:

Following analysis by Tomas Kalibera, it seems the MXE cross-compilation build of curl did provide the certificate bundle, but its location got lost in transfer to target; cache.db gets only 48KB and is empty:

$ echo 55824.4970 2394454.2120 | PROJ_NETWORK=ON PROJ_DEBUG=2 PROJ_LIB=proj ./cs2cs "EPSG:27572" "EPSG:2154"
pj_open_lib(proj.db): call fopen(proj\proj.db) - succeeded
pj_open_lib(proj.ini): call fopen(proj\proj.ini) - succeeded
pj_open_lib(fr_ign_gr3df97a.tif): call fopen(proj\fr_ign_gr3df97a.tif) - failed
Cannot open https://cdn.proj.org/fr_ign_gr3df97a.tif: SSL certificate problem: unable to get local issuer certificate
xyzgridshift: could not find required grid(s).
Pipeline: Bad step definition: proj=xyzgridshift (failed to load datum shift file)
pj_open_lib(fr_ign_gr3df97a.tif): call fopen(proj\fr_ign_gr3df97a.tif) - failed
Cannot open https://cdn.proj.org/fr_ign_gr3df97a.tif: SSL certificate problem: unable to get local issuer certificate
xyzgridshift: could not find required grid(s).
Pipeline: Bad step definition: proj=xyzgridshift (failed to load datum shift file)
Using coordinate operation Inverse of Lambert zone II + Inverse of Transformation from RGF93 to NTF (Paris) altered to use prime meridian of RGF93 + Ballpark geographic offset from NTF (Paris) altered to use prime meridian of RGF93 to RGF93 + Lambert-93
107316.59       6832279.59 0.00

Resolution - to set the required environment variable, and cache.db goes to 64KB and seems viable:

$ echo 55824.4970 2394454.2120 | PROJ_NETWORK=ON PROJ_DEBUG=2 PROJ_LIB=proj CURL_CA_BUNDLE="C:/Program Files/R/R-devel/etc/curl-ca-bundle.crt" ./cs2cs "EPSG:27572" "EPSG:2154"pj_open_lib(proj.db): call fopen(proj\proj.db) - succeeded
pj_open_lib(proj.ini): call fopen(proj\proj.ini) - succeeded
pj_open_lib(fr_ign_gr3df97a.tif): call fopen(proj\fr_ign_gr3df97a.tif) - failed
Using https://cdn.proj.org/fr_ign_gr3df97a.tif
pj_open_lib(fr_ign_gr3df97a.tif): call fopen(proj\fr_ign_gr3df97a.tif) - failed
Using https://cdn.proj.org/fr_ign_gr3df97a.tif
Using coordinate operation Inverse of Lambert zone II + Inverse of Transformation from RGF93 to NTF (Paris) altered to use prime meridian of RGF93 + Ballpark geographic offset from NTF (Paris) altered to use prime meridian of RGF93 to RGF93 + Lambert-93
107316.59       6832279.59 0.00

So the problem was certificate bundling.

[rsbivand]

And:

$ echo 529393.5 181020.6 | PROJ_LIB=share/proj PROJ_NETWORK=ON PROJ_DEBUG=2 CURL_CA_BUNDLE="C:/Program Files/R/R-devel/etc/curl-ca-bundle.crt" bin/cs2cs.exe EPSG:27700 EPSG:32630
pj_open_lib(proj.db): call fopen(share/proj\proj.db) - succeeded
pj_open_lib(proj.ini): call fopen(share/proj\proj.ini) - succeeded
Using coordinate operation Inverse of British National Grid + OSGB 1936 to WGS 84 (9) + UTM zone 30N
pj_open_lib(uk_os_OSTN15_NTv2_OSGBtoETRS.tif): call fopen(share/proj\uk_os_OSTN15_NTv2_OSGBtoETRS.tif) - failed
pj_open_lib(OSTN15_NTv2_OSGBtoETRS.gsb): call fopen(share/proj\OSTN15_NTv2_OSGBtoETRS.gsb) - failed
Using https://cdn.proj.org/uk_os_OSTN15_NTv2_OSGBtoETRS.tif
698673.86       5710795.91 0.00

Resolution until R-devel UCRT is updated is to add: Sys.setenv("CURL_CA_BUNDLE"="C:/Program Files/R/R-devel/etc/curl-ca-bundle.crt") manually.

[rsbivand]

Today's build of R-devel for UCRT has fixed the problem for packages using the CDN.