peirce_q - proj string to WKT convertion not working for all +shapes

[tcwilkinson]

Not all peirce_q methods seem to be convertible to (or findable in) pseudo-WKT2 form from proj string.

Noticed this issue while attempting to test peirce_q from current post-8.2.1 HEAD/master branch with external software which uses PROJ (R-spatial/sf and Rspatial/rgdal).

*Creating CRS using a proj string using +shape=horizontal, +shape=vertical, +shape=shemisphere silently returns METHOD["Peirce Quincuncial (Square)"]. cf. r-spatial/sf#1904 *

For example, run in R, this seems ok-

> library("rgdal")
> rgdal::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)"],              <-- ***CORRECT***
        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]]],
[...snip...]

But this is not correct (and identical to above), as it assumes horizontal is the same as square:-

> rgdal::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   <-- ***INCORRECT***
WKT2 2019 representation:
PROJCRS["unknown",
[...snip...]
    CONVERSION["unknown",
        METHOD["Peirce Quincuncial (Square)"],                                     <-- ***INCORRECT***
        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]],
[...snip...]

Good news +shape=diamond works as expected.

Problem description

This effect is problematic when using external tools which are increasingly converting proj strings to WKT format before accessing further proj/gdal tools/api, as there is no way to specify the other methods (afaik).

• Suspect this is caused by a call to API proj_create using the proj string from rgdal and it is PROJ's proj_create which returns the faulty WKT2 conversion (maybe a test for "diamond" and then otherwise assigns "square") but not tracked down where.

• The current convertion to WKT was primarily to allow matching of square and diamond to ESRI's equivalent methods (introduced during other fixes: peirce_q: rename +type parameter wrongly introduced in 8.2.1 to +shape (fixes #3011) #3014 ), but it seems a pity to restrict only to ESRI options. It also makes it harder to access alternative projections such as Grieger Triptychial, which are based on peirce_q.

Expected Output

IMHO:-

• Perhaps most flexibly, the parameter +shape should be passable directly within proj's WKT (and generated), using suitable a PARAMETER declaration. Then METHOD would be simply be ["Peirce Quincuncial"] for all.
• Or alternatively, additional METHODs would be needed to cover all +shape options in WKT.
• At very least, a proj string of "+proj=peirce_q +lon_0=25 +shape=horizontal" should not silently transform into +shape=square/METHOD["Peirce Quincuncial (Square)"]

e.g. Something along lines of:-

        METHOD["Peirce Quincuncial"],
        PARAMETER["Shape","diamond/square/horizontal/vertical/etc"],
OR
        METHOD["Peirce Quincuncial (Diamond/Square/Horizontal/Vertical)"],

Would be happy to contribute again with guidance, just unclear exactly how best to do so!

Environment Information (fwiw)

• PROJ version (proj) - current HEAD (reports as 8.2.0, but is actually post-8.2.1)
• Operation System Information - currently testing on macOS

Installation method (fwiw)

• source via homebrew (brew install proj --head --build-from-source etc.)
• R packages installed from source to link to

[rouault]

PARAMETER["Shape","diamond/square/horizontal/vertical/etc"],

This would be invalid WKT. The BNF of the WKT CRS spec mandates per http://docs.opengeospatial.org/is/18-010r7/18-010r7.html#77 that the parameter value is numeric.

With #3057, we'll get the generic PROJ WKT form when there's no known (most EPSG driven) official WKT:

$ bin/projinfo "+proj=peirce_q +shape=horizontal +datum=WGS84 +units=m +no_defs +type=crs"
PROJ.4 string:
+proj=peirce_q +shape=horizontal +datum=WGS84 +units=m +no_defs +type=crs

WKT2:2019 string:
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"]],
    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]]]]

$ bin/projinfo '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"]],
>     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]]]]'
PROJ.4 string:
+proj=peirce_q +shape=horizontal +datum=WGS84 +units=m +no_defs +type=crs