Converting from UPS(x, y) to geodetic latitude/longitude

I'm trying to create a program which converts latitude/longitude (wgs84) to UPS (Universal Polar stereographic) coordinates, and then UPS(x, y) to WGS 84. I mean Given UPS(x, y), be able to compute latitude/longitude and vice versa.

Is there any easy way to do this?

• Is this related to Wolfram Mathematica? If so, where are you stuck? – Kuba Apr 6 at 21:42
• I'm not really familiar with this, but you might try looking at GeodesyData in the MMa help index. – Bill Watts Apr 6 at 22:08
• Unfortunately, Mathematica doesn't currently support EPSG codes to use for GeodesyData, but it is still possible to convert these coordinates if you try hard enough (by using a separate library via python, or by implementing the conversion yourself). Would a solution like that work for you? Please, if possible, give more detail in your question, such as example points or your expected use case. – Carl Lange Apr 6 at 22:46

GeoProjectionData has the "UPSNorth" and "UPSSouth" projections. For example, take the same point chosen by Carl:

In:= alert = GeoPosition[{82.5307536, -62.2750895}]
Out= GeoPosition[{82.5308, -62.2751}]

Then you can do:

Out= GeoGridPosition[{1.26494*10^6, 1.61368*10^6}, "UPSNorth"]

Or you can construct a map of the area around the North Pole:

In:= GeoGraphics[GeoMarker[alert], GeoRange -> {{60, 90}, {-180, 180}}, GeoProjection -> "UPSNorth", GeoGridLines -> Automatic, Frame -> True] • I didn't know that! Well, my answer works for arbitrary SRSs so I'll leave it there for now. Nice work! I hope we can get GeoTransformation[{53, -8}, 4326, 5127] sometime soon ;) – Carl Lange Apr 7 at 19:08

You can solve this problem by using a web API to do arbitrary transformations (in this way you can convert to and from any spatial reference system, and are not limited to the systems that WL supports natively)

For instance:

geoposToXYZ[pos_] :=
Module[{ll = QuantityMagnitude[LatitudeLongitude@pos]}, <|
"y" -> ll[], "x" -> ll[]|>]
SetAttributes[geoposToXYZ, Listable]

transform[pointsXYZ_, inEPSG_, outEPSG_] :=
ToExpression /@
Import[URLBuild["http://epsg.io/trans", <|
"data" ->
StringRiffle[{#[["x"]], #[["y"]]} & /@ pointsXYZ, ";", ","],
"s_srs" -> inEPSG,
"t_srs" -> outEPSG
|>], "RawJSON"]

Which you can call by doing, for example:

transform[geoposToXYZ@{Here,
Entity["City", {"Dublin", "Dublin", "Ireland"}]}, 4326, 4326]
{<|"y" -> 53.33, "x" -> -6.25, "z" -> 0.|>, <|"y" -> 53.33, "x" -> -6.25, "z" -> 0.|>}

This particular call will not do any transformation, since we have passed the same spatial reference system code for both. (And as I am in Dublin, both resultant points are the same).

To solve your exact problem, let's say I have a point near the north pole, in this case the northernmost town of Alert, Canada, in the standard Web Mercator projection, which is widely used in GPS devices and so on (and in Mathematica by default).

We know that the EPSG code for the northern UPS zone is 32661, so we can convert:

{<|"y" -> -4.61741·10^6, "x" -> 5.24739·10^7, "z" -> 0.|>}

Well, there's our conversion. We can naturally pass larger numbers of points. To reverse the conversion, simply swap the arguments:

transform[{<|"y" -> -4.61741·10^6, "x" -> 5.24739·10^7, "z" -> 0.|>}, 32661, 4326]
{<|"y" -> 82.5308, "x" -> -62.2751, "z" -> 0.|>}

There are alternatives:

• implement the transformation yourself (only attempt this if you are willing to go a little crazy - the maths involved can be a little complicated)
• install GDAL and use the gdal_transform CLI tool
• install GDAL and use the python bindings and ExternalEvaluate to convert coordinates
• wait for (possibly distant) future versions of mathematica that will support conversion for more EPSG codes
• I have had good luck training a Predict function to convert coordinates, but this requires you to have a dataset of converted coordinates in the first place.