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

GRASS GIS GRASS GIS developer
1

Hi devs.

I get an unexpected range of evi values using "i.vi" (G7), i.e.

  # derive EVI
  i.vi viname=evi red=B.Trimmed.ToAR.3 blue=B.Trimmed.ToAR.1
nir=B.Trimmed.ToAR.4 output=evi

  # range is...
  r.info -r evi

  min=-6912.82161611806
  max=2264.42037461018

while evi should range in [-1,1]. Both the traditional ndvi and the new(er)
evi2 are fine, using the same input Landsat bands.

  r.info -r ndvi

  min=-1.8046124457954
  max=0.751867027194004

and

  r.info -r evi2

  min=-0.341919561754256
  max=0.624114796171297

The range of the input L7 bands also looks fine:

  # blue, used for evi
  r.info -r B.Trimmed.ToAR.1

  min=0.085822088033344
  max=0.632558130168047

  # red
  r.info -r B.Trimmed.ToAR.3

  min=-0.00851692799904866
  max=0.640896679014357

  # nir
  r.info -r B.Trimmed.ToAR.4

  min=-0.0207162855708529
  max=0.979352245320124

The region is correctly set to match one of the bands (e.g. nir) and the code
(evi.c) looks to me (the ignorant) clean:

# cat /geo/osgeo/src/grass_trunk/imagery/i.vi/evi.c
..
  tmp = nirchan + 6.0 * redchan - 7.5 * bluechan + 1.0;
..
    result = 2.5 * (nirchan - redchan) / tmp

What am I doing/interpret wrong?

Regards, Nikos

2

Nikos Alexandris wrote:

I get an unexpected range of evi values using "i.vi" (G7), i.e.

  # derive EVI
  i.vi viname=evi red=B.Trimmed.ToAR.3 blue=B.Trimmed.ToAR.1
nir=B.Trimmed.ToAR.4 output=evi

  # range is...
  r.info -r evi

  min=-6912.82161611806
  max=2264.42037461018

Something must be "wrong" with one of the bands... ? Manually calculating the
index gives the same result.

Will check further - sorry for the noise.

Nikos

3

Hi Nikos,

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

- water pixel
- green vegetation
- asphalt

... 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

# 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"

Likewise for the other "landuses". 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.

thanks
Markus

4

Markus Neteler wrote:

Hi Nikos,

Salut!

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

--%<---
I have to report-back, though, that using corrected (DOS1) bands (==
reflectances), EVI looks OK!, i.e.

  min = -1.27297025550878 max = 1.4896100159326
--->%--

- water pixel
- green vegetation
- asphalt

... 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

# 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"

Likewise for the other "landuses". 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.

Noted (for tonight).

With the exception of using DNs (which is not good anyway for VIs), I
wrongly(?) expected numbers ranging in [-1,1]. For example, when using
"i.landsat.toar method=uncorrected" or any other transformed set of values
(meaning radiance or reflectance).

Anyhow, will try to check tonight.

Thanks, Nikos

5

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
##############################################################################

6

Nikos Alexandris wrote:

> I get an unexpected range of evi values using "i.vi" (G7), i.e.

> # derive EVI
> i.vi viname=evi red=B.Trimmed.ToAR.3 blue=B.Trimmed.ToAR.1
> nir=B.Trimmed.ToAR.4 output=evi

> # range is...
> r.info -r evi
  
> min=-6912.82161611806
> max=2264.42037461018

..

Will check further - sorry for the noise.

Back to this post with one example (and, maybe I have to test this exact
pixel/coordinates as in the other post the way Markus N suggested).

In the beginning

g.region rast=B.Trimmed.ToAR.1 -pa res=30 # wrongly map-calced the map to 15m,
thus using res=

Next, calculating the evi for *simple* Top-of-Atmosphere values, that is
"i.landsat.toar method=uncorrected"

  i.vi viname=evi blue=B.Trimmed.ToAR.1 red=B.Trimmed.ToAR.3
nir=B.Trimmed.ToAR.4 output=evi_ToAR --o

gives for example

  r.info -r evi_ToAR
  min=-5905.44171917482
  max=6952.80543731566

Looking for these pixels, e.g.

  r.stats evi_ToAR -g | grep '\-5905'

  784935 2695215 -5905.4417191748
  786465 2694675 -5905.4417191748
  766185 2685615 -5905.4417191748

and

  r.stats evi_ToAR -g | grep '6952\.'

  783435 2690025 6952.8054373157

Then tracing the values that produced this out-of-range value,

  r.what coordinates=784935,2695215
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

  784935|2695215||0.230260364323039|0.110923862670943|0.0614305088322746|
0.35181236673774|-5905.44171917482

and

  r.what coordinates=783435,2690025
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

  783435|2690025||0.228143006540123|0.106632395012542|0.0712654621906476|
0.296610169491526|6952.80543731566

Doing the math in R is like

blue <- 0.230260364323039
red <- 0.110923862670943
nir <- 0.0614305088322746
evi <- ( 2.5 * ( nir - red) ) / (nir + 6 * red - 7 * blue + 1.0)
evi

[1] -1.07453

and

blue <- 0.228143006540123
red <- 0.106632395012542
nir <- 0.0712654621906476
evi <- ( 2.5 * ( nir - red) ) / (nir + 6 * red - 7 * blue + 1.0)
evi

[1] -0.7751909

What am I doing wrong... ? Any idea(s)?

More info

r.univar evi_ToAR
100%

total null and non-null cells: 46883168
total null cells: 11617784

Of the non-null cells:
----------------------
n: 35265384
minimum: -5905.44
maximum: 6952.81
range: 12858.2
mean: -0.041618
mean of absolute values: 0.0923752
standard deviation: 2.31373
variance: 5.35336
variation coefficient: -5559.45 %
sum: -1467674.39355633

Nikos

7

Nikos Alexandris wrote:

> > I get an unexpected range of evi values using "i.vi" (G7), i.e.

[..]

> > # range is...
> > r.info -r evi

> > min=-6912.82161611806
> > max=2264.42037461018

[..]

Looking for these pixels, e.g.

  r.stats evi_ToAR -g | grep '\-5905'

  784935 2695215 -5905.4417191748
  786465 2694675 -5905.4417191748
  766185 2685615 -5905.4417191748

and

  r.stats evi_ToAR -g | grep '6952\.'

  783435 2690025 6952.8054373157

Then tracing the values that produced this out-of-range value,

  r.what coordinates=784935,2695215
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

  784935|2695215||0.230260364323039|0.110923862670943|0.0614305088322746|
0.35181236673774|-5905.44171917482

and

  r.what coordinates=783435,2690025
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

  783435|2690025||0.228143006540123|0.106632395012542|0.0712654621906476|
0.296610169491526|6952.80543731566

A bit more clear perhaps:

# res=30
g.region -pagc rast=B.Trimmed.ToAR.1 res=30

n=2815500
s=2621340
w=658560
e=875880
nsres=30
ewres=30
rows=6472
cols=7244
cells=46883168
center_easting=767220.000000
center_northing=2718420.000000

# what is there?
r.what coordinates=784935,2695215
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

784935|2695215||0.230260364323039|0.110923862670943|0.0614305088322746|
0.35181236673774|-5905.44171917482

# while res=1...
g.region -pagc rast=B.Trimmed.ToAR.1 res=1

n=2815500
s=2621340
w=658560
e=875880
nsres=1
ewres=1
rows=194160
cols=217320
cells=42194851200
center_easting=767220.000000
center_northing=2718420.000000

# ...gives
r.what coordinates=784935,2695215
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

784935|2695215||0.230260364323039|0.110923862670943|0.0614305088322746|
0.35181236673774|-5905.44171917482

Now, "pointing" to the suspect pixel:

# both res=30 & res=1 give the same result, of course
g.region -pagc e=784935 n=2695215 rows=1 cols=1 res=30

n=2695230
s=2621340
w=658560
e=784950
nsres=73890
ewres=126390
rows=1
cols=1
cells=1
center_easting=721755.000000
center_northing=2658285.000000

# what is there?
r.what coordinates=784935,2695215
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

784935|2695215||0.264138088849702|0.372703389833447|
0.438437054236573|-0.203045685279188|0.0970312073830427

[..]

I guess it has to do with the extent/res definitions...

Thanks, Nikos

8

[all deleted]

I did several tests but I don't want to add more "noise" in this thread.

The "problem" boils down to when using the full extent of the map(s), e.g.
an "out-of-range" max value (in one or in a few pixels) appears. When using
single-pixels, all is fine.

Nikos

9

Feeding the test values and the evi(2) formula to r.mapcalc, the
results are more or less in the expected range, still beyond [-1, 1],
but not much. Obviously, there is a bug in i.vi. Furthermore, i.vi is
slower than r.mapcalc. I recommend to replace the C module i.vi with a
i.vi script that calls r.mapcalc.

Markus M

On Sun, Jun 16, 2013 at 8:13 PM, Nikos Alexandris
<nik@nikosalexandris.net> wrote:

Nikos Alexandris wrote:

> > I get an unexpected range of evi values using "i.vi" (G7), i.e.

[..]

> > # range is...
> > r.info -r evi

> > min=-6912.82161611806
> > max=2264.42037461018

[..]

Looking for these pixels, e.g.

      r.stats evi_ToAR -g | grep '\-5905'

      784935 2695215 -5905.4417191748
      786465 2694675 -5905.4417191748
      766185 2685615 -5905.4417191748

and

      r.stats evi_ToAR -g | grep '6952\.'

      783435 2690025 6952.8054373157

Then tracing the values that produced this out-of-range value,

      r.what coordinates=784935,2695215
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

      784935|2695215||0.230260364323039|0.110923862670943|0.0614305088322746|
0.35181236673774|-5905.44171917482

and

      r.what coordinates=783435,2690025
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

      783435|2690025||0.228143006540123|0.106632395012542|0.0712654621906476|
0.296610169491526|6952.80543731566

A bit more clear perhaps:

# res=30
g.region -pagc rast=B.Trimmed.ToAR.1 res=30

n=2815500
s=2621340
w=658560
e=875880
nsres=30
ewres=30
rows=6472
cols=7244
cells=46883168
center_easting=767220.000000
center_northing=2718420.000000

# what is there?
r.what coordinates=784935,2695215
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

784935|2695215||0.230260364323039|0.110923862670943|0.0614305088322746|
0.35181236673774|-5905.44171917482

# while res=1...
g.region -pagc rast=B.Trimmed.ToAR.1 res=1

n=2815500
s=2621340
w=658560
e=875880
nsres=1
ewres=1
rows=194160
cols=217320
cells=42194851200
center_easting=767220.000000
center_northing=2718420.000000

# ...gives
r.what coordinates=784935,2695215
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

784935|2695215||0.230260364323039|0.110923862670943|0.0614305088322746|
0.35181236673774|-5905.44171917482

Now, "pointing" to the suspect pixel:

# both res=30 & res=1 give the same result, of course
g.region -pagc e=784935 n=2695215 rows=1 cols=1 res=30

n=2695230
s=2621340
w=658560
e=784950
nsres=73890
ewres=126390
rows=1
cols=1
cells=1
center_easting=721755.000000
center_northing=2658285.000000

# what is there?
r.what coordinates=784935,2695215
map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR

784935|2695215||0.264138088849702|0.372703389833447|
0.438437054236573|-0.203045685279188|0.0970312073830427

[..]

I guess it has to do with the extent/res definitions...

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

10

Markus Metz wrote:

Feeding the test values and the evi(2) formula to r.mapcalc, the
results are more or less in the expected range, still beyond [-1, 1],
but not much.

[cut]

Yes, but feeding the "suspect" values in r.mapcalc, still gives, correctly,
large/out of range (regarding EVI's expected range) result(s).

Anyhow, just to have a quick-check on "r.what", should I upload the bands in
question somewhere? Would anyone have the time to explain/check why r.what
gives different results depending on the extent/resolution for the same
coordinates? Which, might be expected, but why does "g.region rows=1 cols=1"
set the resolution to... see below:

Nikos

Nikos Alexandris:

[snip]

> Now, "pointing" to the suspect pixel:
>
> # both res=30 & res=1 give the same result, of course
> g.region -pagc e=784935 n=2695215 rows=1 cols=1 res=30

-------------------------------------------^------^-----^^

>
> n=2695230
> s=2621340
> w=658560
> e=784950

----------vvvvvv

> nsres=73890
> ewres=126390

----------^^^^^^

> rows=1
> cols=1
> cells=1
> center_easting=721755.000000
> center_northing=2658285.000000
>
> # what is there?
> r.what coordinates=784935,2695215
> map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR
>
> 784935|2695215||0.264138088849702|0.372703389833447|
> 0.438437054236573|-0.203045685279188|0.0970312073830427
>
> [..]
>
> I guess it has to do with the extent/res definitions...

11

On Mon, June 17, 2013 22:58, Nikos Alexandris wrote:

Markus Metz wrote:

Feeding the test values and the evi(2) formula to r.mapcalc, the
results are more or less in the expected range, still beyond [-1, 1],
but not much.

[cut]

Yes, but feeding the "suspect" values in r.mapcalc, still gives,
correctly,
large/out of range (regarding EVI's expected range) result(s).

Anyhow, just to have a quick-check on "r.what", should I upload the bands
in
question somewhere? Would anyone have the time to explain/check why
r.what
gives different results depending on the extent/resolution for the same
coordinates? Which, might be expected, but why does "g.region rows=1
cols=1"
set the resolution to... see below:

[...]

> g.region -pagc e=784935 n=2695215 rows=1 cols=1 res=30

-------------------------------------------^------^-----^^

>
> n=2695230
> s=2621340
> w=658560
> e=784950

----------vvvvvv

> nsres=73890
> ewres=126390

----------^^^^^^

rows and cols has precedence over res, but res=30 plus -a still has an
effect, so you get an extension that is in multiples of 30, with a
resolution equal to the total extension, since you asked for one row and
one col.

Moritz

12

On Mon, Jun 17, 2013 at 10:58 PM, Nikos Alexandris
<nik@nikosalexandris.net> wrote:

Markus Metz wrote:

Feeding the test values and the evi(2) formula to r.mapcalc, the
results are more or less in the expected range, still beyond [-1, 1],
but not much.

[cut]

Yes, but feeding the "suspect" values in r.mapcalc, still gives, correctly,
large/out of range (regarding EVI's expected range) result(s).

Just to confirm, the suspect values are

blue <- 0.230260364323039
red <- 0.110923862670943
nir <- 0.0614305088322746
evi <- ( 2.5 * ( nir - red) ) / (nir + 6 * red - 7 * blue + 1.0)
evi

[1] -1.07453

and

blue <- 0.228143006540123
red <- 0.106632395012542
nir <- 0.0712654621906476
evi <- ( 2.5 * ( nir - red) ) / (nir + 6 * red - 7 * blue + 1.0)
evi

[1] -0.7751909

?

Using r.mapcalc I get the same results like in R, not
-5905.44171917482 or 6952.80543731566, respectively. Why are the way
out of range values correct if R produces reasonable results?

Markus M

Anyhow, just to have a quick-check on "r.what", should I upload the bands in
question somewhere? Would anyone have the time to explain/check why r.what
gives different results depending on the extent/resolution for the same
coordinates? Which, might be expected, but why does "g.region rows=1 cols=1"
set the resolution to... see below:

That explained Moritz. I would only add that you should align the
region to the raster (not the raster's resolution, these are two
different things), before using r.what, otherwise the raster lib will
do some internal nn resamplng (if the region's resoltion is larger
than the raster's resolution).

Nikos

Nikos Alexandris:

[snip]

> Now, "pointing" to the suspect pixel:
>
> # both res=30 & res=1 give the same result, of course
> g.region -pagc e=784935 n=2695215 rows=1 cols=1 res=30

-------------------------------------------^------^-----^^

>
> n=2695230
> s=2621340
> w=658560
> e=784950

----------vvvvvv

> nsres=73890
> ewres=126390

----------^^^^^^

> rows=1
> cols=1
> cells=1
> center_easting=721755.000000
> center_northing=2658285.000000
>
> # what is there?
> r.what coordinates=784935,2695215
> map=B.Trimmed.ToAR.1,B.Trimmed.ToAR.3,B.Trimmed.ToAR.4,evi_DN,evi_ToAR
>
> 784935|2695215||0.264138088849702|0.372703389833447|
> 0.438437054236573|-0.203045685279188|0.0970312073830427
>
> [..]
>
> I guess it has to do with the extent/res definitions...

13

Moritz Lennert wrote:

>> > g.region -pagc e=784935 n=2695215 rows=1 cols=1 res=30
> --------------------^^^^^^---^^^^^^^------^------^-----^^

>> > n=2695230
>> > s=2621340
>> > w=658560
>> > e=784950

> ---------vvvvv
>> > nsres=73890
>> > ewres=126390
> ---------^^^^^^

rows and cols has precedence over res, but res=30 plus -a still has an
effect, so you get an extension that is in multiples of 30, with a
resolution equal to the total extension, since you asked for one row and
one col.

So, "e=784935" and "n=2695215" are set as requested but the "-a" flag cut's
off a bit while aligning...

Thanks, N

14

Markus Metz wrote:
..

> Anyhow, just to have a quick-check on "r.what", should I upload the bands
> in question somewhere? Would anyone have the time to explain/check why
> r.what gives different results depending on the extent/resolution for the
> same coordinates? Which, might be expected, but why does "g.region
> rows=1 cols=1" set the resolution to... see below:

That explained Moritz. I would only add that you should align the
region to the raster (not the raster's resolution, these are two
different things), before using r.what, otherwise the raster lib will
do some internal nn resamplng (if the region's resoltion is larger
than the raster's resolution).

Right, thanks N

15

Markus M:

>> Feeding the test values and the evi(2) formula to r.mapcalc, the
>> results are more or less in the expected range, still beyond [-1, 1],
>> but not much.
> ------- [cut] -------

O-K, let's get this from the start (if there is still energy...).

First thing: the region!

# with or without "-a", it sets the same extent here
g.region rast=B.Trimmed.ToAR.1 -pag res=30

n=2815500
s=2621340
w=658560
e=875880
nsres=30
ewres=30
rows=6472
cols=7244
cells=46883168

As described in the first post, reflectances(?, -- using "i.landsat.toar
method=uncorrected" in this case for testing) derived from Landsat7 bands that
pertain to the scene LE71610432005160ASN00, range between...

# blue (required for evi)
r.info -r B.Trimmed.ToAR.1

min=0.0756932461701407
max=0.526690453931338

# red
r.info -r B.Trimmed.ToAR.3

min=0.0186573080153069
max=0.53363342702351

# nir
r.info -r B.Trimmed.ToAR.4

min=0.0122557420404096
max=0.815443599640871

# the EVI index
i.vi viname=evi blue=B.Trimmed.ToAR.1 red=B.Trimmed.ToAR.3
nir=B.Trimmed.ToAR.4 output=evi

# ranges between...
r.info -r evi

min=-5905.44171917482
max=6952.80543731566

# locating these high values
r.stats evi -g | grep '\-5905'

784935 2695215 -5905.4417191748
786465 2694675 -5905.4417191748
766185 2685615 -5905.4417191748

# and
r.stats evi -g | grep '6952\.'

783435 2690025 6952.8054373157

# querying for DN, ToAR and respective EVIs in the end
# for example in "coordinates=783435,2690025"
r.what coordinates=783435,2690025
map=B.Trimmed.1,B.Trimmed.ToAR.1,B.Trimmed.3,B.Trimmed.ToAR.3,B.Trimmed.4,B.Trimmed.ToAR.4,evi_DN,evi

783435|2690025||114|0.228143006540123|56|0.106632395012542|28|
0.0712654621906476|0.296610169491526|6952.80543731566

# doing the math in R
# with 0.228143006540123, 0.106632395012542 and 0.0712654621906476
# gives...

blue <- 0.228143006540123
red <- 0.106632395012542
nir <- 0.0712654621906476
evi <- ( 2.5 * ( nir - red) ) / (nir + 6 * red - 7 * blue + 1.0)
evi

[1] -0.7751909

Right, this is nice!

Using r.mapcalc

how? You mean for the specific coordinates/values (I guess... )

I get the same results like in R, not
-5905.44171917482 or 6952.80543731566, respectively. Why are the way
out of range values correct if R produces reasonable results?

Sorry, too many things to work-out these days -- I wasn't clear I guess. I
meant that "r.mapcalc" can't be wrong! I.e., I manually derive for the whole
extent of the L7 bands in question, an EVI:

r.mapcalc "evi_manual = ( 2.5 * (B.Trimmed.ToAR.4 - B.Trimmed.ToAR.3) ) /
(B.Trimmed.ToAR.4 + 6.0 * B.Trimmed.ToAR.3 - 7.5 * B.Trimmed.ToAR.1 + 1)"

This results in
r.info -r evi_manual

min=-5905.44171917482
max=6952.80543731566

What am I doing wrong?

Thanks, N

16

On Wed, June 19, 2013 23:34, Nikos Alexandris wrote:

Moritz Lennert wrote:

>> > g.region -pagc e=784935 n=2695215 rows=1 cols=1 res=30
> --------------------^^^^^^---^^^^^^^------^------^-----^^

>> > n=2695230
>> > s=2621340
>> > w=658560
>> > e=784950

> ---------vvvvv
>> > nsres=73890
>> > ewres=126390
> ---------^^^^^^

rows and cols has precedence over res, but res=30 plus -a still has an
effect, so you get an extension that is in multiples of 30, with a
resolution equal to the total extension, since you asked for one row and
one col.

So, "e=784935" and "n=2695215" are set as requested but the "-a" flag
cut's
off a bit while aligning...

It doesn't cut off, but it extends the region until you hit a multiple of 30:

n=2695215 is extended northward to n=2695230
e=784935 is extended eastward to e=784950

Moritz

17

Nikos Alexandris wrote:

>> >> > g.region -pagc e=784935 n=2695215 rows=1 cols=1 res=30
>> >
>> > --------------------^^^^^^---^^^^^^^------^------^-----^^
>> >
>> >> > n=2695230
>> >> > s=2621340
>> >> > w=658560
>> >> > e=784950
>> >
>> > ---------vvvvv
>> >
>> >> > nsres=73890
>> >> > ewres=126390
>> >
>> > ---------^^^^^^

ML:

>> rows and cols has precedence over res, but res=30 plus -a still has an
>> effect, so you get an extension that is in multiples of 30, with a
>> resolution equal to the total extension, since you asked for one row and
>> one col.

NA:

> So, "e=784935" and "n=2695215" are set as requested but the "-a" flag
> cut's
> off a bit while aligning...

ML:

It doesn't cut off,

Bien sure! I should put it like "cut's off"... I wanted to say that it simply
modifies.

but it extends the region until you hit a multiple of 30:

n=2695215 is extended northward to n=2695230
e=784935 is extended eastward to e=784950

Milles mercis, N

18

I think I figured it out:

The EVI formula in i.vi is for MODIS. The generic formula is

G * ( nir - red) / (nir + C1 * red - C2 * blue + L)

where G is a gain factor, C1, C2 are coefficients to correct for
aerosol influences in the red band using the blue band and L is the
canopy background adjustment that addresses non-linear, differential
NIR and red radiant transfer through a canopy.

The formula is scale-dependent, i.e. dependent on the range of the
input data. Scaling the input test values to a [0, 10000] range (the
range of the MODIS input bands), I get EVI values in the range [-1,
1]. Furthermore, different coefficients might be needed for different
sensors.

Some more points:

The reference in i.vi for EVI is wrong, it must be

A.Huete, K.Didan, T.Miura, E.P.Rodriguez, X.Gao, L.G. Ferreira. 2002.
Overview of the radiometric and biophysical performance of the MODIS
vegetation indices. Remote Sensing of Environment 83 195-213.

The formula we used in R is wrong:

evi <- ( 2.5 * ( nir - red) ) / (nir + 6 * red - 7 * blue + 1.0)

it must be

evi <- ( 2.5 * ( nir - red) ) / (nir + 6 * red - 7.5 * blue + 1.0)

then I get again 6952.805

g.region rast= res= (-a) does not make sense, unless you want to
resample a raster map.

Markus M

19

Markus Metz wrote:

I think I figured it out:

The EVI formula in i.vi is for MODIS.

That's precise, EVI is MODIS specific. We should clearly describe this in the
manual (I will try to alter the respective text).

The generic formula is

G * ( nir - red) / (nir + C1 * red - C2 * blue + L)

where G is a gain factor, C1, C2 are coefficients to correct for
aerosol influences in the red band using the blue band and L is the
canopy background adjustment that addresses non-linear, differential
NIR and red radiant transfer through a canopy.

Right.

The formula is scale-dependent, i.e. dependent on the range of the
input data. Scaling the input test values to a [0, 10000] range (the
range of the MODIS input bands), I get EVI values in the range [-1,
1].

Cool! Sorry that I misused the index :-/

Furthermore, different coefficients might be needed for different
sensors.

Some more points:

The reference in i.vi for EVI is wrong, it must be

A.Huete, K.Didan, T.Miura, E.P.Rodriguez, X.Gao, L.G. Ferreira. 2002.
Overview of the radiometric and biophysical performance of the MODIS
vegetation indices. Remote Sensing of Environment 83 195-213.

O-K, we need to update this as well -- noted!

The formula we used in R is wrong:

evi <- ( 2.5 * ( nir - red) ) / (nir + 6 * red - 7 * blue + 1.0)

I guess I chewed-up the .5 at some point :D, otherwise I have used it
correctly in the beginning of my tests.

it must be

evi <- ( 2.5 * ( nir - red) ) / (nir + 6 * red - 7.5 * blue + 1.0)

then I get again 6952.805

g.region rast= res= (-a) does not make sense, unless you want to
resample a raster map.

O-K.

Great, thanks!

20

On Thu, Jun 20, 2013 at 6:24 PM, Nikos Alexandris
<nik@nikosalexandris.net> wrote:

Markus Metz wrote:

I think I figured it out:

The EVI formula in i.vi is for MODIS.

That's precise, EVI is MODIS specific. We should clearly describe this in the
manual (I will try to alter the respective text).

From the literature, I gor the impression that EVI can be calculated

from other sensors as well, as long as you get the coefficients right.

The generic formula is

G * ( nir - red) / (nir + C1 * red - C2 * blue + L)

where G is a gain factor, C1, C2 are coefficients to correct for
aerosol influences in the red band using the blue band and L is the
canopy background adjustment that addresses non-linear, differential
NIR and red radiant transfer through a canopy.

Assuming that the input to i.vi should be properly preprocessed bands
with a theoretically maximum range of [0, 1], you could set L to
0.0001 and would get reasonable EVI values, sensor-independent.

Markus M