`polar_xy_to_lonlat` does not return expected latitudes for NSIDC South Polar Stereographic on the Hugh's 1980 ellipsoid

[andypbarrett]

Problem: polar_xy_to_lonlat does not return the expected latitudes for the Southern PS grid on Hugh's 1980 ellipsoid given in Table 7 and Figure 1 of A Guide to NSIDC's Polar Stereographic Projection, when the correct latitude_of_true_scale is used (-70 N).

What I expect to happen: Latitudes of the 6.25 km grid on the NSIDC South Polar Stereographic on the Hugh's 1980 ellipsoid should be -39.23 N for the upper left and right corners and -41.45 N for the lower left and right grid corners. What I get: Latitudes returned by polar_xy_to_lonlat for the grid corners are -82.88 N and -88.36 N. The longitudes appear to be correct.

Testing: A minimum reproducible code example to calculate longitude and latitude of grid corners is given below. I demonstrate the problem with polar_xy_to_lonlat and also use pyproj to perform the same conversion. The polar_xy_to_lonlat example fails, but the pyproj passes. polar_xy_to_lonlat returns the correct corner latitudes when latitude_true_scale=70. is used.

Suggested Action: The code needs to be corrected to return the correct corner geographic coordinates for the correct latitude_true_scale. Looking at tests.py, they test for internal consistency but not for the correct grid geographic coordinates. I suggest also adding a test that checks the code returns the expected values.

import numpy as np

from pyproj import CRS, Transformer

from polar_convert.constants import SOUTH
from polar_convert import polar_xy_to_lonlat

# TEST DATA
# Input
# Corners of South PS 6.25 km grid in km
xmin, ymax = -3950., 4350.  # Upper left
xmax, ymin = 3950., -3950.  # Lower Right
xcorner = [xmin, xmax, xmax, xmin]
ycorner = [ymax, ymax, ymin, ymin]
corners = ['Upper Left', 'Upper Right', 'Lower Right', 'Lower Left']

# Expected lon, lat for corners
# From Table 7 https://nsidc.org/data/user-resources/help-center/guide-nsidcs-polar-stereographic-projection#anchor-2
expected_lonlat = [(317.76, -39.23), (42.24, -39.23), (135.00, -41.45), (225.00, -41.45)]

# polar_xy_to_lonlat requires latitude of true scale, semimajor axis length and eccentricity.
# These are defined below.
#
# Parameters for Southern Hemisphere Projection Based on Hughes 1980 Ellipsoid (Table 2)
# https://nsidc.org/data/user-resources/help-center/guide-nsidcs-polar-stereographic-projection#anchor-0
semimajor_axis = 6378273
flattening = 1 / 298.279411123064
latitude_true_scale = -70.

# Calculate eccentricity
semiminor_axis = semimajor_axis * (1 - flattening)
eccentricity = np.sqrt(1 - semiminor_axis**2/semimajor_axis**2)


# pyproj uses either a proj4 string, WKT or EPSG code to define CRS.  The EPSG for the North and 
# South Polar Stereographic projection on the Hugh's 1980 ellipsoid has been depreciated.  To avoid any 
# issues that may arise from this depreciation, I use a WKT.  WKT provide a more complete definition of 
# the CRS than a proj4 string.
# From https://epsg.io/3412
WKT = """PROJCS["NSIDC Sea Ice Polar Stereographic South (deprecated)",
    GEOGCS["Unspecified datum based upon the Hughes 1980 ellipsoid (deprecated)",
        DATUM["Not_specified_based_on_Hughes_1980_ellipsoid",
            SPHEROID["Hughes 1980",6378273,298.279411123064,
                AUTHORITY["EPSG","7058"]],
            AUTHORITY["EPSG","6054"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4054"]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",-70],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","3412"]]"""

# Test using polar_xy_to_lonlat
print("Using polar_xy_to_lonlat")
for s, x, y, expected, in zip(corners, xcorner, ycorner, expected_lonlat):
    lon, lat = polar_xy_to_lonlat(x, y, latitude_true_scale, semimajor_axis, eccentricity, SOUTH)
    print(f"{s:11} Result: {lat:6.2f} N, {lon:6.2f} E   "  
    f"Expected: {expected[1]:6.2f} N {expected[0]:6.2f}    "
    f"Matches: {np.allclose((lon, lat), expected, atol=1e-2)}")
print()

# Test using pyproj
print("Using pyproj")
nsidc_ps_crs = CRS.from_wkt(WKT)
wgs84_crs = CRS.from_epsg(4326)
transformer = Transformer.from_crs(nsidc_ps_crs, wgs84_crs, always_xy=True)
for s, x, y, expected in zip(corners, xcorner, ycorner, expected_lonlat):
    lon, lat = transformer.transform(x*1e3, y*1e3)
    lon = 360+lon if lon < 0 else lon  # Convert lon from -180 to 180, to 0 to 360.
    print(f"{s:11} Result: {lat:6.2f} N, {lon:6.2f} E   "  
    f"Expected: {expected[1]:6.2f} N {expected[0]:6.2f}    "
    f"Matches: {np.allclose((lon, lat), expected, atol=1e-2)}")

Using polar_xy_to_lonlat and latitude_true_scale=-70.0
Upper Left  Result:  82.49 N, 317.76 E   Expected: -39.23 N 317.76    Matches: False
Upper Right Result:  82.49 N,  42.24 E   Expected: -39.23 N  42.24    Matches: False
Lower Right Result:  82.10 N, 135.00 E   Expected: -41.45 N 135.00    Matches: False
Lower Left  Result:  82.10 N, 225.00 E   Expected: -41.45 N 225.00    Matches: False

Using pyproj
Upper Left  Result: -39.23 N, 317.76 E   Expected: -39.23 N 317.76    Matches: True
Upper Right Result: -39.23 N,  42.24 E   Expected: -39.23 N  42.24    Matches: True
Lower Right Result: -41.45 N, 135.00 E   Expected: -41.45 N 135.00    Matches: True
Lower Left  Result: -41.45 N, 225.00 E   Expected: -41.45 N 225.00    Matches: True

When latitude_true_scale=70.

Using polar_xy_to_lonlat and latitude_true_scale=70.0
Upper Left  Result: -39.23 N, 317.76 E   Expected: -39.23 N 317.76    Matches: True
Upper Right Result: -39.23 N,  42.24 E   Expected: -39.23 N  42.24    Matches: True
Lower Right Result: -41.45 N, 135.00 E   Expected: -41.45 N 135.00    Matches: True
Lower Left  Result: -41.45 N, 225.00 E   Expected: -41.45 N 225.00    Matches: True

Using pyproj
Upper Left  Result: -39.23 N, 317.76 E   Expected: -39.23 N 317.76    Matches: True
Upper Right Result: -39.23 N,  42.24 E   Expected: -39.23 N  42.24    Matches: True
Lower Right Result: -41.45 N, 135.00 E   Expected: -41.45 N 135.00    Matches: True
Lower Left  Result: -41.45 N, 225.00 E   Expected: -41.45 N 225.00    Matches: True