i.atcorr's create_iwave.py broken

[neteler]

I'm trying to add WorldView3 support but at time the needed helper script does not behave...:

The issue (trivial?) is reproducible with landsat8 as well:

cd imagery/i.atcorr/

python create_iwave.py sensors_csv/landsat_8.csv
Getting sensor name from csv file: landsat_8
Traceback (most recent call last):
  File "create_iwave.py", line 250, in <module>
    main()
  File "create_iwave.py", line 237, in main
    write_cpp(bands, values, sensor, os.path.dirname(inputfile))
  File "create_iwave.py", line 156, in write_cpp
    fi, li = interpolate_band(values[:,[0,b+1]])
  File "create_iwave.py", line 94, in interpolate_band
    filter_f = f(np.arange(response[0,0], response[-1,0] + 2.5, 2.5))
  File "/usr/lib64/python2.7/site-packages/scipy/interpolate/polyint.py", line 79, in __call__
    y = self._evaluate(x)
  File "/usr/lib64/python2.7/site-packages/scipy/interpolate/interpolate.py",
line 477, in _evaluate
    out_of_bounds = self._check_bounds(x_new)
  File "/usr/lib64/python2.7/site-packages/scipy/interpolate/interpolate.py",
line 507, in _check_bounds
    raise ValueError("A value in x_new is above the interpolation "
ValueError: A value in x_new is above the interpolation range.

[martinl]

The problem is that the list of values which is returned by np.arange contains values above range defined by response[:,0].

response[0,0] ->
[ 427.  428.  429.  430.  431.  432.  433.  434.  435.  436.  437.  438.
  439.  440.  441.  442.  443.  444.  445.  446.  447.  448.  449.  450.
  451.  452.  453.  454.  455.  456.  457.  458.  459.]
response[:,1] ->
[  7.30000000e-05   6.09000000e-04   1.62800000e-03   3.42100000e-03
   8.01900000e-03   2.47670000e-02   8.56880000e-02   2.54149000e-01
   5.17821000e-01   7.65117000e-01   9.08749000e-01   9.58204000e-01
   9.77393000e-01   9.83790000e-01   9.89052000e-01   9.86713000e-01
   9.93683000e-01   9.93137000e-01   1.00000000e+00   9.96969000e-01
   9.82780000e-01   9.72692000e-01   9.05808000e-01   7.45606000e-01
   4.71329000e-01   2.26412000e-01   9.28600000e-02   3.66030000e-02
   1.45370000e-02   5.82900000e-03   2.41400000e-03   9.84000000e-04
   2.55000000e-04]
np.arange(response[0,0], response[-1,0] + 2.5, 2.5) ->
[ 427.   429.5  432.   434.5  437.   439.5  442.   444.5  447.   449.5
  452.   454.5  457.   459.5]

459.5 is outside of range <427,459>

[Nikos Alexandris]

In ​https://trac.osgeo.org/grass/browser/grass/trunk/imagery/i.atcorr/create_iwave.py?rev=69225#L93,

the class interp1d is fed with response[:,0] and response[:,1] as x and y respectively (​scipy/interpolate/interpolate.py#L335 and Line 337).

Internally, lower and upper bounds of x and y are defined and ​checked when calling later on the interpolation function returned by interp1d. In create_iwave.py ​`2.5` is added to `y`, which, obviously, leads to the ValueError.

See also the ​example in the classe's source code.

[Nikos Alexandris]

Hopefully the attached create_iwave_patch.diff addresses correctly what r44576 was meant to fix.

[martinl]

Replying to nikosa:

Hopefully the attached create_iwave_patch.diff addresses correctly what r44576 was meant to fix.

note: just remove print statements or replace them by grass.message() or grass.debug().

[Nikos Alexandris]

Patch updated to handle, now hopefully, all cases -- whether the *input* step in the wavelength slots is less, equal or greater than the *hardcoded* step 2.5 for the interpolation function.

@martinl I'd like to add g.messages, printed only per the user's request. However, this being just a helper script, I would only like to make it more readable and informative.

[Nikos Alexandris]

By the way, the csv files for PRISM-B, -F and -N require some edit. The header should be one line, the first column's name being "WL(nm)", as far as I understand the instructions in the README file. Maybe edit, as well, the header of the worldview3.csv file.

It would be nice to have a short section for each sensor in the README file.

[neteler]

Replying to nikosa:

Patch updated to handle, now hopefully, all cases -- whether the *input* step in the wavelength slots is less, equal or greater than the *hardcoded* step 2.5 for the interpolation function.

Thanks for that. Unfortunately it still fails:

python create_iwave.py sensors_csv/landsat_8.csv
Getting sensor name from csv file: landsat_8
Number of bands found: 9
Traceback (most recent call last):
  File "create_iwave.py", line 267, in <module>
    main()
  File "create_iwave.py", line 254, in main
    write_cpp(bands, values, sensor, os.path.dirname(inputfile))
  File "create_iwave.py", line 173, in write_cpp
    fi, li = interpolate_band(values[:,[0,b+1]])
  File "create_iwave.py", line 108, in interpolate_band
    filter_f = f(np.arange(start, stop, step))
  File "/usr/lib/python2.7/dist-packages/scipy/interpolate/polyint.py", line 54, in __call__
    y = self._evaluate(x)
  File "/usr/lib/python2.7/dist-packages/scipy/interpolate/interpolate.py", line 448, in _evaluate
    out_of_bounds = self._check_bounds(x_new)
  File "/usr/lib/python2.7/dist-packages/scipy/interpolate/interpolate.py", line 478, in _check_bounds
    raise ValueError("A value in x_new is above the interpolation "
ValueError: A value in x_new is above the interpolation range.

[Nikos Alexandris]

Yet another patch for create_iwave.py

[Nikos Alexandris]

Mea culpa. I tested (manually) all other sensors. Attached a new diff (replaced the 1st patch, please see timeline), tested for

• AVNIR
• PRISM-B_new, PRISM-F_new, PRISM-N_new -- see attached files, maybe headers not correct!
• VGT1, VGT2
• geoeye1
• ikonos
• landsat
• pleiades1a, pleiades1b
• quickbird2
• rapideye
• spot6, spot7
• worldview2, worldview3

Not sure about the files spot.2, spot.3, spot.4 and spot.5. They should be merged I guess in one file for SPOT5's spectral responses.

Please, keep commented print statements as introduced by the last diff. I plan to rewrite this helper script, and, try adding appropriate automatic testing.

[neteler]

Sorry, I'm lost with the attachments. Yesterday I already fixed PRISM*.csv in r69384 (backports: r69385, r69386).

Please *replace* patches if a new version is attached rather than accumulating them.

Attached a new diff (replaced the 1st patch, please see timeline)

... it is unclear to me what to pick now.

[neteler]

Using the latest patch I still get this error

python create_iwave.py sensors_csv/quickbird2.csv 
 > Getting sensor name from csv file: quickbird2
 > Number of bands found: 5
Traceback (most recent call last):
  File "create_iwave.py", line 281, in <module>
    main()
  File "create_iwave.py", line 267, in main
    write_cpp(bands, values, sensor, os.path.dirname(inputfile))
  File "create_iwave.py", line 185, in write_cpp
    fi, li = interpolate_band(values[:,[0,b+1]])
  File "create_iwave.py", line 119, in interpolate_band
    filter_f = f(np.arange(start, stop, step))
  File "/usr/lib64/python2.7/site-packages/scipy/interpolate/polyint.py", line 79, in __call__
    y = self._evaluate(x)
  File "/usr/lib64/python2.7/site-packages/scipy/interpolate/interpolate.py", line 497, in _evaluate
    out_of_bounds = self._check_bounds(x_new)
  File "/usr/lib64/python2.7/site-packages/scipy/interpolate/interpolate.py", line 527, in _check_bounds
    raise ValueError("A value in x_new is above the interpolation "
ValueError: A value in x_new is above the interpolation range.

[Nikos Alexandris]

We need another tester to nail this one down. It works for *all* sensors here. Below testing for QB2 and L8:

→ rm create_iwave.py 
rm: remove regular file ‘create_iwave.py’? y

→ svn up
Updating '.':
Restored 'create_iwave.py'
At revision 69402.

→ patch -p0 -i create_iwave_patch.diff 
patching file create_iwave.py

→ python create_iwave.py sensors_csv/quickbird2.csv 

 > Getting sensor name from csv file: quickbird2
 > Number of bands found: 5
 > Filter functions exported to quickbird2_cpp_template.txt
 > Please check file for possible errors before inserting into IWave.cpp
 > Don't forget to add the necessary data to IWave.h

→ python create_iwave.py sensors_csv/landsat_8.csv 

 > Getting sensor name from csv file: landsat_8
 > Number of bands found: 9
 > Filter functions exported to landsat_8_cpp_template.txt
 > Please check file for possible errors before inserting into IWave.cpp
 > Don't forget to add the necessary data to IWave.h

[neteler]

ok cleanup helped :)

cd imagery/i.atcorr/
rm -f create_iwave.py
svn up
wget https://trac.osgeo.org/grass/raw-attachment/ticket/3027/create_iwave_patch.diff
dos2unix create_iwave_patch.diff
patch -p0 -i create_iwave_patch.diff

# test
python create_iwave.py sensors_csv/landsat_8.csv

 > Getting sensor name from csv file: landsat_8
 > Number of bands found: 9
 > Filter functions exported to landsat_8_cpp_template.txt
 > Please check file for possible errors before inserting into IWave.cpp
 > Don't forget to add the necessary data to IWave.h

# perfect.

# run on all files:
for i in `ls  sensors_csv/*.csv` ; do python create_iwave.py $i ; done

Perhaps all sensors should be regenerated and compared to the current content in iwave.cpp/iwave.h.

Wish: also generate the HTML snippet and iwave.h code as a by-product.

[neteler]

In 69411:

i.atcorr create_iwave.py: major cleanup; fixes #3027 (contributed by nikosa, Nikos Alexandris)

[neteler]

Keeping open for wish and backport.

[mlennert]

I suppose that this can be closed as create_iwave.py was recently fixed in relation to #3482 ?

[neteler]

Replying to mlennert:

I suppose that this can be closed as create_iwave.py was recently fixed in relation to #3482 ?

It was not backported yet (i.atcorr remains outdated in G7.2 and 7.0).

Hence moving target to 7.2 for bundle backport.

[neteler]

Due to too many changes no backport to 7.2 planned: please use GRASS GIS 7.4 where it was fixed.