[
https://issues.apache.org/jira/browse/SIS-532?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17472096#comment-17472096
]
Martin Desruisseaux commented on SIS-532:
-----------------------------------------
I could not reproduce with above numbers, but could reproduce with slightly
different numbers. It may be because trigonometric functions are not guaranteed
to produce the same numbers on all machines. The Javadoc specification for
{{Math.sin(double)}}, {{cos}}, _etc._ explicitly allows for 1 ULP variations
(for allowing the use of different hardware I guess), and I observed in the
past that machine-dependent or JVM-dependent variations does happen.
Thanks for the authorization. I will acknowledge with {{@author Emmanuel
Giasson (Thales)}} Javadoc tag then, if no objection.
> 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.1#820001)
