How to go about adding support for Antarctic Polar Stereo (EPSG:3031) and NSIDC Sea Ice Polar Stereographic North (EPSG:3413)

[alex-s-gardner]

I'd like to add support for two of the most used polar projections for Earth Science (EPSG:3031 & EPSG:3413. This only differ slightly from the Universal Polar Stereographic:

EPSG:3413 Latitude of True Origin = 70 Scale factor at longitude of true origin = 1 Central Meridian = -45

EPSG:3013 Latitude of True Origin = -71 Scale factor at longitude of true origin = 1 Central Meridian = 0

I was hoping to be able to just modify polar_stereographic.jl but it seems that everything is precomputed in transverse_mercator.jl without an easy way to modify Latitude of True Origin or Central Meridian. Before I go down a rabbit hole I was hoping that someone that has looked at this code before might provide some guidance as to the steps I should take to add these two projections.

[asinghvi17]

Any reason you couldn't use Proj.jl? It should handle those projections directly from the EPSG codes. Of course if you need it in pure Julia it's a different story.

[alex-s-gardner]

Ya, Proj.jl does a great job but I was hoping to get improved performance out of a pure Julia implementation (multithreading?). Projections are often a costly part of my large workflows. But maybe Proj can't be improved upon.

[c42f]

someone that has looked at this code before might provide some guidance as to the steps I should take to add these two projections

Can GeographicLib do these? The transformation code in this repo is largely a port of the GeographicLib C++ code, so I'd suggest just porting any code from there where possible. (I can't overstate how meticulous the GeographicLib code is. It's a work of numerical art :-) )

[asinghvi17]

You can multithread using Proj contexts, one per thread. I'll see if I can dig up some code I used.

[alex-s-gardner]

@asinghvi17 did you happen to dig up that code? Would be great to see.

[asinghvi17]

ctxs = [Proj.proj_context_clone() for _ in 1:Threads.nthreads()]
transforms = [Proj.Transformation("+proj=longlat", "EPSG:3031"; ctx) for ctx in ctxs]
# in your loop, access `transforms[Threads.threadid()]` rather than `transform::Proj.Transformation`.

This will work with threads (and you can add your own Proj context in ctxs), but not on the GPU - that would pretty much require a Julia implementation of Proj.