WIP: Add EPSG 3031

[ThatcherC]

Addresses issue #239! I have a minimal working implement of the EPSG 3031 CRS. It produces the same forward and backward values as the example given in EPSG guidance note #7-2, although it is currently failing some of the Pkg test due to a type error somewhere.

Here's the EPSG:3031 to WGS84 geodetic lat/lon conversion working forward and backwards, and getting the same values as Proj.jl:

julia> using CoordRefSystems, Proj

julia> epsg3031 = CoordRefSystems.get(EPSG{3031})
PolarStereographicB{-71°, 0°, WGS84Latest}

julia> convert(epsg3031, LatLon(-89, 10))
PolarStereographicB{WGS84Latest} coordinates
├─ x: 18867.758184287268 m
└─ y: 107004.37396749161 m

julia> convert(LatLon, epsg3031(18867.758184, 107004.373967))
GeodeticLatLon{WGS84Latest} coordinates
├─ lat: -89.00000000000477°
└─ lon: 9.999999999895836°

julia> proj = Proj.Transformation("EPSG:4326", "EPSG:3031")
Transformation pipeline
    description: Antarctic Polar Stereographic
    definition: proj=pipeline step proj=axisswap order=2,1 step proj=unitconvert xy_in=deg xy_out=rad step proj=stere lat_0=-90 lat_ts=-71 lon_0=0 x_0=0 y_0=0 ellps=WGS84
    direction: forward

julia> proj(-89, 10)
(18867.758184287333, 107004.37396749199)

julia> inv(proj)(18867.758184, 107004.373967)
(-89.0000000000049, 9.999999999895838)

Right now I've only implement the Polar Stereographic Variant B projection (there's also A and C variants) and I've ignored the north pole case. Is the north pole case and the A and B variants something that'd be valuable to have implemented as well? Currently I don't have a need for them but they'll look a lot like the Polar Stereographic B case so I could try implementing them if there's interest.

I'd love any advice on where to go next with this! Especially if anyone has insight into how to get the failing tests to pass, and whether I should add any more EPSG:3031-specific test.

On another topic, are there docs on what fx and fy are supposed to be? That'd be nice info to add into a guide for contribution. I found these functions a little confusing, particularly the part where they seem to be set up to not input and output dimension quantities.

[juliohm]

That is a very good start @ThatcherC ! 👏🏽🙂

it is currently failing some of the Pkg test due to a type error somewhere.

That is probably due to the type constructor with/without datum. Try to match the code pattern of other CRS types before moving forward.

Is the north pole case and the A and B variants something that'd be valuable to have implemented as well?

It would be super welcome to nail them down now that you are working on it. Check the Ortographic type for an example of a CRS that is parametrized by a Mode categorical parameter. You could follow the same pattern to distinguish the variants and the poles, something like:

struct PolarStereographic{Variant,Hemisphere,lat₁,Datum,Shift,M<:Met} <: Projected{Datum,Shift}
  x::M
  y::M
end

const PolarStereographicNorthA = PolarStereographic{:A,:North}
const PolarStereographicNorthB = PolarStereographic{:A,:North}

add any more EPSG:3031-specific test.

Make sure to add a test against PROJ as we do with other CRS types here

On another topic, are there docs on what fx and fy are supposed to be?

They are the formulas for the projected coordinates, without units and without scaling (multiplication by major axis). We have additional material in GDSJL Chapter 6: https://juliaearth.github.io/geospatial-data-science-with-julia/06-projections.html

[ThatcherC]

Thanks @juliohm! I've pushed a few updates. Removing the datum from the type constructor definitely helped the tests. I'm now seeing issues with the approximate equality tests and the conversion tests against Proj (both of which I think might be a type stability issue). I'll take a look at the A and C variants too - that reference to the Orthographic type is super helpful!

[juliohm]

Should we rename it to PolarStereographic already? I understand B is just one subtype.

[ThatcherC]

Sure! I'll make types for the A, B, and C variants. Should they be subtypes of Mode like EllipticalMode in src/modes.jl? I could also make an abstract type Variant with subtypes A, B, and C (either in src/modes.jl or in a new file like ``src/variants.jl` with a similar layout).

[juliohm]

It would be nice if you could encapsulate these variants in the same source file where you are defining PolarStereographic, assuming these variants are specific to this CRS.

If you think that these variants are general enough to be used by other CRS types, then we can consider a shared file like src/modes.jl or src/variants.jl

[ThatcherC]

Sounds good - I'll keep them in the polarstereographic file for now. Easy to move them later if need be!

[ThatcherC]

Okay, Variant types have been added. I've updated the code so that the PolarStereographic type with the B Variant works, the PolarStereographicB type works, and the Polar Stereographic type with the A and C variants do not work because they are unimplemented (for now):

julia> convert(LatLon, CoordRefSystems.PolarStereographic{CoordRefSystems.B, -71°}(10,10))
GeodeticLatLon{WGS84Latest} coordinates
├─ lat: -89.99986984060983°
└─ lon: 45.0°

julia> convert(LatLon, PolarStereographicB{-71°}(10,10))
GeodeticLatLon{WGS84Latest} coordinates
├─ lat: -89.99986984060983°
└─ lon: 45.0°

julia> convert(LatLon, CoordRefSystems.PolarStereographic{CoordRefSystems.A, -71°}(10,10))
ERROR: MethodError: no method matching formulas(::Type{CoordRefSystems.PolarStereographic{…}}, ::Type{Float64})
The function `formulas` exists, but no method is defined for this combination of argument types.
...

I'm not making use of the North and South modes because according the EPSG guidance notes that "In the
formulas that follow for these variants B and C, the hemisphere of pole is taken to be that of the hemisphere
of the standard parallel." It'll be a little different for Variant A, but I'd prefer to keep this PR focused on PolarStereographicB and EPSG:3031 for now.

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 93.33%. Comparing base (e0db39b) to head (374ec8e).
Report is 2 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #240      +/-   ##
==========================================
+ Coverage   93.12%   93.33%   +0.21%     
==========================================
  Files          34       35       +1     
  Lines        1643     1695      +52     
==========================================
+ Hits         1530     1582      +52     
  Misses        113      113              

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[eliascarv]

According to EPSG Guidance, each variation has different parameters. In this case, the Variants idea will not work.

Probably the most practical idea will be to put the name of the variation in the type, like PolarStereographicVaraintB.

[eliascarv]

Alternatively, the Polar Stereographic Projection can be implemented generically using the Stereographic Projection, as PROJ does.

This is the PROJ string for the EPSG:3031:

Note that PROJ defines South Pole Stereographic Projection by setting lat_0=-90.

Formula reference: https://neacsu.net/docs/geodesy/snyder/5-azimuthal/sect_21/#formulas-for-the-ellipsoid
PROJ implementation: https://github.com/OSGeo/PROJ/blob/master/src/projections/stere.cpp

[juliohm]

Alternatively, the Polar Stereographic Projection can be implemented generically using the Stereographic Projection, as PROJ does.

I tend to favor these general approaches, and agree with @eliascarv that it would be great if we could introduce a general Stereographic type that covers more cases than just the polar case. I understand that this is much more work than what you had planned @ThatcherC, so no worries if you don't have the time to work on it.

Do you think it would be possible to refactor the PR to accommodate this more general design, but only implement the part that you care about for now? That way we avoid introducing a PolarStereographicB type and removing it in the future after the more general Stereographic type is implemented for all cases.

In the meantime, I would like to thank you all for the great research and work that has been done here. It is amazing to see progress and to learn about new projections as we further develop the package.