[GRASS-user] no output when running i.atcorr

Hi all,

I want to do an atmospheric correction of my Landsat ETM+ images using the i.atcorr module in Grass 6.3 on Ubuntu, but I never get an output. When I use a DEM in ialt, the continuous version or the categorized, of the region, I just get “wavelength less than 0.25 micron let’s take s(1)=s(0.25)”. I don’t insert an altitude value in my icdn file however, since I’ve read that that value would overwrite the DEM. When I do put in an altitude value, I get a list with the parameters for the 6S-algorithm, and then ‘percent complete’, so that seems to be ok, but the algorithm stops running, so I don’t get a corrected image. Does anyone have experience with this module, and help me with some tips, preferably so I can use my DEM instead of one altitude value for the whole image?

This is my icdn file when I use the DEM:
8
2 12 08.500 +35.2203028 -2.8893705
1
0
-1
-1000
64

and this when I use one altitude for the whole region:
8
2 12 08.500 +35.2203028 -2.8893705
1
0
-1
1000
-1000
64

Thanks in advance!

Kind regards,
Annekatrien

Annekatrien,

I am doing my first steps with i.atcorr now.
I have added (with support from Yann Chemin) the specs for
IRS 1C LISS which I need for a project. Also did some manual
cosmetics.

To warm up, I am doing first tests with the North Carolina
test data set [1] which contains sample LANDSAT 7 data
at reduced size.

I observe that the module (or the procedure) is pretty slow,
perhaps 5 pixels/second on my 2.5GHz box... said this, some
comments below:

On Fri, Jan 30, 2009 at 6:01 PM, Annekatrien Debien
<annekatrien.debien@gmail.com> wrote:

Hi all,

I want to do an atmospheric correction of my Landsat ETM+ images using the
i.atcorr module in Grass 6.3 on Ubuntu, but I never get an output.

what does that precisely mean? It keeps you at 0% and calculates or
it exists and no output?

With NC data, I get (I hope I got the icnd.txt file right):

# set region:
g.region rast=lsat7_2002_10 -p

# get central scene LatLong coordinates:
g.region -gb | grep ll_c
ll_clon=-78.70014170
ll_clat=35.71092005

cat icnd.txt
8 - geometrical conditions=Landsat ETM+
5 24 15.50 -78.700 35.710 - month day hh.ddd longitude latitude
("hh.ddd" is a decimal hour GMT)
2 - atmospheric mode=midlatitude summer
1 - aerosols model=continental
50 - visibility [km] (aerosol model concentration)
-.600 - target at 600m above sea level
-1000 - sensor on board a satellite
61 - 1th band of ETM+ Landsat 7

To be sure about the time (Landsat overpass, North Carolina is UTC-5):
http://earthobservatory.nasa.gov/MissionControl/overpass.php

# zoom to be faster for a test
d.zoom

# We can check the sun position with r.sunmask:
r.sunmask -s --v lsat7_2002_10 out=dummy year=2002 \
      month=5 day=24 hour=10 minute=30 sec=0 timezone=-5
Using map center coordinates: 636732.750000 217569.000000
Calculating sun position... (using solpos (V. 11 April 2001) from NREL)
2002.05.24, daynum 144, time: 10:30:00 (decimal time: 10.500000)
long: -78.700150, lat: 35.710922, timezone: -5.000000
Solar position: sun azimuth: 117.002434,
sun angle above horz.(refraction corrected): 63.261234
Sunrise time (without refraction): 05:08:34
Sunset time (without refraction): 19:15:46
No map calculation requested. Finished.

# do the correction
i.atcorr lsat7_2002_10 ialt=elevation icnd=icnd.txt \
             oimg=lsat7_2002_10_atcorr

* ****************************** 6s version 4.2b
****************************** *
* geometrical conditions identity
         *
* -------------------------------
         *
* etm+ observation
         *
*
         *
* month: 5 day: 24
        *
* solar zenith angle: 26.71 deg solar azimuthal angle:
117.09 deg *
* view zenith angle: 0.00 deg view azimuthal angle:
0.00 deg *
* scattering angle: 153.29 deg azimuthal angle difference:
117.09 deg *
*
         *
* atmospheric model description
         *
* -----------------------------
         *
* atmospheric model identity :
         *
* midlatitude summer (uh2o=2.93g/cm2,uo3=.319cm-atm)
         *
*
         *
* aerosols type identity :
         *
* Continental aerosols model
         *
*
         *
* optical condition identity :
         *
* visibility : 50.00 km opt. thick. 550nm :
0.1518 *
*
         *
* spectral condition
         *
* ------------------
         *
* etm+ 1
         *
* value of filter function :
         *
* wl inf= 0.435 mic wl sup= 0.520 mic
        *
*
         *
* target type
         *
* -----------
         *
* homogeneous ground
         *
* constant reflectance over the spectra 0.000
        *
*
         *
* target elevation description
         *
* ----------------------------
         *
* ground pressure [mb] 944.86
        *
* ground altitude [km] 0.600
        *
* gaseous content at target level:
         *
* uh2o= 2.189 g/cm2 uo3= 0.317 cm-atm
        *
*
         *
* atmospheric correction activated
         *
* --------------------------------
         *
0%

... unfortunately the sun angle differs from the r.sunmask output.
What's wrong?

When I
use a DEM in ialt, the continuous version or the categorized, of the region,
I just get "wavelength less than 0.25 micron let's take
s(1)=s(0.25)". I don't insert an altitude value in my icdn file however,
since I've read that that value would overwrite the DEM. When I do put in an
altitude value, I get a list with the parameters for the 6S-algorithm, and
then 'percent complete', so that seems to be ok, but the algorithm stops
running, so I don't get a corrected image. Does anyone have experience with
this module, and help me with some tips, preferably so I can use my DEM
instead of one altitude value for the whole image?

I have just tried on a tiny subset of lsat7_2002_10_atcorr but I see
that the module uses the *entire* image:

GRASS 6.5.svn (nc_spm_07):~ > g.region -p
projection: 99 (Lambert Conformal Conic)
...
rows: 18
cols: 26
cells: 468

GRASS 6.5.svn (nc_spm_07):~ > fg
i.atcorr lsat7_2002_10 ialt=elevation icnd=icnd.txt
oimg=lsat7_2002_10_atcorr
D0/0: Computed r0, c289
D0/0: Computed r0, c290
D0/0: Computed r0, c291 <- running out of current region!

So some work is still needed.

Markus

Hi,
i.atcorr has a test_suite in the source code,
it maybe a good start point.
Yann
------------------------
grass/imagery/i.atcorr/test_suite# ls
ETM4_400x400_atms_corr.hdr
ETM4_400x400.hdr
ETM4_atmospheric_input_GRASS.txt
ETM4.res
ETM4_400x400_atms_corr.raw
ETM4_400x400.raw
ETM4_atmospheric_input.txt
--------------------------

2009/2/13 Markus Neteler <neteler@osgeo.org>:

Annekatrien,

I am doing my first steps with i.atcorr now.
I have added (with support from Yann Chemin) the specs for
IRS 1C LISS which I need for a project. Also did some manual
cosmetics.

To warm up, I am doing first tests with the North Carolina
test data set [1] which contains sample LANDSAT 7 data
at reduced size.

I observe that the module (or the procedure) is pretty slow,
perhaps 5 pixels/second on my 2.5GHz box... said this, some
comments below:

On Fri, Jan 30, 2009 at 6:01 PM, Annekatrien Debien
<annekatrien.debien@gmail.com> wrote:

Hi all,

I want to do an atmospheric correction of my Landsat ETM+ images using the
i.atcorr module in Grass 6.3 on Ubuntu, but I never get an output.

what does that precisely mean? It keeps you at 0% and calculates or
it exists and no output?

With NC data, I get (I hope I got the icnd.txt file right):

# set region:
g.region rast=lsat7_2002_10 -p

# get central scene LatLong coordinates:
g.region -gb | grep ll_c
ll_clon=-78.70014170
ll_clat=35.71092005

cat icnd.txt
8 - geometrical conditions=Landsat ETM+
5 24 15.50 -78.700 35.710 - month day hh.ddd longitude latitude
("hh.ddd" is a decimal hour GMT)
2 - atmospheric mode=midlatitude summer
1 - aerosols model=continental
50 - visibility [km] (aerosol model concentration)
-.600 - target at 600m above sea level
-1000 - sensor on board a satellite
61 - 1th band of ETM+ Landsat 7

To be sure about the time (Landsat overpass, North Carolina is UTC-5):
http://earthobservatory.nasa.gov/MissionControl/overpass.php

# zoom to be faster for a test
d.zoom

# We can check the sun position with r.sunmask:
r.sunmask -s --v lsat7_2002_10 out=dummy year=2002 \
     month=5 day=24 hour=10 minute=30 sec=0 timezone=-5
Using map center coordinates: 636732.750000 217569.000000
Calculating sun position... (using solpos (V. 11 April 2001) from NREL)
2002.05.24, daynum 144, time: 10:30:00 (decimal time: 10.500000)
long: -78.700150, lat: 35.710922, timezone: -5.000000
Solar position: sun azimuth: 117.002434,
sun angle above horz.(refraction corrected): 63.261234
Sunrise time (without refraction): 05:08:34
Sunset time (without refraction): 19:15:46
No map calculation requested. Finished.

# do the correction
i.atcorr lsat7_2002_10 ialt=elevation icnd=icnd.txt \
            oimg=lsat7_2002_10_atcorr

* ****************************** 6s version 4.2b
****************************** *
* geometrical conditions identity
        *
* -------------------------------
        *
* etm+ observation
        *
*
        *
* month: 5 day: 24
       *
* solar zenith angle: 26.71 deg solar azimuthal angle:
117.09 deg *
* view zenith angle: 0.00 deg view azimuthal angle:
0.00 deg *
* scattering angle: 153.29 deg azimuthal angle difference:
117.09 deg *
*
        *
* atmospheric model description
        *
* -----------------------------
        *
* atmospheric model identity :
        *
* midlatitude summer (uh2o=2.93g/cm2,uo3=.319cm-atm)
        *
*
        *
* aerosols type identity :
        *
* Continental aerosols model
        *
*
        *
* optical condition identity :
        *
* visibility : 50.00 km opt. thick. 550nm :
0.1518 *
*
        *
* spectral condition
        *
* ------------------
        *
* etm+ 1
        *
* value of filter function :
        *
* wl inf= 0.435 mic wl sup= 0.520 mic
       *
*
        *
* target type
        *
* -----------
        *
* homogeneous ground
        *
* constant reflectance over the spectra 0.000
       *
*
        *
* target elevation description
        *
* ----------------------------
        *
* ground pressure [mb] 944.86
       *
* ground altitude [km] 0.600
       *
* gaseous content at target level:
        *
* uh2o= 2.189 g/cm2 uo3= 0.317 cm-atm
       *
*
        *
* atmospheric correction activated
        *
* --------------------------------
        *
0%

... unfortunately the sun angle differs from the r.sunmask output.
What's wrong?

When I
use a DEM in ialt, the continuous version or the categorized, of the region,
I just get "wavelength less than 0.25 micron let's take
s(1)=s(0.25)". I don't insert an altitude value in my icdn file however,
since I've read that that value would overwrite the DEM. When I do put in an
altitude value, I get a list with the parameters for the 6S-algorithm, and
then 'percent complete', so that seems to be ok, but the algorithm stops
running, so I don't get a corrected image. Does anyone have experience with
this module, and help me with some tips, preferably so I can use my DEM
instead of one altitude value for the whole image?

I have just tried on a tiny subset of lsat7_2002_10_atcorr but I see
that the module uses the *entire* image:

GRASS 6.5.svn (nc_spm_07):~ > g.region -p
projection: 99 (Lambert Conformal Conic)
...
rows: 18
cols: 26
cells: 468

GRASS 6.5.svn (nc_spm_07):~ > fg
i.atcorr lsat7_2002_10 ialt=elevation icnd=icnd.txt
oimg=lsat7_2002_10_atcorr
D0/0: Computed r0, c289
D0/0: Computed r0, c290
D0/0: Computed r0, c291 <- running out of current region!

So some work is still needed.

Markus

--
Yann Chemin
International Center of Water for Food Security
Office: http://www.icwater.org
Perso: http://www.freewebs.com/ychemin
YiKingDo: http://yikingdo.unblog.fr/

Hi Annekatrien,

news regarding i.atcorr!

On Fri, Jan 30, 2009 at 7:01 PM, Annekatrien Debien
<annekatrien.debien@gmail.com> wrote:

Hi all,

I want to do an atmospheric correction of my Landsat ETM+ images using the
i.atcorr module in Grass 6.3 on Ubuntu, but I never get an output. When I
use a DEM in ialt, the continuous version or the categorized, of the region,
I just get "wavelength less than 0.25 micron let's take
s(1)=s(0.25)". I don't insert an altitude value in my icdn file however,
since I've read that that value would overwrite the DEM.

@Yann: do you have any pointer for us?
Since the DEM *is* read in in the code, is the documentation statement true?
I don't understand the code well.

When I do put in an
altitude value, I get a list with the parameters for the 6S-algorithm, and
then 'percent complete', so that seems to be ok, but the algorithm stops
running, so I don't get a corrected image.

As discussed, it just runs "forever", say it takes really long. But there is
hope:

With great help from Yann Chemin, I have submitted to all branches/trunk
some changes:
- cache of 1024 instead of 128
- bugfix in example (channel 1 and 4 was mixed; fixed target height)
- example update (suggesting and using integer DEM)

The use of an integer DEM seems to be the trick.
It is relatively fast now with the North Carolina Landsat (earlier it
didn't finished in several hours and I gave up):

g.region rast=lsat7_2002_40 -p
r.info lsat7_2002_40

# using an integer DEM greatly accelerates the i.atcorr computations
r.mapcalc "elev_int = round(elevation)"

# find mean elevation (target above sea level)
r.univar elev_int

# create control file for channel 4 (NIR)
echo "8 - geometrical conditions=Landsat ETM+
5 24 14.30 -78.691 35.749 - month day hh.ddd longitude latitude
("hh.ddd" is a decimal hour GMT)
2 - atmospheric mode=midlatitude summer
1 - aerosols model=continental
50 - visibility [km] (aerosol model concentration)
-.110 - target at 110m above sea level
-1000 - sensor on board a satellite
64 - 4th band of ETM+ Landsat 7" > icnd.txt

# run atmospheric correction (note the -o for optimization)
i.atcorr -o lsat7_2002_40 ialt=elev_int icnd=icnd.txt oimg=lsat7_2002_40_atcorr

This takes 20 seconds (!) on a "normal" laptop. Without -o it takes
long, after
10 minutes it is still at 0%...

Cheers
Markus