Reproducible reprojection in `sf`

Question

I have two machines that give different results for the same CRS transformation. Both are using renv with the same lockfile, and hence they have the same version of sf, namely 1.0-7, so it puzzles me that my results on one machine aren't reproducible on the other.

Here's an example. On system 1:

> sf::st_transform(sf::st_as_sf(data.frame(-82, 29), coords = c(1, 2), crs = 4326), crs = 2163)
Simple feature collection with 1 feature and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 1752402 ymin: -1590928 xmax: 1752402 ymax: -1590928
Projected CRS: US National Atlas Equal Area
                  geometry
1 POINT (1752402 -1590928)

On system 2:

> sf::st_transform(sf::st_as_sf(data.frame(-82, 29), coords = c(1, 2), crs = 4326), crs = 2163)
Simple feature collection with 1 feature and 0 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 1755246 ymin: -1588142 xmax: 1755246 ymax: -1588142
Projected CRS: NAD27 / US National Atlas Equal Area
                  geometry
1 POINT (1755246 -1588142)

I guess what's going on is that one of the depended-upon system-level libraries, such as GDAL, is at a different version, and that's what makes a difference. My questions are:

1. Which of these two different results is right?
2. How do I correct the system giving the wrong result?

Here's sessionInfo() for system 1:

R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.4 LTS
Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C

 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8

 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8

 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C

 [9] LC_ADDRESS=C               LC_TELEPHONE=C

[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods

[7] base
loaded via a namespace (and not attached):
 [1] Rcpp_1.0.9         magrittr_2.0.3     units_0.8-0

 [4] tidyselect_1.1.2   R6_2.5.1           rlang_1.0.6

 [7] fansi_1.0.3        dplyr_1.0.9        tools_4.2.1

[10] grid_4.2.1         KernSmooth_2.23-20 utf8_1.2.2

[13] cli_3.4.1          e1071_1.7-11       DBI_1.1.3

[16] ellipsis_0.3.2     class_7.3-20       tibble_3.1.7

[19] lifecycle_1.0.3    sf_1.0-7           crayon_1.5.1

[22] purrr_0.3.4        vctrs_0.5.0        glue_1.6.2

[25] proxy_0.4-27       compiler_4.2.1     pillar_1.7.0

[28] generics_0.1.3     classInt_0.4-7     renv_0.16.0

[31] pkgconfig_2.0.3

For system 2:

R version 4.2.2 (2022-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.2 LTS
Matrix products: default
BLAS:   /usr/local/lib/R/lib/libRblas.so
LAPACK: /usr/local/lib/R/lib/libRlapack.so
locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C

 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8

 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8

 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C

 [9] LC_ADDRESS=C               LC_TELEPHONE=C

[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods

[7] base
loaded via a namespace (and not attached):
 [1] Rcpp_1.0.9         magrittr_2.0.3     units_0.8-0

 [4] tidyselect_1.1.2   R6_2.5.1           rlang_1.0.6

 [7] fansi_1.0.3        dplyr_1.0.9        tools_4.2.2

[10] grid_4.2.2         KernSmooth_2.23-20 utf8_1.2.2

[13] cli_3.4.1          e1071_1.7-11       DBI_1.1.3

[16] ellipsis_0.3.2     class_7.3-20       tibble_3.1.7

[19] lifecycle_1.0.3    sf_1.0-7           crayon_1.5.1

[22] purrr_0.3.4        vctrs_0.5.0        glue_1.6.2

[25] proxy_0.4-27       compiler_4.2.2     pillar_1.7.0

[28] generics_0.1.3     classInt_0.4-7     renv_0.16.0

[31] pkgconfig_2.0.3

Answer 1

On my computer I'm seeing the same behavior as your system 2. For example, trying to transform a layer in R into EPSG:2163 actually returns a result in EPSG:9311.

These answers, especially the 2nd, IMHO indicate that indeed EPSG:2163 was deprecated and EPSG:9311 is (System 2) is the "right" way to go:

• Choosing projection for Continental United States?
• Deprecated CRS (EPSG:2163) gets reinterpreted as another CRS (EPSG:9311) by GDAL

To get the same (EPSG:9311) behavior on system 1 too, I'd try to update GDAL. (By the way, would be interesting to check the GDAL versions on your system 1 & 2, just to make sure that the version on system 2 is indeed newer.)

---

As for why a CRS gets deprecated, I guess that to avoid confusion and corruption of layers when there are multiple CRS with only a slight difference, serving the same purpose for the same exact area. For example, when someone sends you a vector layer with missing CRS information and says "it's in US National Atlas", which one would you set - EPSG:2163 or EPSG:9311? If you happen to set the wrong one, then the layer would be "corrupt", i.e., the coordinates would not be in the right place and would appear shifted with respect to other layers.

This exact situation exists with the Israeli ITM CRS, with two versions EPSG:2039 and EPSG:6991 (for example see here). Whenever you get a layer in "ITM", with missing CRS, you have to guess, or investigate, which CRS to use, often through trial and error. Sometimes you can tell that the layer is shifted, when there is an obvious reference (such as in a road layer), but sometimes there is no reference and no way to tell. Moreover, there are many layers, even from Israeli governmental sources, where the CRS definition does not match the coordinates. I even encountered layers in EPSG:4326 with turned out to be shifted few tens of meters, due the fact that they were originally in ITM but converted to WGS84 though the wrong "version".