Polar stereographic projection with latitude of standard parallel

Hello, thank you for any guidance on this matter.

I want to do a polar stereographic grid projection with distortion minimized at 60N. At first I tried using the -Js option, but it was not clear if this could be done. Instead, I tried a PROJ string, like such:

gmt grdproject rtopo_2.0.4/bed_topography.nc -R-130/25/0/52r -J"+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-94 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -GNorth_America_base_topo.nc -D80000= -Fe -V -r -C

This PROJ string is based off of EPSG:4313. The +lat_ts option should set the standard parallel. However, the resulting WKT text does not indicate the standard parallel was applied:

grdproject [INFORMATION]: Proj4 string to be converted to WKT:
	+proj=stere +lat_0=90 +lat_ts=60 +lon_0=-94 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
grdproject [INFORMATION]: WKT converted from proj4 string:
PROJCS["unknown",
    GEOGCS["unknown",
        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",60],
    PARAMETER["central_meridian",-94],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",SOUTH],
    AXIS["Northing",SOUTH]]


For comparison, this is the WKT text for EPSG:4313 (which has a standard parallel at 70N, note the “Latitude of standard parallel” parameter)

WKT
PROJCRS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",
    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["US NSIDC Sea Ice polar stereographic north",
        METHOD["Polar Stereographic (variant B)",
            ID["EPSG",9829]],
        PARAMETER["Latitude of standard parallel",70,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8832]],
        PARAMETER["Longitude of origin",-45,
            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["easting (X)",south,
            MERIDIAN[45,
                ANGLEUNIT["degree",0.0174532925199433]],
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",south,
            MERIDIAN[135,
                ANGLEUNIT["degree",0.0174532925199433]],
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - N hemisphere - north of 60°N"],
        BBOX[60,-180,90,180]],
    ID["EPSG",3413]]
Proj4
+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

Not sure of this but your second listing uses the WKT2 whilst the first uses WKT1, and it says

PARAMETER["latitude_of_origin",60],

So maybe (likely) it was applied (that’s a work of GDAL/PROJ that are called by GMT) and is the WKT1 report that is less explicit.

And BTW, EPSG:4313 is for Belgium

julia> println(epsg2wkt(4313))
GEOGCS["BD72",
    DATUM["Reseau_National_Belge_1972",
        SPHEROID["International 1924",6378388,297,
            AUTHORITY["EPSG","7022"]],
        AUTHORITY["EPSG","6313"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4313"]]