[GRASS-dev] [GRASS GIS] #2918: Infinite loop in i.landsat.toar with method=DOS4

[GRASS GIS]

#2918: Infinite loop in i.landsat.toar with method=DOS4
---------------------------+-------------------------
 Reporter:  DmitryKolesov  |      Owner:  grass-dev@…
     Type:  defect         |     Status:  new
 Priority:  normal         |  Milestone:  7.0.4
Component:  Imagery        |    Version:  unspecified
 Keywords:                 |        CPU:  Unspecified
 Platform:  All            |
---------------------------+-------------------------
 == Description: ==
 I.landsat.toar with method=DOS4 fails to stop. For example if run the
 module using Landsat8 scene (id=LC81680222013351LGN00) as input, the
 module doesn't finish the calculation.

 == Details: ==
 The problem cycle is the loop on lines 109--115
 (http://trac.osgeo.org/grass/browser/grass/branches/releasebranch_7_0/imagery/i.landsat.toar/landsat.c#L109).
 This is an iterative process for approximate atmospheric transmittance
 coefficients TAUv and TAUz.

 Math details of the approximation are described in <<Song, C., Woodcock,
 C. E., Seto, K. C., Lenney, M. P., & Macomber, S. A. (2001).
 Classification and change detection using Landsat TM data: when and how to
 correct atmospheric effects?. Remote sensing of Environment, 75(2),
 230-244.>> (https://scholar.google.ru/scholar?cluster=2327996098064583125,
 see equations 5 -- 10).

 But the process can't stabilize if solar elevation angle is small (for
 example for winter scenes in North semisphere): in this case Tz (line 113,
 http://trac.osgeo.org/grass/browser/grass/branches/releasebranch_7_0/imagery/i.landsat.toar/landsat.c#L113)
 is negative because of division:

 {{{
 Tz = 1 - (4. * pi_d2 * Lp) / (lsat->band[i].esun * sin_e)
 }}}

 There pi_d2, lsat->band[i].esun are constants, Lp is finite positive
 number, but sin_e can be very small (if solar elevation angle is small).
 As result, Tz is negative, log(Tz) is NaN
 (http://trac.osgeo.org/grass/browser/grass/branches/releasebranch_7_0/imagery/i.landsat.toar/landsat.c#L114),
 so Tv is Nan and the exit condition isn't reachable.


 == Proposed solution: ==
 I see two possibilities if Tz becomes NaN:
  * emit fatal-error like "DOS4 method isn't applicable for the scene. Use
 other method or use other sun_elevation parameter." and exit.
  * emit warning, then reset Tz and Tv, using default values (1.0) for the
 coefficients
 (http://trac.osgeo.org/grass/browser/grass/branches/releasebranch_7_0/imagery/i.landsat.toar/landsat.c#L122).

 I think I could write a patch, but I need an advise about solution.

--
Ticket URL: <http://trac.osgeo.org/grass/ticket/2918>
GRASS GIS <https://grass.osgeo.org>