[GRASS-user] Estimating Albedo from Landsat

Hi all,

I'm trying to estimate albedo from landsat5 images but I'm getting
very strange result. Water bodies, for instance, are giving me an
albedo of 0.15 to 0.17, unusualy high.
Here is the procedure I'm doing

1) Convert TM images from grey level to radiance using the Lmax / Lmin
formula (parameters used given bellow [1])
2) Apply 6S atmospheric correction, using standard tropical
atmosphere, continental aerosol mode
3) Calculate shorwave albedo using Liang (2000) formula (based on
bands 1; 3; 4; 5 and 7)

Since the albedo of water areas is comming out too high, I'm guessing
I messed up the i.atcorr procedure. I don't quite understand what
should be the input image (DN or Radiances?) and what is the output
image (reflectance scaled to 255?)

Thanks
Daniel

[1]
Band LMIN LMAX
1 -1,52 193
2 -2,84 365
3 -1,17 264
4 -1,51 221
5 -0,37 30,2
7 -0,15 16,5

Hey Daniel

I think the input map to be given in i.atcorr should not be the DN. Conversion to radiance or reflectance should be done before implementing the 6s algorithm. Have you tried i.landsat.toar to obtain radiance values?

Success!

JM

2010/8/2 Daniel Victoria <daniel.victoria@gmail.com>

Hi all,

I’m trying to estimate albedo from landsat5 images but I’m getting
very strange result. Water bodies, for instance, are giving me an
albedo of 0.15 to 0.17, unusualy high.
Here is the procedure I’m doing

  1. Convert TM images from grey level to radiance using the Lmax / Lmin
    formula (parameters used given bellow [1])
  2. Apply 6S atmospheric correction, using standard tropical
    atmosphere, continental aerosol mode
  3. Calculate shorwave albedo using Liang (2000) formula (based on
    bands 1; 3; 4; 5 and 7)

Since the albedo of water areas is comming out too high, I’m guessing
I messed up the i.atcorr procedure. I don’t quite understand what
should be the input image (DN or Radiances?) and what is the output
image (reflectance scaled to 255?)

Thanks
Daniel

[1]
Band LMIN LMAX
1 -1,52 193
2 -2,84 365
3 -1,17 264
4 -1,51 221
5 -0,37 30,2
7 -0,15 16,5


grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

Hi Jose,

I haven't tryed i.landsat.toar addon but from what I read in the wiki
(http://grass.osgeo.org/wiki/Atmospheric_correction) I could either
use i.landsat.toar or a DN - Radiance conversion formula that is
given...

Will try to recalculate the images, I might have done something wrong

Daniel

2010/8/2 José Miguel Barrios <jmbarriosg@gmail.com>:

Hey Daniel

I think the input map to be given in i.atcorr should not be the DN.
Conversion to radiance or reflectance should be done before implementing the
6s algorithm. Have you tried i.landsat.toar to obtain radiance values?

Success!

JM

2010/8/2 Daniel Victoria <daniel.victoria@gmail.com>

Hi all,

I'm trying to estimate albedo from landsat5 images but I'm getting
very strange result. Water bodies, for instance, are giving me an
albedo of 0.15 to 0.17, unusualy high.
Here is the procedure I'm doing

1) Convert TM images from grey level to radiance using the Lmax / Lmin
formula (parameters used given bellow [1])
2) Apply 6S atmospheric correction, using standard tropical
atmosphere, continental aerosol mode
3) Calculate shorwave albedo using Liang (2000) formula (based on
bands 1; 3; 4; 5 and 7)

Since the albedo of water areas is comming out too high, I'm guessing
I messed up the i.atcorr procedure. I don't quite understand what
should be the input image (DN or Radiances?) and what is the output
image (reflectance scaled to 255?)

Thanks
Daniel

[1]
Band LMIN LMAX
1 -1,52 193
2 -2,84 365
3 -1,17 264
4 -1,51 221
5 -0,37 30,2
7 -0,15 16,5
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

Daniel Victoria:

I haven't tryed i.landsat.toar addon but from what I read in the wiki
(http://grass.osgeo.org/wiki/Atmospheric_correction) I could either
use i.landsat.toar or a DN - Radiance conversion formula that is
given...

Will try to recalculate the images, I might have done something wrong

Hi Daniel!

Did you finally used i.landsat.toar for the Landsat5 - TM bands?

I am trying to use the module on the LT51830332007248MOR00 image(s) and have
several (expected) problems, like:

(note that the metadata file I downloaded via GLOVis has a ".meta" extension)

1. the "LT51830332007248MOR00.meta" isn't read correctly or the module is
looking for a string other than the "date_acquired = 20070905" string which is
the only "date-string" included in the ".meta" file

2. does not recognize the solar elevation value

(the above "cases" are mentioned in the manual, so manual feed is required
here)

3. after manually providing the "solar_elevation" and the "product_date" and
using the "-v" flag, some random (?) segfaults occur.

4. without the "-v" flag it runs but seems to take too much time? OK, the
image set is large ( 7654 cols x 7127 rows x 7 bands). --- Is it really taking
so long or is there something wrong?

Nikos

I have to go out in the field tomorrow, but I have put together a
script, semi-automatically, to do this and calculate NDVI. I can pass
it along for what it is worth.

On Sun, Sep 12, 2010 at 7:35 PM, Nikos Alexandris
<nikos.alexandris@felis.uni-freiburg.de> wrote:

Daniel Victoria:

I haven't tryed i.landsat.toar addon but from what I read in the wiki
(http://grass.osgeo.org/wiki/Atmospheric_correction) I could either
use i.landsat.toar or a DN - Radiance conversion formula that is
given...

Will try to recalculate the images, I might have done something wrong

Hi Daniel!

Did you finally used i.landsat.toar for the Landsat5 - TM bands?

I am trying to use the module on the LT51830332007248MOR00 image(s) and have
several (expected) problems, like:

(note that the metadata file I downloaded via GLOVis has a ".meta" extension)

1. the "LT51830332007248MOR00.meta" isn't read correctly or the module is
looking for a string other than the "date_acquired = 20070905" string which is
the only "date-string" included in the ".meta" file

2. does not recognize the solar elevation value

(the above "cases" are mentioned in the manual, so manual feed is required
here)

3. after manually providing the "solar_elevation" and the "product_date" and
using the "-v" flag, some random (?) segfaults occur.

4. without the "-v" flag it runs but seems to take too much time? OK, the
image set is large ( 7654 cols x 7127 rows x 7 bands). --- Is it really taking
so long or is there something wrong?

Nikos
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

--
Stephen Sefick
____________________________________
| Auburn University |
| Department of Biological Sciences |
| 331 Funchess Hall |
| Auburn, Alabama |
| 36849 |
|___________________________________|
| sas0025@auburn.edu |
| http://www.auburn.edu/~sas0025 |
|___________________________________|

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods. We are mammals, and have not exhausted the
annoying little problems of being mammals.

                            \-K\. Mullis

"A big computer, a complex algorithm and a long time does not equal science."

                          \-Robert Gentleman

Hi Nikos,

I did not use i.landsat.toar. I converted the DN values to radiance
manually, using r.mapcalc and the formula and coefficients provided.

Are you going to use i.atcorr? I'm still trying to figure out which is
the correct procedure for i.atcorr:

1) Convert DN to radiance and run i.atcorr (wiki example)
2) Use DN values in i.atcorr (man page exemple)

And also, what are the output units since they appear to be scaled from 1 to 255

Cheers
Daniel

2010/9/12 Nikos Alexandris <nikos.alexandris@felis.uni-freiburg.de>:

Daniel Victoria:

I haven't tryed i.landsat.toar addon but from what I read in the wiki
(http://grass.osgeo.org/wiki/Atmospheric_correction) I could either
use i.landsat.toar or a DN - Radiance conversion formula that is
given...

Will try to recalculate the images, I might have done something wrong

Hi Daniel!

Did you finally used i.landsat.toar for the Landsat5 - TM bands?

I am trying to use the module on the LT51830332007248MOR00 image(s) and have
several (expected) problems, like:

(note that the metadata file I downloaded via GLOVis has a ".meta" extension)

1. the "LT51830332007248MOR00.meta" isn't read correctly or the module is
looking for a string other than the "date_acquired = 20070905" string which is
the only "date-string" included in the ".meta" file

2. does not recognize the solar elevation value

(the above "cases" are mentioned in the manual, so manual feed is required
here)

3. after manually providing the "solar_elevation" and the "product_date" and
using the "-v" flag, some random (?) segfaults occur.

4. without the "-v" flag it runs but seems to take too much time? OK, the
image set is large ( 7654 cols x 7127 rows x 7 bands). --- Is it really taking
so long or is there something wrong?

Nikos

Those are digital numbers.

On Mon, Sep 13, 2010 at 6:38 AM, Daniel Victoria
<daniel.victoria@gmail.com> wrote:

Hi Nikos,

I did not use i.landsat.toar. I converted the DN values to radiance
manually, using r.mapcalc and the formula and coefficients provided.

Are you going to use i.atcorr? I'm still trying to figure out which is
the correct procedure for i.atcorr:

1) Convert DN to radiance and run i.atcorr (wiki example)
2) Use DN values in i.atcorr (man page exemple)

And also, what are the output units since they appear to be scaled from 1 to 255

Cheers
Daniel

2010/9/12 Nikos Alexandris <nikos.alexandris@felis.uni-freiburg.de>:

Daniel Victoria:

I haven't tryed i.landsat.toar addon but from what I read in the wiki
(http://grass.osgeo.org/wiki/Atmospheric_correction) I could either
use i.landsat.toar or a DN - Radiance conversion formula that is
given...

Will try to recalculate the images, I might have done something wrong

Hi Daniel!

Did you finally used i.landsat.toar for the Landsat5 - TM bands?

I am trying to use the module on the LT51830332007248MOR00 image(s) and have
several (expected) problems, like:

(note that the metadata file I downloaded via GLOVis has a ".meta" extension)

1. the "LT51830332007248MOR00.meta" isn't read correctly or the module is
looking for a string other than the "date_acquired = 20070905" string which is
the only "date-string" included in the ".meta" file

2. does not recognize the solar elevation value

(the above "cases" are mentioned in the manual, so manual feed is required
here)

3. after manually providing the "solar_elevation" and the "product_date" and
using the "-v" flag, some random (?) segfaults occur.

4. without the "-v" flag it runs but seems to take too much time? OK, the
image set is large ( 7654 cols x 7127 rows x 7 bands). --- Is it really taking
so long or is there something wrong?

Nikos

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

--
Stephen Sefick
____________________________________
| Auburn University |
| Department of Biological Sciences |
| 331 Funchess Hall |
| Auburn, Alabama |
| 36849 |
|___________________________________|
| sas0025@auburn.edu |
| http://www.auburn.edu/~sas0025 |
|___________________________________|

Let's not spend our time and resources thinking about things that are
so little or so large that all they really do for us is puff us up and
make us feel like gods. We are mammals, and have not exhausted the
annoying little problems of being mammals.

                            \-K\. Mullis

"A big computer, a complex algorithm and a long time does not equal science."

                          \-Robert Gentleman

Nikos:

Did you finally used i.landsat.toar for the Landsat5 - TM bands?

I haven't tried with landsat 5 data from glovis yet, have some scenes
on order, expect them to arrive in the next few days. It seems to be
working for me with landsat7 data. Band 61 temperatures seem plausible
anyway & metadata file is read cleanly.

I am trying to use the module on the LT51830332007248MOR00
image(s) and have several (expected) problems, like:

(note that the metadata file I downloaded via GLOVis has a
".meta" extension)

1. the "LT51830332007248MOR00.meta" isn't read correctly or
the module is looking for a string other than the "date_acquired =
20070905" string which is the only "date-string" included in the ".meta"
file

2. does not recognize the solar elevation value

(the above "cases" are mentioned in the manual, so manual
feed is required here)

3. after manually providing the "solar_elevation" and the
"product_date" and using the "-v" flag, some random (?) segfaults occur.

I am a bit suspicious of the if(flag->answer) bits of the code, if pointers
and boolean integers are getting mixed up.

I think over the next week or two I'll continue to make small cleanups and
improvements to the i.landsat.toar and i.landsat.acca code.

4. without the "-v" flag it runs but seems to take too much
time? OK, the image set is large ( 7654 cols x 7127 rows x 7 bands). ---
Is it really taking so long or is there something wrong?

how much time? for me it takes a few minutes.

Hamish

Daniel Victoria wrote:

I did not use i.landsat.toar. I converted the DN values to radiance
manually, using r.mapcalc and the formula and coefficients provided.

ok

Are you going to use i.atcorr?

Currently I only require to automatically detect clouds and nothing more than
that. Since I don't require (at the moment) surface (or top-of-canopy)
reflectances, I guess I don't need to run "i.atcorr".

I will use i.landsat.toar (with the default option "method=uncorrected") and
i.landsat.acca afterwards.

I'm still trying to figure out which is the correct procedure for i.atcorr:

1) Convert DN to radiance and run i.atcorr (wiki example)
2) Use DN values in i.atcorr (man page exemple)

IIUC,

the "i.atcorr" module (which implements the 6S algorithm) expects either (toa)
radiance or (toa) reflectance. Option "2" above is not correct.

And also, what are the output units since they appear to be scaled from 1
to 255

I am puzzled as well :frowning: -- I guess it's a (re-)scaling related puzzle? Will
post again when I'll get the answer to this.

FWIW, I am convinced by writers (e.g. Lillesand and Kiefer) supporting that
the term "correction" is not exactly correct in the case of retouching
remotely sensed data with the aim to "remove" atmospheric effects.

Best regards, Nikos

[...]

On Monday 13 of September 2010 04:20:33 stephen sefick wrote:

I have to go out in the field tomorrow, but I have put together a
script, semi-automatically, to do this and calculate NDVI. I can pass
it along for what it is worth.

Sure, and if it's generally helpful we might put that in the wiki (with your
permission of course).

Thanks, Nikos

Nikos:

> Did you finally used i.landsat.toar for the Landsat5 - TM bands?

Hamish:

I haven't tried with landsat 5 data from glovis yet, have some scenes
on order, expect them to arrive in the next few days. It seems to be
working for me with landsat7 data. Band 61 temperatures seem plausible
anyway & metadata file is read cleanly.

And example of metadata (of the image I am using):
<http://edcsns17.cr.usgs.gov/cgi-
bin/EarthExplorer/phtml/results/display_full.phtml?cseq=3119&row=3&fgdc=N&entity_id=LT51830332007248MOR00>.
I guess it's structure is quite different than the "meta" of ETM+ data.

> I am trying to use the module on the LT51830332007248MOR00
> image(s) and have several (expected) problems, like:

> (note that the metadata file I downloaded via GLOVis has a
> ".meta" extension)

> 1. the "LT51830332007248MOR00.meta" isn't read correctly or
> the module is looking for a string other than the "date_acquired =
> 20070905" string which is the only "date-string" included in the ".meta"
> file

> 2. does not recognize the solar elevation value

> (the above "cases" are mentioned in the manual, so manual
> feed is required here)

> 3. after manually providing the "solar_elevation" and the
> "product_date" and using the "-v" flag, some random (?) segfaults occur.

I am a bit suspicious of the if(flag->answer) bits of the code, if pointers
and boolean integers are getting mixed up.

I think over the next week or two I'll continue to make small cleanups and
improvements to the i.landsat.toar and i.landsat.acca code.

> 4. without the "-v" flag it runs but seems to take too much
> time? OK, the image set is large ( 7654 cols x 7127 rows x 7 bands). ---
> Is it really taking so long or is there something wrong?

how much time? for me it takes a few minutes.

Oh man, I guess I am doing something wrong [ again :frowning: ]. It takes too much
time here. I broke the process after >30 min and set up a new (smaller)
computational region (of my interest). I'll give it another try and report
back.

Thanks for all of your work Hamish, Nikos

Hi,

I've just made some small changes to i.landsat.acca and i.landsat.toar in
svn which shouldn't change the behaviour of the module at all, just do
some more error checks and make better use of libgis (writes out more
metadata etc. e.g. added units=Kelvin to band 6 metadata, not sure about
W/m^2 for others or how safe it is to include non-ascii chars in the units
field? (micron "u", etc.)).

there were some troubles if e.g. production date was not set, it would
continue anyway with a very old date. (it would be nice to have a check
that production date is not before acq. date as well, but that's not
in it yet.)
also if the date format was not right it would continue anyway with a
bogus value and crazy value for solar distance, etc.

I've got an old landsat-5 TM scene here from 2004 which seemed to run ok.
i.landsat.toar didn't recognize the header.dat metadata file, so I had
to enter some values by hand.

band 6 temperatures do not look very plausible (30 degC ocean and -6 degC
land) on a sunny autumn day.

i.landsat.acca then does a really bad job finding the clouds, but I think
that's just because the toar.6 band is so badly wrong. ISTR doing the
temperature calc with r.mapcalc from Markus and Helena's book some years
ago, I don't remember it being so badly wrong then. there's some commented
out landsat-5 code in i.landsat.acca though, maybe it is as yet unfinished?

It would be good to compare with i.atcorr, but I haven't figured out the
exact recipe for that module yet.

Probably we should concentrate on the North Carolina dataset for common
tests... its landsat mapset has:

lsat5_1987_10
lsat5_1987_20
lsat5_1987_30
lsat5_1987_40
lsat5_1987_50
lsat5_1987_60
lsat5_1987_70

not enough metadata for toar, might have to redownload the scene from
GloVis:
IMAGE_ID=P016R35_5T871014,PATH=16,ROW=35,DATE=10/14/87,PLATFORM=LANDSAT5

and
lsat7_2000_10
lsat7_2000_20
lsat7_2000_30
lsat7_2000_40
lsat7_2000_50
lsat7_2000_61
lsat7_2000_70
lsat7_2000_80

metadata from hist/ file:
SPACECRAFT_ID=Landsat7,SENSOR_ID=ETM+,ACQUISITION_DATE=2000-03-31,WRS_P
ATH=16,CPF_FILE_NAME=L7CPF20000101_20000331_12,LMAX_BAND1=191.600,LMIN_
BAND1=-6.200,LMAX_BAND2=196.500,LMIN_BAND2=-6.400,LMAX_BAND3=152.900,LM
IN_BAND3=-5.000,LMAX_BAND4=241.100,LMIN_BAND4=-5.100,LMAX_BAND5=31.060,
LMIN_BAND5=-1.000,LMAX_BAND61=17.040,LMIN_BAND61=0.000,LMAX_BAND62=12.6
50,LMIN_BAND62=3.200,LMAX_BAND7=10.800,LMIN_BAND7=-0.350,LMAX_BAND8=243
.100,LMIN_BAND8=-4.700,QCALMAX_BAND1=255.0,QCALMIN_BAND1=1.0,QCALMAX_BA
ND2=255.0,QCALMIN_BAND2=1.0,QCALMAX_BAND3=255.0,QCALMIN_BAND3=1.0,QCALM
AX_BAND4=255.0,QCALMIN_BAND4=1.0,QCALMAX_BAND5=255.0,QCALMIN_BAND5=1.0,
QCALMAX_BAND61=255.0,QCALMIN_BAND61=1.0,QCALMAX_BAND62=255.0,QCALMIN_BA
ND62=1.0,QCALMAX_BAND7=255.0,QCALMIN_BAND7=1.0,QCALMAX_BAND8=255.0,QCAL
MIN_BAND8=1.0,SUN_AZIMUTH=139.6033279,SUN_ELEVATION=51.524652

should we be calculating the sat_zenith= param from SUN_AZIMUTH in the
metadata file, or is that something else? I just used the default 8.2 deg.

> > 1. the "LT51830332007248MOR00.meta" isn't read correctly or
> > the module is looking for a string other than the "date_acquired =
> > 20070905" string which is the only "date-string" included in the
> > ".meta" file

the program wants "yyyy-mm-dd". (the test for that could still be improved)

> > 2. does not recognize the solar elevation value

for me that worked (at least it was reported as given on the command line
when the -v show settings flag was used)

> > 3. after manually providing the "solar_elevation" and the
> > "product_date" and using the "-v" flag, some random (?) segfaults
> > occur.

do you still get those with the new svn version?

> > 4. without the "-v" flag it runs but seems to take too much
> > time? OK, the image set is large ( 7654 cols x 7127 rows x 7 bands).

were any other values bogus? e.g. earth-sun distance not near 1AU?
probably a memory bug if the time variable was in the millions, etc.

Oh man, I guess I am doing something wrong [ again :frowning: ].
It takes too much time here. I broke the process after >30 min
and set up a new (smaller) computational region (of my interest).

AFAIU at least i.landsat.acca ignores the region settings (why?), and
I'm not sure if i.landsat.toar does or not, I suspect it ignores it as
well and works on the full image.

I'll give it another try and report back.

thanks

Thanks for all of your work Hamish, Nikos

really I've just made some minor tests and tidying of the code, the credit
belongs to others. I'm still learning these modules as much as anybody
else ..

Hamish

ps- the grass7 patches will be badly out of date now, as they are a pain
to update for each rev, I'll wait at least until the libgis/libraster
G_fatal_error() and other messages are sync'd with the standardized
versions before cutting a new one of those. or at least until a little
dust settles.

Hamish:

Probably we should concentrate on the North Carolina
dataset for common tests... its landsat mapset has:

...

lsat7_2000_10
lsat7_2000_20
lsat7_2000_30
lsat7_2000_40
lsat7_2000_50
lsat7_2000_61
lsat7_2000_70
lsat7_2000_80

metadata from hist/ file:
SPACECRAFT_ID=Landsat7,SENSOR_ID=ETM+,ACQUISITION_DATE=2000-03-31,WRS_P
ATH=16,CPF_FILE_NAME=L7CPF20000101_20000331_12,LMAX_BAND1=191.600,LMIN_
BAND1=-6.200,LMAX_BAND2=196.500,LMIN_BAND2=-6.400,LMAX_BAND3=152.900,LM
IN_BAND3=-5.000,LMAX_BAND4=241.100,LMIN_BAND4=-5.100,LMAX_BAND5=31.060,
LMIN_BAND5=-1.000,LMAX_BAND61=17.040,LMIN_BAND61=0.000,LMAX_BAND62=12.6
50,LMIN_BAND62=3.200,LMAX_BAND7=10.800,LMIN_BAND7=-0.350,LMAX_BAND8=243
.100,LMIN_BAND8=-4.700,QCALMAX_BAND1=255.0,QCALMIN_BAND1=1.0,QCALMAX_BA
ND2=255.0,QCALMIN_BAND2=1.0,QCALMAX_BAND3=255.0,QCALMIN_BAND3=1.0,QCALM
AX_BAND4=255.0,QCALMIN_BAND4=1.0,QCALMAX_BAND5=255.0,QCALMIN_BAND5=1.0,
QCALMAX_BAND61=255.0,QCALMIN_BAND61=1.0,QCALMAX_BAND62=255.0,QCALMIN_BA
ND62=1.0,QCALMAX_BAND7=255.0,QCALMIN_BAND7=1.0,QCALMAX_BAND8=255.0,QCAL
MIN_BAND8=1.0,SUN_AZIMUTH=139.6033279,SUN_ELEVATION=51.524652

GRASS> i.landsat.toar band=lsat7_2000 sensor=7 date=2000-03-31 product_date=2000-07-02 -v solar_elevation=51.524652 gain=HHHLHLHHL

production date and gains are based on above coeffs & examining the book
values.

GRASS> i.landsat.acca lsat7_2000.toar -2fs out=lsat7_2000.clouds

looks a bit over enthusiastic in the output map.

The stats line looks reasonable-- Cloud cover : 3.332 %

but output mask shows 66% of cells as clouds:
G65> r.stats lsat5_1987.clouds -pl
6 Cold cloud 62.96%
9 Warm cloud 3.07%
* no data 33.97%

?!

(fwiw, I've had very good results using i.landsat.acca with landsat7 over
water, except for high thin cirrus clouds, but it's known those are missed)

Auto-cloud determination for the NC 1987 Landsat5 scene is no good.

Hamish

Hi Hamish!

I was about to report things related to what you describe below but did not do
so when I realised that I don't really have the "correct" metadata file for
the scene I was given. The "meta" file available at Glovis is not really the
complete set of meta-information (is my guess -- missing "sat_zenith=" value
for example).

I ordered the scene again just to get the metadata stuff (since I can't
contact at the moment the person who provided the Landsat5 TM image) and wait
for it.

Hamish wrote:

I've just made some small changes to i.landsat.acca and i.landsat.toar in
svn which shouldn't change the behaviour of the module at all, just do
some more error checks...

The "long and never ending process" was due to missing parameters I guess.
Adding some arbitrary "sat_zenith=" value (actually I used the sun_azimuth
given in the ".meta" file :-)), along with "product_date=2007-09-05
date=2007-09-05 solar_elevation=52.48237006", did the trick and completed the
task within a short time.

FWIW,

the output values (even if the sat_zenith value I used is non-sense) seem to
be as expected for radiance/reflectance values (e.g. band 1,
min=-0.00308741863579021, max=0.39202091888652) .

...and make better use of libgis (writes out more
metadata etc. e.g. added units=Kelvin to band 6 metadata, not sure about
W/m^2 for others or how safe it is to include non-ascii chars in the units
field? (micron "u", etc.)).

Is the following the problem (you expect-ed)? I see strange characters, e.g.:

-------------------
BAND 1 (code 1)
   calibrated digital number (DN): 1.0 to 255.0
   calibration constants (L): -1.520 to 193.000
   at-sensor radiance = 0.76583 � DN + -2.28583
   mean solar exoatmospheric irradiance (ESUN): 1983.000
   at-sensor reflectance = radiance / 492.32067
-------------------
BAND 2 (code 2)
   calibrated digital number (DN): 1.0 to 255.0
   calibration constants (L): -2.840 to 365.000
   at-sensor radiance = 1.44819 � DN + -4.28819
   mean solar exoatmospheric irradiance (ESUN): 1796.000
   at-sensor reflectance = radiance / 445.89406

there were some troubles if e.g. production date was not set, it would
continue anyway with a very old date. (it would be nice to have a check
that production date is not before acq. date as well, but that's not
in it yet.)
also if the date format was not right it would continue anyway with a
bogus value and crazy value for solar distance, etc.

I've got an old landsat-5 TM scene here from 2004 which seemed to run ok.
i.landsat.toar didn't recognize the header.dat metadata file, so I had
to enter some values by hand.

band 6 temperatures do not look very plausible (30 degC ocean and -6 degC
land) on a sunny autumn day.

i.landsat.acca then does a really bad job finding the clouds, but I think
that's just because the toar.6 band is so badly wrong.

(
Are the "raw" values given in band 6 in Kelvin or are they the result of some
(re-)scaling to 0-255 grey values?
)

ISTR doing the
temperature calc with r.mapcalc from Markus and Helena's book some years
ago, I don't remember it being so badly wrong then. there's some commented
out landsat-5 code in i.landsat.acca though, maybe it is as yet unfinished?

It would be good to compare with i.atcorr, but I haven't figured out the
exact recipe for that module yet.

Probably we should concentrate on the North Carolina dataset for common
tests... its landsat mapset has:

lsat5_1987_10
lsat5_1987_20
lsat5_1987_30
lsat5_1987_40
lsat5_1987_50
lsat5_1987_60
lsat5_1987_70

not enough metadata for toar, might have to redownload the scene from
GloVis:
IMAGE_ID=P016R35_5T871014,PATH=16,ROW=35,DATE=10/14/87,PLATFORM=LANDSAT5

and
lsat7_2000_10
lsat7_2000_20
lsat7_2000_30
lsat7_2000_40
lsat7_2000_50
lsat7_2000_61
lsat7_2000_70
lsat7_2000_80

metadata from hist/ file:
SPACECRAFT_ID=Landsat7,SENSOR_ID=ETM+,ACQUISITION_DATE=2000-03-31,WRS_P
ATH=16,CPF_FILE_NAME=L7CPF20000101_20000331_12,LMAX_BAND1=191.600,LMIN_
BAND1=-6.200,LMAX_BAND2=196.500,LMIN_BAND2=-6.400,LMAX_BAND3=152.900,LM
IN_BAND3=-5.000,LMAX_BAND4=241.100,LMIN_BAND4=-5.100,LMAX_BAND5=31.060,
LMIN_BAND5=-1.000,LMAX_BAND61=17.040,LMIN_BAND61=0.000,LMAX_BAND62=12.6
50,LMIN_BAND62=3.200,LMAX_BAND7=10.800,LMIN_BAND7=-0.350,LMAX_BAND8=243
.100,LMIN_BAND8=-4.700,QCALMAX_BAND1=255.0,QCALMIN_BAND1=1.0,QCALMAX_BA
ND2=255.0,QCALMIN_BAND2=1.0,QCALMAX_BAND3=255.0,QCALMIN_BAND3=1.0,QCALM
AX_BAND4=255.0,QCALMIN_BAND4=1.0,QCALMAX_BAND5=255.0,QCALMIN_BAND5=1.0,
QCALMAX_BAND61=255.0,QCALMIN_BAND61=1.0,QCALMAX_BAND62=255.0,QCALMIN_BA
ND62=1.0,QCALMAX_BAND7=255.0,QCALMIN_BAND7=1.0,QCALMAX_BAND8=255.0,QCAL
MIN_BAND8=1.0,SUN_AZIMUTH=139.6033279,SUN_ELEVATION=51.524652

should we be calculating the sat_zenith= param from SUN_AZIMUTH in the
metadata file, or is that something else? I just used the default 8.2 deg.

Good question! As described above, for the sake of testing I just fed the
"sat_zenith" parameter the "sun_azimuth" value and the module run without
complaints. My guess is that this parameter is important and will make a
difference depending on the date and time of the acquisition.

> > > 1. the "LT51830332007248MOR00.meta" isn't read correctly or
> > > the module is looking for a string other than the "date_acquired =
> > > 20070905" string which is the only "date-string" included in the
> > > ".meta" file

the program wants "yyyy-mm-dd". (the test for that could still be improved)

yep, finally fed it manually

> > > 2. does not recognize the solar elevation value

for me that worked (at least it was reported as given on the command line
when the -v show settings flag was used)

> > > 3. after manually providing the "solar_elevation" and the
> > > "product_date" and using the "-v" flag, some random (?) segfaults
> > > occur.

do you still get those with the new svn version?

In my last test(s) no segfault occured... with or without "-v" !??

> > > 4. without the "-v" flag it runs but seems to take too much
> > > time? OK, the image set is large ( 7654 cols x 7127 rows x 7 bands).

were any other values bogus? e.g. earth-sun distance not near 1AU?

no I don't think so. I can try to repeat it but does it worth the time
spending? I mean, after feeding the module with the "expected" values it runs
fine.

probably a memory bug if the time variable was in the millions, etc.

> Oh man, I guess I am doing something wrong [ again :frowning: ].
> It takes too much time here. I broke the process after >30 min
> and set up a new (smaller) computational region (of my interest).

AFAIU at least i.landsat.acca ignores the region settings (why?), and
I'm not sure if i.landsat.toar does or not, I suspect it ignores it as
well and works on the full image.

confirmed, the "toar" module ignores the region settings. Working on

rows=4906
cols=3707

gives

Rows: 7127
Columns: 7654

> Thanks for all of your work Hamish, Nikos

really I've just made some minor tests and tidying of the code, the credit
belongs to others. I'm still learning these modules as much as anybody
else ..

Well, to me your work (as well as the work of all others) is invaluable. Even
though "nice words" are inexpensive, I still think they are important (when
they are honest). From my side, I just need some timing to contribute back -
all seems to go wrong currently :-p

ps- the grass7 patches will be badly out of date now, as they are a pain
to update for each rev, I'll wait at least until the libgis/libraster
G_fatal_error() and other messages are sync'd with the standardized
versions before cutting a new one of those. or at least until a little
dust settles.

Hamish wrote:

> should we be calculating the sat_zenith= param from
> SUN_AZIMUTH in the metadata file, or is that something
> else? I just used the default 8.2 deg.

Nikos:

Good question! As described above, for the sake of testing
I just fed the "sat_zenith" parameter the "sun_azimuth" value
and the module run without complaints.

hmmm, the sun, earth, and landsat satellite are three different
bodies all moving independently. so time of day of the overpass
(sun angle) is independent of the position of the satellite
(sat_zenith) unless you can calculate the orbit.

Hamish

Nikos wrote:

... I realised that I don't really have the "correct" metadata file for
the scene I was given. The "meta" file available at Glovis is not really
the complete set of meta-information (is my guess -- missing
"sat_zenith=" value for example).

fwiw with Landsat-7 downloads from GloVis I get files like
L7*_MTL.txt which works with the i.landsat.toar metfile= option, it
includes SUN_ELEVATION and SUN_AZIMUTH but nothing obvious about zenith
that I can see.

looking in the source code, landsat_met.c (which parses the metadata file)
isn't looking for satellite zenith- it only comes from the command line,
and in landsat.c it seems that only the DOS2b, DOS3, and DOS4 correction
methods actually use it. So for uncorrected method (which I guess is the
appropriate one if the task is cloud detection??) that variable doesn't
matter.

the output values (even if the sat_zenith value I used is
non-sense) seem to be as expected for radiance/reflectance values

good.

(e.g. band 1,
min=-0.00308741863579021, max=0.39202091888652) .

I think this already got answered a couple of days ago, but is the
reflectance on a scale from 0.0-1.0 (0-100%)?

and what are the W/m^2-like units for radiance used for the output maps?

Is the following the problem (you expect-ed)? I see strange
characters, e.g.:

-------------------
BAND 1 (code 1)
   calibrated digital number (DN): 1.0 to 255.0
   calibration constants (L): -1.520 to 193.000
   at-sensor radiance = 0.76583 � DN + -2.28583
   mean solar exoatmospheric irradiance (ESUN): 1983.000
   at-sensor reflectance = radiance / 492.32067
-------------------
BAND 2 (code 2)
   calibrated digital number (DN): 1.0 to 255.0
   calibration constants (L): -2.840 to 365.000
   at-sensor radiance = 1.44819 � DN + -4.28819
   mean solar exoatmospheric irradiance (ESUN): 1796.000
   at-sensor reflectance = radiance / 445.89406

no, that is something else, it's a non-ASCII char in the fprintf()
string. Maybe it will show up in the email: it's the dot multiplier "·".
for me in "less main.c" it shows up as <B7>, in nedit and vi I see the dot.
I guess that should be %c + the code for it, or just "*" ?
.. displaying it will be terminal and locale dependent ..

In my last test(s) no segfault occured...

good

with or without "-v" !??

both should work.

Even though "nice words" are inexpensive,

I don't think enough people really appreciate that cost/benefit ..

I still think they are important (when they are honest). From my side,
I just need some timing to contribute back -
all seems to go wrong currently :-p

best of luck & enjoy your time,
Hamish

Hamish:

> > should we be calculating the sat_zenith= param from
> > SUN_AZIMUTH in the metadata file, or is that something
> > else? I just used the default 8.2 deg.

Nikos:

> Good question! As described above, for the sake of testing
> I just fed the "sat_zenith" parameter the "sun_azimuth" value
> and the module run without complaints.

Hamish:

hmmm, the sun, earth, and landsat satellite are three different
bodies all moving independently. so time of day of the overpass
(sun angle) is independent of the position of the satellite
(sat_zenith) unless you can calculate the orbit.

Right. It was only for testing to see if the module runs at all.

Hamish:

[...]

band 6 temperatures do not look very plausible (30 degC ocean and -6 degC
land) on a sunny autumn day.

here

<http://landsathandbook.gsfc.nasa.gov/handbook/handbook_htmls/chapter11/chapter11.html&gt;

something relevant I guess ("11.3.3 Band 6 Conversion to Temperature")

Nikos

[...]

Nikos wrote:

> ... I realised that I don't really have the "correct" metadata file for
> the scene I was given. The "meta" file available at Glovis is not really
> the complete set of meta-information (is my guess -- missing
> "sat_zenith=" value for example).

Hamish wrote:

fwiw with Landsat-7 downloads from GloVis I get files like
L7*_MTL.txt which works with the i.landsat.toar metfile= option, it
includes SUN_ELEVATION and SUN_AZIMUTH but nothing obvious about zenith
that I can see.

Check that out (got some metadata files finally -- no ".met" file though) --
just a part of one of the ".lea" files, i.e. the lea_01.txt:

--%<---

------- Map proj. ancillary rec. number 1 -------
record sequence number = 3
1-st record subtype code = 36
record type code = 36
2-nd subtype code = 18
3-rd subtype code = 9
length of this record = 4320
input line nom. number of scene pixel = 7452
input image nom. number of scene lines = 6884
nominal scale of input inter-pixel... = 30.0000000
nominal scale of input inter-line... = 30.0000000
image skew at scene centre = 2.9918563
UTM datum and zone number for input... = GRS80 34
nom. WRS northing of centre in metres = 0.0000000
nom. WRS easting of centre in metres = 0.0000000
northing of input image centre in metres= 4309451.5000000
easting of input image centre in metres = 723735.5625000
vert. offset of scene centre to WRS... = 0.0000000
horiz. offset of scene centre to WRS... = 0.0000000
orient. of input image centre in degrees= -8.8923416
blanks =
num. of pixels per line of proc. image = 6336
num. of lines per line of proc. image = 6884
scale of proc. inter-pix. distance... = 30.0000000
scale of proc. inter-line distance... = 30.0000000
UTM zone number for processed image = 34
line number in proc. image at WRS... = 3442.0000000
pixel number in proc. image at WRS... = 3726.0000000
orientation of proc. image centre... = -8.8923416
nomimal satellite orbital inclination = 98.1999969
nomimal ascending node = 103.0374985
nomimal altitude in metres = 712278.3125000
nom. ground speed in metres per second = 6746.4267578
satellite heading in degrees = 10.5169640
spare =
cross-track field of view in degrees = 15.3900003
sensor scan rate in scans per second = 13.9934511
sensor active sampli. rate in samples...= 104047.4453125
sun elevation angle at WRS centre... = 52.4986238
sun azimuth angle at WRS centre... = 142.9811328
top left latitude = 39.8559723
top left longitude = 22.3088074
top right latitude = 39.7971458
top right longitude = 24.9181290
bottom left latitude = 37.9956665
bottom left longitude = 22.2750587
bottom right latitude = 37.9405975
bottom right longitude = 24.8173637
first ten zero fill =

------- Radiom. ancillary rec. number 1 (total 1) -------
record sequence number = 4
1-st record subtype code = 63
record type code = 36
2-nd subtype code = 18
3-rd subtype code = 9
length of this record = 4320
band number = 1
lower radiance limit = -15
upper radiance limit = 1521
blanks =
offset coefficient = -1.5000000000e+00
gain coefficient = 6.0235294118e-01

[...]

radiance to reflectance conv. factor = 2.0569763e-03
--%<---

Is "satellite heading in degrees = 10.5169640" the "sat_zenith" value ?

looking in the source code, landsat_met.c (which parses the metadata file)
isn't looking for satellite zenith- it only comes from the command line,
and in landsat.c it seems that only the DOS2b, DOS3, and DOS4 correction
methods actually use it. So for uncorrected method (which I guess is the
appropriate one if the task is cloud detection??) ...

I think this is the appropriate way to go for "acca"

...that variable doesn't matter.

ok.

# so, here it goes... seems to run fine
i.landsat.toar band_prefix=landsat_postfire_east_sterea_ellas
method=uncorrected sensor=5 date=2007-09-05 solar_elevation=52.4986238
product_date=2007-09-05

> the output values (even if the sat_zenith value I used is
> non-sense) seem to be as expected for radiance/reflectance values

good.

> (e.g. band 1,
> min=-0.00308741863579021, max=0.39202091888652) .

I think this already got answered a couple of days ago, but is the
reflectance on a scale from 0.0-1.0 (0-100%)?

and what are the W/m^2-like units for radiance used for the output maps?

# hmmm...
r.info landsat_postfire_east_sterea_ellas.toar.6 -r
min=203.36601783578
max=319.083636050686

No way the min value is true ( 203 K = -70.15 C -- irrational for Greece - for
the moment! -)

The max value (319 K = 45.85 C) seems to fit for Greece but not for September!
Maybe for July/August.

> Is the following the problem (you expect-ed)? I see strange
> characters, e.g.:
>
> -------------------
>
> BAND 1 (code 1)
>
> calibrated digital number (DN): 1.0 to 255.0
> calibration constants (L): -1.520 to 193.000
> at-sensor radiance = 0.76583 � DN + -2.28583
> mean solar exoatmospheric irradiance (ESUN): 1983.000
> at-sensor reflectance = radiance / 492.32067
>
> -------------------
>
> BAND 2 (code 2)
>
> calibrated digital number (DN): 1.0 to 255.0
> calibration constants (L): -2.840 to 365.000
> at-sensor radiance = 1.44819 � DN + -4.28819
> mean solar exoatmospheric irradiance (ESUN): 1796.000
> at-sensor reflectance = radiance / 445.89406

no, that is something else, it's a non-ASCII char in the fprintf()
string. Maybe it will show up in the email: it's the dot multiplier "·".

Yep, it shows in the e-mail and I now understand that I was to fast "reading"
those lines before.

for me in "less main.c" it shows up as <B7>, in nedit and vi I see the dot.
I guess that should be %c + the code for it, or just "*" ?

or simply "x" ?

.. displaying it will be terminal and locale dependent ..

Nikos:

# so, here it goes... seems to run fine
i.landsat.toar band_prefix=landsat_postfire_east_sterea_ellas
method=uncorrected sensor=5 date=2007-09-05 solar_elevation=52.4986238
product_date=2007-09-05

...

r.info landsat_postfire_east_sterea_ellas.toar.6 -r
min=203.36601783578
max=319.083636050686

No way the min value is true ( 203 K = -70.15 C -- irrational for Greece -
for the moment! -)

The max value (319 K = 45.85 C) seems to fit for Greece but not for
September! Maybe for July/August.

looking at the image/histogram again those min and max are outliers I guess.
The main body of information starts at ~280 K (6.85 C) and stops at ~307 K
(33.85 C).

I am not sure how to explain the min < 270 K (those values lie mostly within
water bodies), but the max > 310 K areas are burned surfaces (a few days old
burn scars) which will explain the 45 C (at surface) in a September day in
Greece.

So, those values seem rational after all.

Nikos