[GRASS-user] Negative raster values from i.landsat.dehaze

James_Travers · May 2, 2012, 9:10pm

Hi,

I have been trying to get some results from GRASS for a couple of weeks now but am hitting a wall. What I am trying to do is to take an area of land from a series of Landsat 7 ETM+ images and calculate NDVI over the area for the purpose of comparison over a time-period. This obviously requires atmospheric correction, so the command sequence I am using after importing the images for a given mapset is:

Generate radiance values

i.landsat.toar input_prefix=“B” output_prefix=“R” metfile=“./L71174038_03820080324/L71174038_03820080324_MTL.txt” sensor=“tm7”

Atmospheric correction

i.tasscap -7 band1=R1 band2=R2 band3=R3 band4=R4 band5=R5 band7=R7 outprefix=TC
i.landsat.dehaze band1=R1 band2=R2 band3=R3 band4=R4 band5=R5 band7=R7 tasscap4=TC.4 outprefix=D

Generate NDVI map

r.mapcalc ‘ndvi=if((D.3 == 0 && D.4 == 0), -1, (D.4 - D.3) / (D.4 + D.3))’

Then I mask the NDVI map with an area raster which I already created. There is probably a better way of doing this.

Mask NDVI map with individual regions

r.mapcalc ‘area_ndvi = ndvi * area_mask@PERMANENT’

Finally I want an average value of the NDVI over the area. There is definitely a better way of doing this but I couldn’t find a grass command that spits out a single value. Instead I used r.stats to give the area of each group of values and then wrote a Perl script to average these.

Process statistics

r.stats -a --q area_ndvi | ./stats.pl

Ok, so here come the questions:

Is there a better way to simply get an average of all values over a given area to replace my stats.pl script?
and this is the big one. When I perform i.landsat.dehaze I get rasters with primarily negative values. This renders the NDVI calculation meaningless as they are ratios of +ve and -ve values. I am presuming this is a mistake as negative reflectance has little meaning anyway - did I use the wrong sequence of commands? Below are info outputs for the red-band after the TOAR and dehaze operations respectively:

Many thanks for any help on this,
Jamie

Timestamp: none

Type of Map: raster Number of Categories: 255
Data Type: DCELL
Rows: 6961
Columns: 7911
Total Cells: 55068471
Projection: UTM (zone 36)
N: 3619200 S: 3410370 Res: 30
E: 843630 W: 606300 Res: 30
Range of data: min = -0.0129146013370688 max = 0.605436510681785

Data Description:
generated by i.landsat.toar

Comments:
Reflectance of Landsat-7 ETM+ (method uncorrected)
----------------------------------------------------------------
Acquisition date … 2008-03-24
Production date … 2012-02-23

Earth-sun distance (d) … 0.99699779
Digital number (DN) range … 1 to 255
Calibration constants (Lmin to Lmax) … -5.000 to +234.400
DN to Radiance (gain and bias) … +0.94252 and -5.94252
Mean solar irradiance (ESUN) … 1551.000
Reflectance = Radiance divided by … 387.15868
-----------------------------------------------------------------

i.landsat.toar input_prefix=“B” output_prefix=“R” metfile="./L711740\
38_03820080324/L71174038_03820080324_MTL.txt" sensor=“tm7” method="u\
ncorrected" percent=0.01 pixel=1000 sat_zenith=8.2000 rayleigh=0.0

±---------------------------------------------------------------------------+
±---------------------------------------------------------------------------+
Layer: D.3 Date: Wed May 2 20:22:21 2012
Mapset: L71174038_03820080324 Login of Creator: jamie
Location: area
DataBase: /home/jamie/grassdata
Title: ( D.3 )
Timestamp: none
----------------------------------------------------------------------------

Type of Map: raster Number of Categories: 255
Data Type: DCELL
Rows: 334
Columns: 166
Total Cells: 55444
Projection: UTM (zone 36)
N: 3493000 S: 3483000 Res: 29.94011976
E: 649000 W: 644000 Res: 30.12048193
Range of data: min = -0.463376848126162 max = 0.19725822381175

Data Description:
generated by r.mapcalc

Comments:
R3 - (TC.4 - 0.134975) * -3.08586

±---------------------------------------------------------------------------+

James_Travers · May 3, 2012, 8:05pm

Ok, I have done a bit more digging and found r.univar for the stats problem - don’t know how I missed that one!

On the atmospheric correction I have ditched i.landsat.dehaze in favour of the presumed higher accuracy of i.atcorr. So, with this I have created a control file with the relevant values for e.g. Band-3:

#--------- start of control file -----
8 - geometrical conditions=Landsat ETM+
3 24 08.020 35.400 32.680 - month day hh.ddd longitude latitude (“hh.ddd” is in GMT decimal hours)
2 - atmospheric mode=midlatitude summer
2 - aerosols model=maritime
50 - visibility [km] (aerosol model concentration)
-0.080 - mean target elevation above sea level [km]
-1000 - sensor on board a satellite
63 - 3rd band of ETM+ Landsat 7
#--------- end of control file -------

When I then run i.atcorr:

$ i.atcorr -r -a iimg=R3 icnd=icnd3.txt oimg=C3D
ERROR: Raster map <dem_float> not found

The error is refers to the default altitude raster-map, but the command help specifies that the ialt argument is optional and I don’t have a DEM for the area. Any ideas what is going on?

Thanks,
Jamie

From: James Travers traversjames@yahoo.ie
To: “grass-user@lists.osgeo.org” grass-user@lists.osgeo.org
Sent: Wednesday, 2 May 2012, 22:10:39
Subject: [GRASS-user] Negative raster values from i.landsat.dehaze

Hi,

I have been trying to get some results from GRASS for a couple of weeks now but am hitting a wall. What I am trying to do is to take an area of land from a series of Landsat 7 ETM+ images and calculate NDVI over the area for the purpose of comparison over a time-period. This obviously requires atmospheric correction, so the command sequence I am using after importing the images for a given mapset is:

Generate radiance values

i.landsat.toar input_prefix=“B” output_prefix=“R” metfile=“./L71174038_03820080324/L71174038_03820080324_MTL.txt” sensor=“tm7”

Atmospheric correction

i.tasscap -7 band1=R1 band2=R2 band3=R3 band4=R4 band5=R5 band7=R7 outprefix=TC
i.landsat.dehaze band1=R1 band2=R2 band3=R3 band4=R4 band5=R5 band7=R7 tasscap4=TC.4 outprefix=D

Generate NDVI map

r.mapcalc ‘ndvi=if((D.3 == 0 && D.4 == 0), -1, (D.4 - D.3) / (D.4 + D.3))’

Then I mask the NDVI map with an area raster which I already created. There is probably a better way of doing this.

Mask NDVI map with individual regions

r.mapcalc ‘area_ndvi = ndvi * area_mask@PERMANENT’

Finally I want an average value of the NDVI over the area. There is definitely a better way of doing this but I couldn’t find a grass command that spits out a single value. Instead I used r.stats to give the area of each group of values and then wrote a Perl script to average these.

Process statistics

r.stats -a --q area_ndvi | ./stats.pl

Ok, so here come the questions:

Is there a better way to simply get an average of all values over a given area to replace my stats.pl script?
and this is the big one. When I perform i.landsat.dehaze I get rasters with primarily negative values. This renders the NDVI calculation meaningless as they are ratios of +ve and -ve values. I am presuming this is a mistake as negative reflectance has little meaning anyway - did I use the wrong sequence of commands? Below are info outputs for the red-band after the TOAR and dehaze operations respectively:

Many thanks for any help on this,
Jamie

Timestamp: none

Type of Map: raster Number of Categories: 255
Data Type: DCELL
Rows: 6961
Columns: 7911
Total Cells: 55068471
Projection: UTM (zone 36)
N: 3619200 S: 3410370 Res: 30
E: 843630 W: 606300 Res: 30
Range of data: min = -0.0129146013370688 max = 0.605436510681785

Data Description:
generated by i.landsat.toar

Comments:
Reflectance of Landsat-7 ETM+ (method uncorrected)
----------------------------------------------------------------
Acquisition date … 2008-03-24
Production date … 2012-02-23

Earth-sun distance (d) … 0.99699779
Digital number (DN) range … 1 to 255
Calibration constants (Lmin to Lmax) … -5.000 to +234.400
DN to Radiance (gain and bias) … +0.94252 and -5.94252
Mean solar irradiance (ESUN) … 1551.000
Reflectance = Radiance divided by … 387.15868
-----------------------------------------------------------------

i.landsat.toar input_prefix=“B” output_prefix=“R” metfile="./L711740\
38_03820080324/L71174038_03820080324_MTL.txt" sensor=“tm7” method="u\
ncorrected" percent=0.01 pixel=1000 sat_zenith=8.2000 rayleigh=0.0

±---------------------------------------------------------------------------+
±---------------------------------------------------------------------------+
Layer: D.3 Date: Wed May 2 20:22:21 2012
Mapset: L71174038_03820080324 Login of Creator: jamie
Location: area
DataBase: /home/jamie/grassdata
Title: ( D.3 )
Timestamp: none
----------------------------------------------------------------------------

Type of Map: raster Number of Categories: 255
Data Type: DCELL
Rows: 334
Columns: 166
Total Cells: 55444
Projection: UTM (zone 36)
N: 3493000 S: 3483000 Res: 29.94011976
E: 649000 W: 644000 Res: 30.12048193
Range of data: min = -0.463376848126162 max = 0.19725822381175

Data Description:
generated by r.mapcalc

Comments:
R3 - (TC.4 - 0.134975) * -3.08586

±---------------------------------------------------------------------------+

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

neteler · May 3, 2012, 11:01pm

Hi James,

On Thu, May 3, 2012 at 10:05 PM, James Travers <traversjames@yahoo.ie> wrote:

Ok, I have done a bit more digging and found r.univar for the stats problem
- don't know how I missed that one!

glad you found it. If you see any potential to improve the documentation,
please make suggestions.

On the atmospheric correction I have ditched i.landsat.dehaze in favour of
the presumed higher accuracy of i.atcorr.

i.landsat.dehaze is "just" a simple statistical approach.

So, with this I have created a
control file with the relevant values for e.g. Band-3:

#--------- start of control file -----
8 - geometrical conditions=Landsat ETM+
3 24 08.020 35.400 32.680 - month day hh.ddd longitude latitude ("hh.ddd" is in GMT decimal hours)
2 - atmospheric mode=midlatitude summer
2 - aerosols model=maritime
50 - visibility [km] (aerosol model concentration)
-0.080 - mean target elevation above sea level [km]
-1000 - sensor on board a satellite
63 - 3rd band of ETM+ Landsat 7
#--------- end of control file -------
When I then run i.atcorr:

$ i.atcorr -r -a iimg=R3 icnd=icnd3.txt oimg=C3D
ERROR: Raster map <dem_float> not found

I just tried with the example of the manual page in the North Carolina
sample dataset and it worked ok. But: don't put the comment lines
into the control file! I fell into this trap in the beginning and i.atcorr
offered weird error messages.

The error is refers to the default altitude raster-map, but the command help
specifies that the ialt argument is optional and I don't have a DEM for the
area. Any ideas what is going on?

I suspect that you have the two comment lines in the file which
unfortunately confuse i.atcorr.

Markus

Hamish · May 4, 2012, 10:51am

James wrote:

> Ok, I have done a bit more digging and found r.univar
> for the stats problem - don't know how I missed that one!

MarkusN:

glad you found it. If you see any potential to improve the
documentation, please make suggestions.

historical footnote: I wrote r.univar(.c) because I didn't notice that
the earlier r.univar(.sh) existed until I was near done. I guess the
same thinness in the docs still exists. Happy mistake.

Markus:

I just tried with the example of the manual page in the
North Carolina sample dataset and it worked ok. But: don't
put the comment lines into the control file! I fell into this
trap in the beginning and i.atcorr
offered weird error messages.

...

I suspect that you have the two comment lines in the file
which unfortunately confuse i.atcorr.

imagery/i.atcorr/GeomCond.cpp's GeomCond::parse():
cin >> igeom;
cin.ignore(numeric_limits<int>::max(),'\n'); /* read the rest of the scraps, like comments */
...

can a C++'er recommend ways to silently skip empty lines and lines
beginning with '#'?

thanks,
Hamish

James_Travers · May 4, 2012, 11:16am

Hi Hamish & Markus,

Thanks for the quick replies. I did have a look around the documentation for the functionality provided by r.univar, but I think the fact that I couldn’t see it in any of the relevant examples put me off track.

The comment solution makes sense, and I don’t think it is a file that needs extensive commenting, but coming from a UNIX background the hash marks screamed comment at me - maybe angle brackets would be a better choice in the documentation:

<--------- start of control file ----->

I will try again later with them removed anyway.

Cheers,
Jamie

From: Hamish hamish_b@yahoo.com
To: James Travers traversjames@yahoo.ie; Markus Neteler neteler@osgeo.org
Cc: GRASS user list grass-user@lists.osgeo.org
Sent: Friday, 4 May 2012, 11:51:44
Subject: Re: [GRASS-user] Negative raster values from i.landsat.dehaze

James wrote:

Ok, I have done a bit more digging and found r.univar
for the stats problem - don’t know how I missed that one!
MarkusN:
glad you found it. If you see any potential to improve the
documentation, please make suggestions.

historical footnote: I wrote r.univar(.c) because I didn’t notice that
the earlier r.univar(.sh) existed until I was near done. I guess the
same thinness in the docs still exists. Happy mistake.

Markus:

I just tried with the example of the manual page in the
North Carolina sample dataset and it worked ok. But: don’t
put the comment lines into the control file! I fell into this
trap in the beginning and i.atcorr
offered weird error messages.
…
I suspect that you have the two comment lines in the file
which unfortunately confuse i.atcorr.

imagery/i.atcorr/GeomCond.cpp’s GeomCond::parse():
cin >> igeom;
cin.ignore(numeric_limits::max(),‘\n’); /* read the rest of the scraps, like comments */
…

can a C++'er recommend ways to silently skip empty lines and lines
beginning with ‘#’?

thanks,
Hamish

James_Travers · May 5, 2012, 6:42pm

Hi again,

I stripped the comments from the icnd file for the i.atcorr command, but still get the same error:

$ i.atcorr -r -a iimg=R3 icnd=icnd3.txt oimg=C3
ERROR: Raster map <dem_float> not found

I tried then to generate a DEM by just applying a fixed 80m value to a raster - I have no idea if this creates a valid elevation map, by the way. This allows i.atcorr to run but the resulting map has all zeros in the non-null areas:

$ r.mapcalc ‘dem=80’
100%
$ i.atcorr -r -a iimg=R3 ialt=dem icnd=icnd3.txt oimg=C3
100%
$ r.univar C3
100%
total null and non-null cells: 55444
total null cells: 17574

Of the non-null cells:

n: 37870
minimum: 0
maximum: 0
range: 0
mean: 0
mean of absolute values: 0
standard deviation: 0
variance: 0
variation coefficient: -nan %
sum: 0

Any ideas what is wrong here or how I can get i.atcorr to run without a valid elevation map?

Thanks,
Jamie

From: James Travers traversjames@yahoo.ie
To: Hamish hamish_b@yahoo.com; Markus Neteler neteler@osgeo.org
Cc: GRASS user list grass-user@lists.osgeo.org
Sent: Friday, 4 May 2012, 12:16:45
Subject: Re: [GRASS-user] Negative raster values from i.landsat.dehaze

Hi Hamish & Markus,

Thanks for the quick replies. I did have a look around the documentation for the functionality provided by r.univar, but I think the fact that I couldn’t see it in any of the relevant examples put me off track.

The comment solution makes sense, and I don’t think it is a file that needs extensive commenting, but coming from a UNIX background the hash marks screamed comment at me - maybe angle brackets would be a better choice in the documentation:

<--------- start of control file ----->

I will try again later with them removed anyway.

Cheers,
Jamie

From: Hamish hamish_b@yahoo.com
To: James Travers traversjames@yahoo.ie; Markus Neteler neteler@osgeo.org
Cc: GRASS user list grass-user@lists.osgeo.org
Sent: Friday, 4 May 2012, 11:51:44
Subject: Re: [GRASS-user] Negative raster values from i.landsat.dehaze

James wrote:

Ok, I have done a bit more digging and found r.univar
for the stats problem - don’t know how I missed that one!
MarkusN:
glad you found it. If you see any potential to improve the
documentation, please make suggestions.

historical footnote: I wrote r.univar(.c) because I didn’t notice that
the earlier r.univar(.sh) existed until I was near done. I guess the
same thinness in the docs still exists. Happy mistake.

Markus:

I just tried with the example of the manual page in the
North Carolina sample dataset and it worked ok. But: don’t
put the comment lines into the control file! I fell into this
trap in the beginning and i.atcorr
offered weird error messages.
…
I suspect that you have the two comment lines in the file
which unfortunately confuse i.atcorr.

imagery/i.atcorr/GeomCond.cpp’s GeomCond::parse():
cin >> igeom;
cin.ignore(numeric_limits::max(),‘\n’); /* read the rest of the scraps, like comments */
…

can a C++'er recommend ways to silently skip empty lines and lines
beginning with ‘#’?

thanks,
Hamish

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