2

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
2
  • System 2 seems to use 9311 instead. Looks related - gis.stackexchange.com/questions/377099/…
    – margusl
    Commented Mar 22, 2023 at 21:36
  • @margusl Yep, looks like that's the culprit. Setting OSR_USE_NON_DEPRECATED=NO gets me system 1's behavior. I'll leave this question open in the hope that someone can answer "Which of these two different results is right?". In particular, I don't understand why a CRS was deprecated, nor which CRS is the true US National Atlas projection. Commented Mar 23, 2023 at 12:42

1 Answer 1

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:

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".

1
  • 1
    If you register for EPSG, you can look at deprecated objects. 2163 says "No longer supported by EPSG because datum information is required for unambiguous spatial referencing." USGS used an authalic sphere based on Clarke 1866 so that was the GeoCRS of 2163. EPSG then defined one based on NAD27 as a replacement. They won't give the same results and wgs84 to nad27 may be using a transformation as well.
    – mkennedy
    Commented Apr 2, 2023 at 14:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.