[jira] [Closed] (SIS-532) NaN from unhandled case in reverse Oblique Mercator calculations

["Martin Desruisseaux (Jira)" <ji...@apache.org>]

     [ https://issues.apache.org/jira/browse/SIS-532?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel ]

Martin Desruisseaux closed SIS-532.
-----------------------------------

> NaN from unhandled case in reverse Oblique Mercator calculations
> ----------------------------------------------------------------
>
>                 Key: SIS-532
>                 URL: https://issues.apache.org/jira/browse/SIS-532
>             Project: Spatial Information Systems
>          Issue Type: Bug
>          Components: Referencing
>    Affects Versions: 1.1
>            Reporter: Emmanuel Giasson
>            Assignee: Martin Desruisseaux
>            Priority: Major
>             Fix For: 1.2
>
>
> Reversing 2D coordinates into long,lat that happen to be very near or right on North/South pole may produce "NaN" latitude, instead of +/- 90 degrees, because of an unhandled limit case and/or unavoidable rounding errors on the last bits of finite precision numbers. 
> In release 1.1, org/apache/sis/referencing/operation/projection/ObliqueMercator.java, lines 404-406
> {{  final double U  = (V*cosγ0 + S*sinγ0) / T;}}
> {{  final double λ  = -atan2( (S*cosγ0 - V*sinγ0), cos(y ) ) / B;}}
> {{  final double φ  = φ(pow(H / sqrt((1 + U) / (1 - U)), 1/B));}}
> Value for U is typically within ]-1, 1[ but can be
>  * exactly +/- 1.0 (theoretical bounds) which causes a division by zero on (1-U) or H/sqrt(...)
>  * slightly beyond +/- 1.0 (like 1.0000000000000002) owing to previous operations on finite-precision numbers. This then causes a NaN on the square root of a negative number.
> Your code follows EPSG guidance notes and indeed, the latter makes no mention of that limit case (I looked at p.62 of 2019/09 revision). But it is noted in Snyder, 1987 ([https://pubs.usgs.gov/pp/1395/report.pdf),] from which EPSG guidance notes were most probably copied. See p.75, after equation 9-47. In the current case, we also have to cater for the "beyond +/- 1.0" case because of imprecisions.
> Checking if abs(U) >= 1.0 then forcing {{φ}} to +/- pi/2 (and {{λ}} to any valid value) should thus be enough to fix this issue.
> Should you need an actual case to reproduce the above issue, I stumbled upon it with the following transform:
> Inverse_MT[
>   Param_MT["Hotine Oblique Mercator (variant A)",
>     Parameter["semi_major", 6378137.0, Unit["metre", 1]],
>     Parameter["semi_minor", 6356752.314245179, Unit["metre", 1]],
>     Parameter["Latitude of projection centre", 89.8, Unit["degree", 0.017453292519943295]],
>     Parameter["Longitude of projection centre", 179.8, Unit["degree", 0.017453292519943295]],
>     Parameter["Azimuth of initial line", -174.84239553505424, Unit["degree", 0.017453292519943295]]]]
> While trying to reverse POINT(-905672.3035667514 -1.0011576308445137E7)
> into long,lat, which gave POINT(41.98971287679859 NaN)
>  



--
This message was sent by Atlassian Jira
(v8.20.10#820010)