Mapping polar stereographic projection to CF

[dopplershift]

Code Sample, a copy-pastable example if possible

Running this code with PyProj 3.0.1

# Your code here
proj_parms = {
    'proj': 'stere',
    'lat_0': 90.0,
    'lat_ts': 60.0,
    'lon_0': -150.0,
    'R': 6371229.0
}
crs = pyproj.CRS.from_dict(proj_parms)
cf = crs.to_cf()

Yields:

{'crs_wkt': 'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["unknown",ELLIPSOID["unknown",6371229,0,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Polar Stereographic (variant B)",ID["EPSG",9829]],PARAMETER["Latitude of standard parallel",60,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8832]],PARAMETER["Longitude of origin",-150,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)",south,MERIDIAN[90,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",south,MERIDIAN[180,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]],ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]',
 'semi_major_axis': 6371229.0,
 'semi_minor_axis': 6371229.0,
 'inverse_flattening': 0.0,
 'reference_ellipsoid_name': 'unknown',
 'longitude_of_prime_meridian': 0.0,
 'prime_meridian_name': 'Greenwich',
 'geographic_crs_name': 'unknown',
 'horizontal_datum_name': 'unknown',
 'projected_crs_name': 'unknown',
 'grid_mapping_name': 'polar_stereographic',
 'standard_parallel': 60.0,
 'straight_vertical_longitude_from_pole': -150.0,
 'false_easting': 0.0,
 'false_northing': 0.0}

Problem description

Based on the CF Convention 1.8 document for polar stereographic, this should contain a latitude_of_projection_origin attribute that has a value of +/-90.

Expected Output

The mapping for the code above should be:

{
...
'latitude_of_projection_origin': 90.0,
...
}

[snowman2]

Conversion reference: https://github.com/pyproj4/pyproj/blob/6127f7a6db680b159edac62dccf54538065ae123/pyproj/crs/_cf1x8.py#L539

[snowman2]

The WKT has variant B in it.

[snowman2]

Looks like the WKT is missing the information:

>>> print(crs.coordinate_operation.to_wkt(pretty=True))
CONVERSION["unknown",
    METHOD["Polar Stereographic (variant B)",
        ID["EPSG",9829]],
    PARAMETER["Latitude of standard parallel",60,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8832]],
    PARAMETER["Longitude of origin",-150,
        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]]]

It is in the PROJ string:

>> print(crs.coordinate_operation.to_proj4(pretty=True))
+proj=pipeline
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-150 +x_0=0 +y_0=0 +R=6371229

[snowman2]

Looks like it does not exist in the PROJ proj_create_conversion_polar_stereographic_variant_b: https://github.com/OSGeo/PROJ/blob/7b407e36e650aeae986218a0e213b2d8248c008d/src/iso19111/c_api.cpp#L6596-L6600

[snowman2]

Looks like it isn't supposed to be there based on this: https://pro.arcgis.com/en/pro-app/latest/help/mapping/properties/stereographic.htm

Latitude Of Origin is only on variant A.

[snowman2]

Related #592

[snowman2]

@dopplershift I am thinking this is potentially a bug with the CF convention documentation.

[dopplershift]

Hrmmm...without it, how do you properly map the polar_stereographic to e.g. PROJ? Don't you need to be able to give that to +proj=stere?

[snowman2]

See #534. Pyproj uses the WKT/PROJ JSON to do the conversions.

[snowman2]

This is what is used internally: https://pyproj4.github.io/pyproj/stable/build_crs.html

[snowman2]

@dopplershift I believe pyproj is handling it correctly.

If I do:

from pyproj import CRS
cf_dict = {'semi_major_axis': 6371229.0, 'semi_minor_axis': 6371229.0, 'inverse_flattening': 0.0, 'reference_ellipsoid_name': 'unknown', 'longitude_of_prime_meridian': 0.0, 'prime_meridian_name': 'Greenwich', 'geographic_crs_name': 'unknown', 'horizontal_datum_name': 'unknown', 'projected_crs_name': 'unknown', 'grid_mapping_name': 'polar_stereographic', 'standard_parallel': 60.0, 'straight_vertical_longitude_from_pole': -150.0, 'false_easting': 0.0, 'false_northing': 0.0}
crs = CRS.from_cf(cf_dict)

I get:

<Projected CRS: {"$schema": "https://proj.org/schemas/v0.2/projjso ...>
Name: unknown
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown
- method: Polar Stereographic (variant B)
Datum: unknown
- Ellipsoid: unknown
- Prime Meridian: Greenwich

Which is correct.

What I suggest is to do what is done internally by PROJ:

crs.to_dict()

which outputs:

{'proj': 'stere', 'lat_0': 90, 'lat_ts': 60, 'lon_0': -150, 'x_0': 0, 'y_0': 0, 'R': 6371229, 'units': 'm', 'no_defs': None, 'type': 'crs'}

I think in this scenario, lat_0 of 90 is output by PROJ to ensure that it is polar stereographic even though lat_0 wasn't in the input. I don't think lat_0 has any other purpose for variant B in the PROJ parameters.

If you are using CF to map to PROJ parameters, I think you can safely do this to make cartopy happy with the Variant B scenario (until it can update to support WKT):

lat_0 = cf_dict.get("latitude_of_projection_origin", 90) 

If you find something that indicates that the logic in pyproj is incorrect, please reopen or post to this issue thread to discuss this further.