Missing params in WKT string generated by st_crs()

tcwilkinson commented 2 years ago

Describe the bug Essential parameter of the (recently updated) peirce_q projection do not seem to appear in the WKT form after st_crs(). This has obvious consequences for translation etc.

This may well be an upstream issue, but can't get the same WKT using proj cmd line tools so not sure how st_crs() generates the WKT string (presumably via C api?).

Any help tracking this down, much appreciated.

To Reproduce

Will need the current master branch version of PROJ (>8.2.1) in order to test with "+shape"

> st_crs("+proj=peirce_q +lon_0=25")
Coordinate Reference System:
  User input: +proj=peirce_q +lon_0=25 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["PROJ peirce_q"],
        PARAMETER["lon_0",25,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

OK that's fine, but add in the parameter +shape=square [in HEAD branch of proj, in 8.2.1 the value was (wrongly) +type=square], and you end up with exactly the same, except the user input:-

> st_crs("+proj=peirce_q +lon_0=25 +shape=square")
Coordinate Reference System:
  User input: +proj=peirce_q +lon_0=25 +shape=square 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["PROJ peirce_q"],
        PARAMETER["lon_0",25,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

For comparison, PROJ via command line projinfo provides:

$ projinfo -o WKT2:2019 "+proj=peirce_q +lon0=25"
WKT2:2019 string:
CONVERSION["PROJ-based coordinate operation",
    METHOD["PROJ-based operation method: +proj=peirce_q +lon0=25"]]

$ projinfo -o WKT2:2019 "+proj=peirce_q +lon0=25 +shape=square"
WKT2:2019 string:
CONVERSION["PROJ-based coordinate operation",
    METHOD["PROJ-based operation method: +proj=peirce_q +lon0=25 +shape=square"]]

Additional context The context is trying to get the new peirce_q projection to work for a mapping project, hence am currently using the master branch (HEAD) build of PROJ. (Have been contributing to PROJ library to get this to work, and trying to find + iron out any issues before next PROJ release.)

```r > sessionInfo() R version 4.1.2 (2021-11-01) Platform: x86_64-apple-darwin17.0 (64-bit) Running under: macOS Monterey 12.1 Matrix products: default LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib locale: [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] ggnewscale_0.4.5 ggrepel_0.9.1 smoothr_0.2.2 hillshader_0.1.0 [5] stars_0.5-6 abind_1.4-5 terra_1.5-12 raster_3.5-11 [9] sp_1.4-6 lwgeom_0.2-8 sf_1.0-6 forcats_0.5.1 [13] stringr_1.4.0 dplyr_1.0.7 purrr_0.3.4 readr_2.1.1 [17] tidyr_1.1.4 tibble_3.1.6 ggplot2_3.3.5 tidyverse_1.3.1 loaded via a namespace (and not attached): [1] httr_1.4.2 jsonlite_1.7.3 foreach_1.5.1 [4] modelr_0.1.8 assertthat_0.2.1 cellranger_1.1.0 [7] progress_1.2.2 pillar_1.6.4 backports_1.4.1 [10] lattice_0.20-45 glue_1.6.0 digest_0.6.29 [13] rvest_1.0.2 colorspace_2.0-2 htmltools_0.5.2 [16] pkgconfig_2.0.3 broom_0.7.11 haven_2.4.3 [19] scales_1.1.1 tzdb_0.2.0 proxy_0.4-26 [22] farver_2.1.0 generics_0.1.1 ellipsis_0.3.2 [25] withr_2.4.3 cli_3.1.0 magrittr_2.0.1 [28] crayon_1.4.2 readxl_1.3.1 fs_1.5.2 [31] fansi_1.0.2 doParallel_1.0.16 xml2_1.3.3 [34] class_7.3-20 textshaping_0.3.6 tools_4.1.2 [37] prettyunits_1.1.1 hms_1.1.1 lifecycle_1.0.1 [40] munsell_0.5.0 reprex_2.0.1 compiler_4.1.2 [43] e1071_1.7-9 systemfonts_1.0.2 rlang_0.4.12 [46] classInt_0.4-3 units_0.7-2 grid_4.1.2 [49] iterators_1.0.13 rstudioapi_0.13 htmlwidgets_1.5.4 [52] gtable_0.3.0 codetools_0.2-18 DBI_1.1.2 [55] R6_2.5.1 lubridate_1.8.0 rgdal_1.5-28 [58] knitr_1.37 fastmap_1.1.0 utf8_1.2.2 [61] ragg_1.1.3 KernSmooth_2.23-20 stringi_1.7.6 [64] parallel_4.1.2 Rcpp_1.0.8 vctrs_0.3.8 [67] rgl_0.108.3 xfun_0.29 dbplyr_2.1.1 [70] rayshader_0.26.8 tidyselect_1.1.1 ``` ```r > sf_extSoftVersion() GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H "3.10.2" "3.4.1" "8.2.0" "true" "true" PROJ "8.2.0" ``` NOTE: that although PROJ is reported as `8.2.0` above, it is actually >8.2.1, as the version I am testing with is built from source at HEAD.

rsbivand commented 2 years ago

Please confirm that you rebuilt GDAL with new PROJ, as the two are connected. What does sf_extSoftVersion() say? You will need to re-install sf avoiding pre-installation of proj.db on your platform - if you have access to Linux, use rather that than macOS. On Windows and macOS, sf and other packages using PROJ, GDAL, GEOS package basic shared metadata, so very possibly your PROJ HEAD proj.db is simply not seen by sf. My PROJ 8.2.1 shows errors with both sf and rgdal - the latter does not go via GDAL.

> CRS("+proj=peirce_q +lon_0=25 +shape=square")
proj_create: Error 1027 (Invalid value for an argument): peirce_q: peirce_q: invalid value for 'type' parameter
proj_create: Error 1027 (Invalid value for an argument): peirce_q: peirce_q: invalid value for 'type' parameter
Error in CRS("+proj=peirce_q +lon_0=25 +shape=square") : NA
> CRS("+proj=peirce_q +lon_0=25 +type=square")
proj_normalize_for_visualization: Cannot retrieve source or target CRS
proj_normalize_for_visualization: Cannot retrieve source or target CRS
Error in CRS("+proj=peirce_q +lon_0=25 +type=square") : NA

tcwilkinson commented 2 years ago

Thanks for checking this out.

Yes, I've rebuilt GDAL (though to stable, not HEAD/master) with the new PROJ; and re-installed sf and rgdal from within R.
sf_extSoftVersion() as above under Details.
I wondered about the sf not being able to find the HEAD proj.db but since sf_extSoftVersion show 8.2.0 (which I take to mean it is linked to HEAD version, as before build it was linked to current stable, i.e. 8.2.1), am confused as to where sf is looking.

For me, CRS() gives different WKT results from st_crs(), with the HEAD branch (see below). How do each of these generate their WKT strings?

> CRS("+proj=peirce_q +lon_0=25 +shape=square")
Coordinate Reference System:
Deprecated Proj.4 representation:
 +proj=peirce_q +shape=square +lat_0=90 +lon_0=25 +k_0=1 +x_0=0 +y_0=0
+datum=WGS84 +units=m +no_defs 
WKT2 2019 representation:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Peirce Quincuncial (Square)"],
        PARAMETER["Latitude of natural origin",90,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

> st_crs("+proj=peirce_q +lon_0=25 +shape=square")
Coordinate Reference System:
  User input: +proj=peirce_q +lon_0=25 +shape=square 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["PROJ peirce_q"],
        PARAMETER["lon_0",25,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

CRS() above seems to acknowledge the square param as METHOD["Peirce Quincuncial (Square)"], but then it gives the same response for another shape (+shape=horizontal), which is wrong (see below). Again this may well be upstream somewhere, but I don't know where to look to chase it up. Any help gratefully received!

> CRS("+proj=peirce_q +lon_0=25 +shape=horizontal")
Coordinate Reference System:
Deprecated Proj.4 representation:
 +proj=peirce_q +shape=square +lat_0=90 +lon_0=25 +k_0=1 +x_0=0 +y_0=0
+datum=WGS84 +units=m +no_defs 
WKT2 2019 representation:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Peirce Quincuncial (Square)"],
        PARAMETER["Latitude of natural origin",90,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

ps. Have not tried on linux yet... (Some time needed to set up an equivalent working vm...)

rsbivand commented 2 years ago

Briefly, sf_proj_search_paths() shows where the runtime in sf is looking. sf goes through GDAL, so is perhaps losing contact with the updated definition there.

tcwilkinson commented 2 years ago

OK, makes sense, thanks.

sf_proj_search_paths() was searching remnant homebrew/cellar installation of version 7.2.1 I had... but updating these paths made no difference to outputted WKT from sf::st_crs() so problem may well be trip via gdal as you say.

Do you have any idea why i) CRS("+proj=peirce_q +lon_0=25 +shape=square") and ii) CRS("+proj=peirce_q +lon_0=25 +shape=horizontal") both produce the same WKT+proj output (with (ii) being wrong). Is this WKT string created in rgdal, or returned from a call to the PROJ api ? ...I can chase up over on PROJ if the second...

edzer commented 2 years ago

Is this WKT string created in rgdal, or returned from a call to the PROJ api ? ...I can chase up over on PROJ if the second...

It's created by GDAL, which calls PROJ in turn. See https://github.com/r-spatial/sf/blob/main/src/gdal.cpp#L117-L122

rsbivand commented 2 years ago

Building against GDAL current HEAD and PROJ current HEAD:

> st_crs("+proj=peirce_q +lon_0=25 +shape=square")
Coordinate Reference System:
  User input: +proj=peirce_q +lon_0=25 +shape=square 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Peirce Quincuncial (Square)"],
        PARAMETER["Latitude of natural origin",90,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

And yes, +shape=square and +shape=horizonal give the same WKT2, +shape=diamond does seem to feed through. With GDAL current HEAD:

library(sf)
sq <- st_crs("+proj=peirce_q +lon_0=25 +shape=square")
library(rnaturalearth)
sf_world <- ne_download( scale = 110, type = 'countries' , returnclass="sf")
plot(st_geometry(st_transform(sf_world, sq)))

rsbivand commented 2 years ago

dm <- st_crs("+proj=peirce_q +lon_0=25 +shape=diamond")
plot(st_geometry(st_transform(sf_world, dm)))

edzer commented 2 years ago

Looks like an upstream (PROJ) problem that has been solved in PROJ 9.0.0RC1:

> library(sf)
Linking to GEOS 3.10.0, GDAL 3.4.1, PROJ 9.0.0; sf_use_s2() is TRUE
> x1 = st_crs("+proj=peirce_q +lon_0=25 +shape=square")
> x2 = st_crs("+proj=peirce_q +lon_0=25 +shape=horizontal")
> all.equal(x1, x2)
[1] "Component \"input\": 1 string mismatch"
[2] "Component \"wkt\": 1 string mismatch"  
> x1
Coordinate Reference System:
  User input: +proj=peirce_q +lon_0=25 +shape=square 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Peirce Quincuncial (Square)"],
        PARAMETER["Latitude of natural origin",90,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
> x2
Coordinate Reference System:
  User input: +proj=peirce_q +lon_0=25 +shape=horizontal 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["PROJ peirce_q shape=horizontal"],
        PARAMETER["lon_0",25,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
    REMARK["PROJ CRS string: +proj=peirce_q +lon_0=25 +shape=horizontal"]]

r-spatial / sf

Missing params in WKT string generated by st_crs() #1904