Add support for Universal Polar Stereographic

[DanielJDufour]

For polar regions, mgrs switches from UTM to Universal Polar Stereographic:
https://en.wikipedia.org/wiki/Military_Grid_Reference_System#Polar_regions

For previous discussion see: #11

[Tronald]

US ARMY TEC-SR-7 1996 is the latest reference I have in terms of UPS/ MGRS POLAR conversions. We implemented it in CoordinateSharp, but the precision loss is extremely poor during conversions we have found. It ranges from 0-33m between the 80/84-88th parallels and up to 2.2km above the 88th parallel as you near the poles. Other services we tested against seem to exhibit similar results.

Would love to work on this with you in terms of formulas if you have any advancements or additional resources on this. The polars aren't used often, so there seems to be little in terms of public resources on conversions into and out of the system.

[DanielJDufour]

Hi, @Tronald . Have you taken a look at https://github.com/codice/usng.js/? I haven't tested their UPS conversions yet, but seems promising. I'd also love to contribute to any efforts to provide more consistent and reliable UPS conversions. Personally, I'm trying to get IBCAO 4.0 on a web map and I need a UPS solution to do that. Thanks for reaching out!

[Tronald]

I just reviewed the package and I believe they are experiencing problems with UPS to LL conversions.

The LL to UPS results appear to be consistent with ours at a glance (within a meter). We gauge our precision results by getting the differences of conversion-to-converted back values though.

Precision = (LL -> UPS) - (UPS -> LL)

Their convert back logic doesn't seem to work (not sure what's going on), so I can't even gauge the precision loss. I could just be using it wrong though (I'm weak when it comes to JS).

RunKit Example

[DanielJDufour]

Hmmm... I wonder what could be the cause of it. Have you checked out GEOTRANS and see how they are doing it? I published a container to dockerhub under geodesy/geotrans that you can use to get easy access to the code. You can run it with docker pull geodesy/geotrans; docker run --rm -i -t geodesy/geotrans bash. It has some test files for MGRS and UPS conversion as well as the source code for UPS conversions. You might have already looked at it, but figured couldn't hurt to ask.

[Tronald]

I actually hadn't looked at the GEOTRANS source, so that would be a great resource to have thank you. GEOTRANS 3.8 suffers precision loss as well. It is slightly different than CoordinateSharp (CS), but I wonder if the truncation is being handled different with GEOTRANS 3.8.

Example Problem

WGS 1984 Ellipsoid

CoordinateSharp

LL INITIAL: N 88º 36' 36" E 95º 41' 9.6"
CONVERTED UPS: Z 2153571mE 2015290mN
CONVERTED BACK LL: N 88º 36' 36.003" E 95º 41' 8.9"

PRECISION LOSS: +0.003 seconds N / - 0.7 seconds E

GEOTRANS

LL INITIAL: N 88º 36' 36" E 95º 41' 9.6"
CONVERTED UPS: 2153571 2015291
CONVERTED BACK LL: N 88º 36' 36" E 95º 41' 10.2"

PRECISION LOSS: 0 seconds N / +0.6 seconds E.

As you can see the precision loss between the two is close. It looks as if the full Northing float value after conversion in CS is 2015290.529... Truncating vs Rounding will change that meter square.

This makes me wonder if GEOTRANS 3.8 is rounding vs truncating. I tested against Earthpoint.us and their result was the truncated Northing value that CS exhibits.

If I adjust the Northing by one meter (2015291mN) in CS and convert back to LL the results are almost identical to GEOTRANS.

My unit tests are heavy for UPS and MGRS Polar, but they are ran against Earthpoint with a 1 meter delta to UPS and a .0003 second delta back to LL. I will look into making one to run against GEOTRANS to determine the delta between the formulas of each tool.

If the results show a similar delta and precision loss during conversions, then I think there may just be a limitation between the two systems when converting. I wonder if GEOTRANS has specified conversion precision loss in their documentation anywhere?

[DanielJDufour]

Hi, @Tronald . This is very interesting research. Thank you!

Here's some code from geotrans3.8/SpreadsheetTester/linux_64/SpreadsheetTester.cpp that speaks to precision:

BOOLEAN Compare_MGRS_Result(char* expected_str, char* results_str)
{
  long expected_zone, results_zone;
  long expected_letters[3], results_letters[3];
  double expected_easting, results_easting;
  double expected_northing, results_northing;
  int expected_length = strlen(expected_str);
  int results_length = strlen(results_str);

  if (expected_length != results_length)
    return FALSE;
  else
  {
    if (!Break_MGRS_Str(expected_str, &expected_zone, expected_letters, &expected_easting, &expected_northing))
    {
      if (!Break_MGRS_Str(results_str, &results_zone, results_letters, &results_easting, &results_northing))
      {
        if (expected_zone != results_zone)
          return FALSE;

        if (expected_letters[0] != results_letters[0] ||
            expected_letters[1] != results_letters[1] ||
            expected_letters[2] != results_letters[2])
          return FALSE;

        if (fabs(expected_easting - results_easting) > 1.0e0)
          return FALSE;

        if (fabs(expected_northing - results_northing) > 1.0e0)
          return FALSE;
      }
    }
  }

It appears that the conversions fails GEOTRANS tests if the discrepancy is greater than 1.0e0, which as you know is equivalent to 1 meter. Hope that helps a bit.

It's not a part of the standard AFAIK https://en.wikipedia.org/wiki/Military_Grid_Reference_System, but if the precision loss is an issue when converting to a grid-based reference system, we could discuss providing an option for a result with greater than 5 digits of precision? Obviously, there'd be a pre-processing step required to make this integrate with other libraries. (I'll also want to double-check with the proj4js owners before adding this functionality).

If you're interested, it might also be useful to create a separate library focused just on UPS that this library could support. My thinking is that there's some applications for UPS that don't need all the other MGRS stuff.

Thank you!

[Tronald]

Great Info.

The 1 meter delta makes sense to me. They are also using double E/N values which means they are testing at the centimeter level to avoid precision loss. I wonder where they obtained their test case values? Is the list of test case values provided in the source, I would love to have it (haven't had a chance to look)?

We have the double values available option in CoordinateSharp (as of a year ago) as other users presented similar issues to us. It certainly corrected them so I think it would be smart to include an option for all the grid based systems proj4js includes.

Upon further research it looks like TEC-SR-7 references DMA TM 8358.2. The manual has more detailed formulas and describes convergence and scale factor distortion a bit. This distortion may be what we are seeing and it certainly appears larger with the UPS system. Uncorrected distortion mixed with truncating and rounding rules would certainly make for the precision losses.

I need to run tests in my own library at the centimeter precision level to see how much it corrects. Thanks for bouncing ideas back and forth. I think it will help both our cases.

[spatialillusions]

Since more and more webmap libraries support globe projections, it is easy to move a point closer to the poles and in my case that made my code break. If it would be possible to support UPS in some way it would be great!
I will dig down a bit more in this and see what sort of solution I will do.

[DanielJDufour]

@spatialillusions , I would be very interested in reviewing a PR that adds UPS support to this library. Thank you for taking the initiative on this!

[cffk]

Here are my thoughts on this. This is triggered by the issue submitted
by @spatialillusions at geographiclib/geographiclib-js#5. I am,
otherwise, a disinterested observer.

I have authored code to tackle this problem in two parallel
implementations

• C++
  code: https://github.com/geographiclib/geographiclib
  documentation: https://geographiclib.sourceforge.io/C++/doc
  classes: TransverseMercator PolarStereographic UTMUPS MGRS
  utility: GeoConvert

• Octave/MATLAB
  code: https://github.com/geographiclib/geographiclib-octave
  doc: https://github.com/geographiclib/geographiclib-octave#readme
  functions: {tranmerc,polarst,utmups,mgrs}_{fwd,inv}

A novel feature of these implementations is bringing together UTM and
UPS into one class/function with UPS being designated as zone 0. The
conversions offered by these functions allow the zone and UTM/UPS
selection to be overridden, if, for some reason, you need the
designation in a neighboring zone. (This is crucial in many
applications.)

I also recommend having the projections (transverse Mercator and polar
stereographic) handled separately. It's probably a bad idea to have the
implementation of transverse Mercator wrapped up in the MGRS
machinations as functions LLtoUTM and UTMtoLL. The problems are:

• The chosen algorithm is awful and is known to be inaccurate at high
  latitudes. Perhaps this explains the inconsistency seen in
  Inconsistent Conversion of MGRS to lat/long and vice versa #85.

• it's duplicating functionality in proj4js/proj4js (which uses an
  accurate algorithm).

A few other matters:

• My code supported MSGS precisions up to 1 um.

• I describe and enforce careful rules for the allowed latitude letters
  in cases where a band boundary intersects a 100 km square.

• Round trip calculations are typically consistent.

I'm happy to answer questions (but I won't be contributing code).

[DanielJDufour]

@cffk , thanks for sharing!