Inaccurate and asymmetric (North/South) transformation utm to lat lon

camillechasset commented 7 months ago

There is apparently an error (wrong parenthesis position) in the transformation from UTM to lat / lon, at this line: https://github.com/Turbo87/utm/blob/418d9c477003b1aa18be551f4c46ff689c9c2dc2/utm/conversion.py#L172

Original computation:

    latitude = (p_rad - (p_tan / r) *
                (d2 / 2 -
                 d4 / 24 * (5 + 3 * p_tan2 + 10 * c - 4 * c2 - 9 * E_P2)) +
                 d6 / 720 * (61 + 90 * p_tan2 + 298 * c + 45 * p_tan4 - 252 * E_P2 - 3 * c2))

Proposed change:

    latitude = p_rad - (p_tan / r) *
                (d2 / 2 -
                 d4 / 24 * (5 + 3 * p_tan2 + 10 * c - 4 * c2 - 9 * E_P2) +
                 d6 / 720 * (61 + 90 * p_tan2 + 298 * c + 45 * p_tan4 - 252 * E_P2 - 3 * c2))

One simple thing to notice is that the original correction term in d6 had the same sign in the North and South hemisphere whereas it is expected to change. I attach some plots showing the performance improvement when implementing this change, evaluated by doing the double conversion from utm to lat/lon and then back:

Residual errors excluding Norway exceptions, before fix (left) and after (right), error reduced from up to 4.5mm to 0.8mm (decreased by a factor ~6). Note the North / South asymmetry.
Same, including Norway exception. Error reduced from up to 6.5cm to 6.5mm (decreased by a factor ~10).
Same, including Svalbard exception. Error reduced from up to 8cm to 13mm (decreased by a factor ~6).

I didn't manage to find the source of the code, but on this page you can see a similar series decomposition, and all terms of the latitude expression series are multiplied by beta1*t1, which in your code is (p_tan / r). That's why I believe the change is needed.

camillechasset commented 6 months ago

There is apparently an error (wrong parenthesis position) in the transformation from UTM to lat / lon, at this line:

https://github.com/Turbo87/utm/blob/418d9c477003b1aa18be551f4c46ff689c9c2dc2/utm/conversion.py#L172

Original computation:
    latitude = (p_rad - (p_tan / r) *
                (d2 / 2 -
                 d4 / 24 * (5 + 3 * p_tan2 + 10 * c - 4 * c2 - 9 * E_P2)) +
                 d6 / 720 * (61 + 90 * p_tan2 + 298 * c + 45 * p_tan4 - 252 * E_P2 - 3 * c2))
Proposed change:
    latitude = p_rad - (p_tan / r) *
                (d2 / 2 -
                 d4 / 24 * (5 + 3 * p_tan2 + 10 * c - 4 * c2 - 9 * E_P2) +
                 d6 / 720 * (61 + 90 * p_tan2 + 298 * c + 45 * p_tan4 - 252 * E_P2 - 3 * c2))
One simple thing to notice is that the original correction term in d6 had the same sign in the North and South hemisphere whereas it is expected to change. I attach some plots showing the performance improvement when implementing this change, evaluated by doing the double conversion from utm to lat/lon and then back:

Residual errors excluding Norway exceptions, before fix (left) and after (right), error reduced from up to 4.5mm to 0.8mm (decreased by a factor ~6). Note the North / South asymmetry.

Same, including Norway exception. Error reduced from up to 6.5cm to 6.5mm (decreased by a factor ~10).

Same, including Svalbard exception. Error reduced from up to 8cm to 13mm (decreased by a factor ~6).

I didn't manage to find the source of the code, but on this page you can see a similar series decomposition, and all terms of the latitude expression series are multiplied by beta1*t1, which in your code is (p_tan / r). That's why I believe the change is needed.

@Turbo87 any thoughts about this? :pray:

tvkn commented 4 months ago

Hey @Turbo87! Can we help by opening a PR to fix this? Thanks!

Turbo87 / utm

Inaccurate and asymmetric (North/South) transformation utm to lat lon #120