[
https://issues.apache.org/jira/browse/SIS-532?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=17472068#comment-17472068
]
Emmanuel Giasson commented on SIS-532:
--------------------------------------
Wow! That was a very fast response to that issue! It's a bit strange though
that you could not reproduce the results. FYI, my transform object was created
the following way, maybe a bit different that what you did:
{{ DefaultMathTransformFactory factory = new
DefaultMathTransformFactory();}}
{{ ParameterValueGroup values =
factory.getDefaultParameters("Hotine_Oblique_Mercator");}}
{{ values.parameter("...").setValue(...); // for all parameters}}
{{ MathTransform transform =
factory.createParameterizedTransform(values, null);}}
{{ MathTransform reverseTransform = transform.inverse();}}
Then the points to "unproject" were DirectPosition2D objects, with crs=null
within.
But in any case, thanks for fixing it!
Yes, you may add my name, thank you. Should that matter, this very slight
contribution was done as part of my work at "Thales".
> 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)
