[
https://issues.apache.org/jira/browse/SIS-532?page=com.atlassian.jira.plugin.system.issuetabpanels:all-tabpanel
]
Emmanuel Giasson updated SIS-532:
---------------------------------
Description:
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 can be
* exactly +/- 1.0 (theoretical bounds) which causes a division by zero on
(1-U)
* 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.
When U tends toward +/- 1.0, the reciprocal with H tends towards 0, giving {{φ
= pi/2.}} This limit case is not mentioned in EPSG guidance notes (I looked at
p.62 of 2019/09 revision), but it is 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.
Checking if abs(U) >= 1.0 then forcing {{φ}} to pi/2 (and {{λ}} to any valid
value) should be enough to fix the bug.
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)
was:
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 can be
* exactly +/- 1.0 (theoretical bounds) which causes a division by zero on
(1-U)
* 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.
When U tends toward +/- 1.0, the reciprocal with H tends towards 0, giving {{φ
= pi/2.}} This limit case is not mentioned in EPSG guidance notes (I looked at
p.62 of 2019/09 revision), but it is 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.
Checking if abs(U) >= 1.0 then forcing {{φ}} to pi/2 (and {{λ}} to any valid
value) should be enough to fix the bug.
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)
> 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
> Priority: Major
>
> 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 can be
> * exactly +/- 1.0 (theoretical bounds) which causes a division by zero on
> (1-U)
> * 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.
> When U tends toward +/- 1.0, the reciprocal with H tends towards 0, giving
> {{φ = pi/2.}} This limit case is not mentioned in EPSG guidance notes (I
> looked at p.62 of 2019/09 revision), but it is 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.
> Checking if abs(U) >= 1.0 then forcing {{φ}} to pi/2 (and {{λ}} to any valid
> value) should be enough to fix the bug.
> 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)
