Write WKT2 string instead of older WKT1 when setting -JEPSG:XXXX

[weiji14]

Description of the desired feature

In modules like surface where a NetCDF grid is produced, the default is to write a OGC WKT1 string when -JEPSG:XXXX is used (see also https://docs.generic-mapping-tools.org/6.1/gmt.html#jproj-full). However, there can be some differences between the older OGC WKT1 format and the newer WKT2 (2018) format, GMT should probably go with the newer format (set using FORMAT=WKT?). Below is the diff of WKT1/WKT2 for Antarctic Polar Stereographic EPSG:3031 (the latitude of origin should be -71S):

 PROJCS["unknown",
     GEOGCS["unknown",
-        DATUM["unknown",
-            SPHEROID["unknown",6378137,298.257220143428]],
+        DATUM["WGS_1984",
+            SPHEROID["WGS 84",6378137,298.257223563,
+                AUTHORITY["EPSG","7030"]],
+            AUTHORITY["EPSG","6326"]],
         PRIMEM["Greenwich",0,
             AUTHORITY["EPSG","8901"]],
         UNIT["degree",0.0174532925199433,
             AUTHORITY["EPSG","9122"]]],
     PROJECTION["Polar_Stereographic"],
-    PARAMETER["latitude_of_origin",-90],
+    PARAMETER["latitude_of_origin",-71],
     PARAMETER["central_meridian",0],
-    PARAMETER["scale_factor",1],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["metre",1,

Relevant code in GMT seems to be here:

gmt/src/gmt_gdalwrite.c

Lines 416 to 420 in b169019

	else if (OSRImportFromProj4(hSRS_2, prhs->P.ProjRefPROJ4) == CE_None) {
	char *pszPrettyWkt = NULL;
	OSRExportToPrettyWkt(hSRS_2, &pszPrettyWkt, false);
	projWKT = pszPrettyWkt;
	}

References:

• https://gdal.org/api/ogrspatialref.html#_CPPv4NK19OGRSpatialReference11exportToWktEPPcPPCKc
• https://gdal.org/programs/gdalsrsinfo.html#cmdoption-gdalsrsinfo-o

Raised originally at GenericMappingTools/pygmt#655 (comment).

Are you willing to help implement and maintain this feature? Yes/No

[joa-quim]

I don't think your diff represents the difference between wkt1 and wkt2. At least looking at the output of a wkt2

gdalsrsinfo EPSG:3031 -o wkt2

PROJCRS["WGS 84 / Antarctic Polar Stereographic",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Antarctic Polar Stereographic",
        METHOD["Polar Stereographic (variant B)",
            ID["EPSG",9829]],
        PARAMETER["Latitude of standard parallel",-71,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8832]],
        PARAMETER["Longitude of origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8833]],
        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)",north,
            MERIDIAN[90,
                ANGLEUNIT["degree",0.0174532925199433]],
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            MERIDIAN[0,
                ANGLEUNIT["degree",0.0174532925199433]],
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Antarctic Digital Database and small scale topographic mapping."],
        AREA["Antarctica."],
        BBOX[-90,-180,-60,180]],
    ID["EPSG",3031]]

[stale]

This issue has been automatically marked as stale because it has not had activity in the last 90 days. It will be closed if no further activity occurs within 7 days. Thank you for your contributions.

[joa-quim]

Closing this as I'm not going to do it and there are issues in GDAL indicating that WKT2 will be the default in a future (soon?) GDAL 4.0