[GRASS-dev] Unexpected EVI range from "i.vi"

Nikos Alexandris nik at nikosalexandris.net
Tue Jun 11 15:02:51 PDT 2013 

Markus Neteler wrote:

> if you have a i.landsat.toar corrected scene, please extract the
> values for

Landsat scene: LE71610432005160ASN00

> - water pixel

e.g. WaterCoordinates=765525.446097,2756869.46097

# DNs
r.what  
map=B.Trimmed.1,B.Trimmed.2,B.Trimmed.3,B.Trimmed.4,B.Trimmed.5,B.Trimmed.7,B.Trimmed.8 
coordinates=${WaterCoordinates}

765525.446097|2756869.46097||62|39|30|18|14|12|26

# uncorrected
r.what  
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.2,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,B.Trimmed.ToAR.5,B.Trimmed.ToAR.7,B.Trimmed.ToAR.8 
coordinates=765525.446097,2756869.46097

765525.446097|2756869.46097||0.118040401828469|0.0759433089269006|
0.0508433154533196|0.0384822843294042|0.0232458453463901|
0.0164118682715267|0.0508219851689881

# DOS1
r.what 
map=B.Trimmed.ToCR.DOS1.1,B.Trimmed.ToCR.DOS1.2,B.Trimmed.ToCR.DOS1.3,B.Trimmed.ToCR.DOS1.4,B.Trimmed.ToCR.DOS1.5,B.Trimmed.ToCR.DOS1.7,B.Trimmed.ToCR.DOS1.8 
coordinates=765525.446097,2756869.46097

765525.446097|2756869.46097||0.0328866250195922|0.0295187172229351|
0.0254622287384144|0.0257491107226925|0.056707690279387|
0.0271800219947907|0.0311700429776789


> - green vegetation

e.g. VegetationCoordinates=732730.388415,2691278.87464

# DNs
r.what  
map=B.Trimmed.1,B.Trimmed.2,B.Trimmed.3,B.Trimmed.4,B.Trimmed.5,B.Trimmed.7,B.Trimmed.8 
coordinates=732730.388415,2691278.87464

732730.388415|2691278.87464||123|132|138|128|35|18|140

# uncorrected
r.what  
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.2,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,B.Trimmed.ToAR.5,B.Trimmed.ToAR.7,B.Trimmed.ToAR.8 
coordinates=732730.388415,2691278.87464

732730.388415|2691278.87464||0.247199226586371|0.291862470817408|
0.282582569007011|0.399097240803081|0.0860689419202952|
0.0335775069404968|0.337889173730889

# DOS1
r.what 
map=B.Trimmed.ToCR.DOS1.1,B.Trimmed.ToCR.DOS1.2,B.Trimmed.ToCR.DOS1.3,B.Trimmed.ToCR.DOS1.4,B.Trimmed.ToCR.DOS1.5,B.Trimmed.ToCR.DOS1.7,B.Trimmed.ToCR.DOS1.8 
coordinates=732730.388415,2691278.87464

732730.388415|2691278.87464||0.188007083485717|0.288838817470502|
0.303782346029874|0.458849655596738|0.132158574576858|
0.0477960483885396|0.37593931432845


> - asphalt

e.g. AsphaltCoordinates=857330.701903,2702202.81571
### quick choice, maybe I didn't pick a good-one...

# DNs
r.what  
map=B.Trimmed.1,B.Trimmed.2,B.Trimmed.3,B.Trimmed.4,B.Trimmed.5,B.Trimmed.7,B.Trimmed.8 
coordinates=857330.701903,2702202.81571

857330.701903|2702202.81571||77|70|77|81|97|67|80

# uncorrected
r.what  
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.2,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,B.Trimmed.ToAR.5,B.Trimmed.ToAR.7,B.Trimmed.ToAR.8 
coordinates=857330.701903,2702202.81571

857330.701903|2702202.81571||0.149800768572215|0.147916362890403|
0.151692805425759|0.245016304855237|0.271546655614681|
0.173763556070419|0.186801179750941

# DOS1
r.what 
map=B.Trimmed.ToCR.DOS1.1,B.Trimmed.ToCR.DOS1.2,B.Trimmed.ToCR.DOS1.3,B.Trimmed.ToCR.DOS1.4,B.Trimmed.ToCR.DOS1.5,B.Trimmed.ToCR.DOS1.7,B.Trimmed.ToCR.DOS1.8 
coordinates=857330.701903,2702202.81571

857330.701903|2702202.81571||0.0710310000522459|0.115958750638791|
0.146583020522661|0.2737976046051|0.354918328217012|0.216160263937489|
0.194481803091202

> ... all bands except for TIR.
> 
> I want then to synthetic channels of one pixel to do
> the i.vi testing:
> 
> # prepare comp. region
> g.region rast=lsat7_2002_10 rows=1 cols=1 -p
g.region rast=B.Trimmed.1 rows=1 cols=1 -p

projection: 1 (UTM)
zone:       39
datum:      wgs84
ellipsoid:  wgs84
north:      2815470
south:      2621340
west:       658560
east:       875880
nsres:      194130
ewres:      217320
rows:       1
cols:       1
cells:      1

> # create synthetic channels (for the letters, we need your values)
> r.mapcalc "water_b7 = aa"
> r.mapcalc "water_b5 = bb"
> r.mapcalc "water_b4 = cc"
> r.mapcalc "water_b3 = dd"
> r.mapcalc "water_b2 = ee"
> r.mapcalc "water_b1 = ff"

Something like...

# for the water pixel
g.region e=765525.446097 n=2756869.46097 rows=1 cols=1 -p
..

for BAND in `g.mlist rast pat=B.Trimmed.ToCR.DOS1.[1234567]`; do r.mapcalc --o 
"Water_${BAND} = `r.what map=${BAND} coordinates=${WaterCoordinates} 
separator=space | cut -d" " -f4`"; done

 
> Likewise for the other "landuses".

# for the vegetation pixel
g.region e=732730.388415 n=2691278.87464 rows=1 cols=1 -p
..

for BAND in `g.mlist rast pat=B.Trimmed.ToCR.DOS1.[1234567]`; do r.mapcalc --o 
"Vegetation_${BAND} = `r.what map=${BAND} coordinates=${VegetationCoordinates} 
separator=space | cut -d" " -f4`"; done

# for the Asphalt pixel
g.region e=857330.701903 n=2702202.81571 rows=1 cols=1 -p
..

for BAND in `g.mlist rast pat=B.Trimmed.ToCR.DOS1.[1234567]`; do r.mapcalc --o 
"Asphalt_${BAND} = `r.what map=${BAND} coordinates=${AsphaltCoordinates} 
separator=space | cut -d" " -f4`"; done


> This allows us to create a tiny
> test environment to check the output. I am just lacking a prepared
> Landsat scene (and the time to process it). Since you already have it,
> you may extract the values easily with r.what or the wxGUI.

And a small bash script to loop over all VIs and over all One-Pixel-Samples 
[*].  Finally, the results are (min == max in this "rows=1 cols=1" case!):

for BAND in `g.mlist rast pat=Water_B.Trimmed.ToCR.DOS1.[1234567]`; do r.info 
-r ${BAND} | grep min; done

raster map(s) available in mapset <LE71610432005160ASN00>:
min=0.0481443750326537
min=0.0490374344458703
min=0.046078533722967
min=0.0454354991260582
min=0.0782650857929503
min=0.0409240395906233

and

g.mlist rast pat=Water*[i2r]raster map(s) available in mapset 
<LE71610432005160ASN00>:
Water_B.Trimmed.ToCR.DOS1.2 ## this line erased below!
Water_arvi
Water_dvi
Water_evi
Water_evi2
Water_gari
Water_gemi
Water_gvi
Water_ipvi
Water_msavi
Water_ndvi
Water_pvi
Water_savi
Water_sr
Water_vari
Water_wdvi

and

for BAND in `g.mlist rast pat=Water*[i2r]`; do r.info -r ${BAND} | grep min; 
done
raster map(s) available in mapset <LE71610432005160ASN00>:

min=-0.00104585171625509
min=-0.0048671942618601
min=-0.0128935704898675 < evi
min=-0.011069571524947 < evi2
min=-0.197757064195647
min=0.191012824457828
min=-0.0196695823724785
min=0.456824639712233
min=1.63041663127611e-322
min=-0.0863507205755333 < ndvi
min=1.30982050129005
min=-0.0131222955034138
min=0.841026072077743
min=0.189696175424249
min=-0.0048671942618601

For vegetation

min=-0.109202427083821
min=0.091932325419388
min=0.0906057329745037 < evi
min=0.0890674653350201 < evi2
min=-0.0101650844110309
min=0.0907088977409898
min=-0.0233998614337842
min=0.54750980461334
min=2.27270197086973e-322
min=0.0950196092266791 < ndvi
min=1.88445187922905
min=0.0939677302110911
min=1.20999263673654
min=-0.658097689747127
min=0.091932325419388

and for Asphalt


min=-0.112886677245694
min=0.087995047738715
min=0.0868600935258041 < evi
min=0.0853831642974438 < evi2
min=-0.0115058252232626
min=0.0862210729152634
min=-0.0252559818827325
min=0.54566087158789
min=2.27270197086973e-322
min=0.0913217431757802 < ndvi
min=1.87044523200757
min=0.0901852442168926
min=1.20099907198164
min=-0.652697895300491
min=0.087995047738715


Salutations, Nikos

[*]:

##############################################################################
# loop over all vegetation indices
for VI in \
  arvi \
  dvi \
  evi \
  evi2 \
  gvi \
  gari \
  gemi \
  ipvi \
  msavi \
  ndvi \
  pvi \
  savi \
  sr \
  vari \
  wdvi

do 

  for SAMPLE in \
    Water \
    Vegetation \
    Asphalt
  do

    # set resolution and extent
    g.region rast=${SAMPLE}_B.Trimmed.ToCR.DOS1.1 rows=1 cols=1 -pa

    # derive vegetation indice
    echo "Deriving ${VI}..."
    i.vi viname=${VI} --o \
    blue=${SAMPLE}_B.Trimmed.ToCR.DOS1.1\
    green=${SAMPLE}_B.Trimmed.ToCR.DOS1.2 \
    red=${SAMPLE}_B.Trimmed.ToCR.DOS1.3 \
    nir=${SAMPLE}_B.Trimmed.ToCR.DOS1.4 \
    chan5=${SAMPLE}_B.Trimmed.ToCR.DOS1.5 \
    chan7=${SAMPLE}_B.Trimmed.ToCR.DOS1.7 \
    output=${SAMPLE}_${VI}

  done

done
##############################################################################

More information about the grass-dev mailing list