pyproj4 / pyproj

Python interface to PROJ (cartographic projections and coordinate transformations library)
https://pyproj4.github.io/pyproj
MIT License
1.04k stars 210 forks source link

to_CF for polar_stereographic, with latitude_true_scale, should include latitude_of_projection_origin #1363

Closed autydp closed 8 months ago

autydp commented 8 months ago

Issue is in CRS.to_CF(\<CRS-Object>) method

If Projection-Type (Grid_Mapping_Name) is polar_stereographic, and in the case that latitude_true_scale is defined, the CF attributes is missing the attribute: latitude_of_projection_origin.

While the latitude_of_projection_origin for polar_stereographic might be inferred by the latitude_true_scale, when that is defined, this suggests the latitude_of_projection_origin is not critical and should be left off. There are, however, end-user applications, such as Panoply, that do not work correctly without this attribute.

It is recommended that latitude_of_projection_origin always be included for polar_stereographic projection in the conversion of a CRS to CF grid_mapping attributes. When converting to CF, this is known and can be provided correctly, without conflict.

Attached (see zip archive) are two images and two NetCDF files. The NetCDF files have a grid_mapping variable generated from CRS . to_CF(), for a polar_stereographic projection, with latitude-true-scale set to 70° - one supplemented with latitude_of_projection_origin, and the other as-generated without latitude_of_projection_origin. There are also two image files that show the output of Panoply for these two cases - but shown without the spatial subsetting for better visual reference. (The data files themselves are subset to minimize their size).

While Panoply’s behavior is perhaps not exemplary, it shows what can happen when latitude_of_projection_origin is missing. It is better if CRS.To_CF() can always include this attribute to avoid any confusion.

For PyProj (pyproj/pyproj/crs/_cf1x8.py#_polarstereographic\_to_cf):

def _polar_stereographic__to_cf(conversion):
    """
    http://cfconventions.org/cf-conventions/cf-conventions.html#polar-stereographic
    """
    params = _to_dict(conversion)
    if conversion.method_name.lower().endswith("(variant b)"):
        return {
            "grid_mapping_name": "polar_stereographic",
            ***"latitude_of_projection_origin":***
            ***   -90 if params["latitude_of_standard_parallel"] < 0 else 90,***
            "standard_parallel": params["latitude_of_standard_parallel"],
            "straight_vertical_longitude_from_pole": params["longitude_of_origin"],
            "false_easting": params["false_easting"],
            "false_northing": params["false_northing"],
        }
    return {
        "grid_mapping_name": "polar_stereographic",
        "latitude_of_projection_origin": params["latitude_of_natural_origin"],
        "straight_vertical_longitude_from_pole": params["longitude_of_natural_origin"],
        "false_easting": params["false_easting"],
        "false_northing": params["false_northing"],
        "scale_factor_at_projection_origin": params["scale_factor_at_natural_origin"],
    }

LatPrjOrigin.zip

djhoese commented 8 months ago

Related and possible duplicate: https://github.com/pyproj4/pyproj/issues/1292

autydp commented 8 months ago

Yes, I think we can mark as duplicate, and apologies for not seeing this before submitting my ticket. I will add this info as a comment to #1292.

I do feel this amplifies the issue a bit as we have use-cases here which do not behave well in e.g., Panoply.