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
±---------------------------------------------------------------------------+
| Layer: R3 Date: Wed May 2 20:15:30 2012 |
| Mapset: L71174038_03820080324 Login of Creator: jamie |
| Location: area |
| DataBase: /home/jamie/grassdata |
| Title: ( R3 ) |
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 |
±---------------------------------------------------------------------------+ |