Calculating EVI from Landsat 8 - strange results

190 views
Skip to first unread message

Benjamin Mewes

unread,
Feb 2, 2017, 2:34:36 AM2/2/17
to RSGISLib Support
Dear all,

as I want to calculate the EVI (Enhances Vegetation Index), I tried to accomplish it via RSGISLIB and a Landsat 8 scene after transfering the scene to radiance (rad) by using Arcsi.py.

The calculation of the NDVI works like a charm for these _rad.kea files, but the EVI after Huete et al. results in very strange rasters which have values > 1000! If I force QGIS only to show values between 0. and 1. the results seem to be fine. Any idea what might go wrong?

Here is my code:

bandDefns = []
bandDefns
.append(imagecalc.BandDefn('nir', inputStack, imageutils.getBandNames(inputStack).index("NIR")))
bandDefns
.append(imagecalc.BandDefn('red', inputStack, imageutils.getBandNames(inputStack).index("Red")))
bandDefns
.append(imagecalc.BandDefn('blue', inputStack, imageutils.getBandNames(inputStack).index("Blue")))

outputImage
= "L8_evi.kea"
datatype_evi
= rsgislib.TYPE_32FLOAT
expression_evi
= "2.5*(nir-red)/(nir+(6.0 * red) - (7.5 * blue)+1.0)"
imagecalc
.bandMath(outputImage, expression_evi, gdalformat, datatype_evi, bandDefns)
imageutils
.popImageStats(outputImage, usenodataval=True,nodataval=999, calcpyramids=True)



Pete Bunting [pfb]

unread,
Feb 2, 2017, 4:20:22 AM2/2/17
to Benjamin Mewes, RSGISLib Support
Hi Benjamin, 

I assume these values are in no data regions of the input image? If so, then add an if statement to your band maths expression. 

The bandmath expression using muparser (http://beltoforion.de/article.php?a=muparser) to parse the expression. All syntax is here:


For this expression, check if all inputs are 0 (therefore no data) and therefore output 0 into the output image otherwise do calculation.

expression_evi = “(nir==0)&&(red==0)&&(blue==0)?0:2.5*(nir-red)/(nir+(6.0 * red) - (7.5 * blue)+1.0)"

Note the inline c syntax 

Operator Description Remarks
?: if then else operator C++ style syntax


Best wishes, 

Pete

****************************************************
* Dr Pete Bunting
* Reader in Remote Sensing
* Earth Observation and Ecosystem Dynamics Group
* Department of Geography and Earth Sciences
* Aberystwyth University
* Aberystwyth
* Ceredigion
* SY23 3DB
* UK

* Ph: +44 (0) 1970 622615
* Mob: +44 (0) 7917 842743
* Email: p...@aber.ac.uk
* ORCID: http://orcid.org/0000-0002-7435-0148
****************************************************
"Please consider the environment before printing this email or any documents attached”

--
You received this message because you are subscribed to the Google Groups "RSGISLib Support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to rsgislib-suppo...@googlegroups.com.
To post to this group, send email to rsgislib...@googlegroups.com.
Visit this group at https://groups.google.com/group/rsgislib-support.
For more options, visit https://groups.google.com/d/optout.



--------------------------------------------------------------------
Un o’r 4 prifysgol uchaf yn y DU a’r orau yng Nghymru am fodlonrwydd myfyrwyr.
(Arolwg Cenedlaethol y Myfyrwyr 2016)
www.aber.ac.uk

Top 4 UK university and best in Wales for student satisfaction
(National Student Survey 2016)
www.aber.ac.uk

Benjamin Mewes

unread,
Feb 2, 2017, 7:30:01 AM2/2/17
to RSGISLib Support, bennam...@gmail.com, p...@aber.ac.uk
Hi Pete,

thanks for your quick answer. You were right about some assumptions, but it didn't help eventually. Now, the EVI has NAN values inside (leading QGis to a crash). What makes me curious: NDVI and Green Index are properly calculated now, but EVI and consequently LAI lead to NAN in my grids. I tried the subparts of the expression showing that (nir-red) leads to nan in the results, which even more strange because NDVI calculates correctly. Is there some kind of "bad" data leading to a nan?

Pete Bunting [pfb]

unread,
Feb 2, 2017, 7:41:06 AM2/2/17
to RSGISLib Support, Benjamin Mewes, bennam...@gmail.com
Hi Benjamin,

That does sound strange, as you say if it works for ndvi it seems very odd it doesn't work for others.

Few things to check:

1. Double check input image is correct with appropriate values.
2. I see you are using  the band names list to define the band numbers. Is that doing what you think? The list will be from 0 while bands are numbered from 1. Maybe try hard coding the band to ensure there is not an error there.

If that doesn't reveal anything then can you send me the image and script to test?

Many thanks,

Pete



From: Benjamin Mewes <bennam...@gmail.com>
Sent: Thursday, February 2, 2017 12:30:00 PM
To: RSGISLib Support
Cc: bennam...@gmail.com; Pete Bunting [pfb]
Subject: Re: [rsgislib-support] Calculating EVI from Landsat 8 - strange results
 
Hi Pete,

thanks for your quick answer. You were right about some assumptions, but it didn't help eventually. Now, the EVI has NAN values inside (leading QGis to a crash). What makes me curious: NDVI and Green Index are properly calculated now, but EVI and consequently LAI lead to NAN in my grids. I tried the subparts of the expression showing that (nir-red) leads to nan in the results, which even more strange because NDVI calculates correctly. Is there some kind of "bad" data leading to a nan?


Benjamin Mewes

unread,
Feb 2, 2017, 8:03:34 AM2/2/17
to RSGISLib Support, bennam...@gmail.com, p...@aber.ac.uk
Hi Pete,

@1) the input image is correct.
@2) I tried to hard code it:

bandDefns = []
bandDefns
.append(imagecalc.BandDefn('red', inputStack, 4))
bandDefns
.append(imagecalc.BandDefn('nir', inputStack, 5))
bandDefns
.append(imagecalc.BandDefn('blue', inputStack, 2))


Once again it failed...

I uploaded the data for you here:
https://www.dropbox.com/sh/w1z8248nepmie3a/AACX_htd_G3BgOSsuqs_EDmZa?dl=0

Benjamin Mewes

unread,
Feb 3, 2017, 3:38:54 AM2/3/17
to RSGISLib Support, bennam...@gmail.com, p...@aber.ac.uk
Hi Pete,

i found out, that just some pixel values within my mask are Inf or -Inf the others are computed correctly. I will try to eliminate these Inf values.

Pete Bunting [pfb]

unread,
Feb 3, 2017, 3:43:26 AM2/3/17
to RSGISLib Support, Benjamin Mewes, bennam...@gmail.com
Hi Benjamin,

There is a function to create a mask of finite image pixel values.

rsgislib.imageutils.genFiniteMask(inimage=string, outimage=string, format=string)

Best wishes

Pete

From: Benjamin Mewes <bennam...@gmail.com>
Sent: Friday, February 3, 2017 8:38:54 AM

Benjamin Mewes

unread,
Feb 3, 2017, 4:32:25 AM2/3/17
to RSGISLib Support, bennam...@gmail.com, p...@aber.ac.uk
Hi,

thanks, that did it!


Pete Bunting [pfb]

unread,
Feb 3, 2017, 4:38:42 AM2/3/17
to Benjamin Mewes, RSGISLib Support
Hi Benjamin, 

The attached script works OK for me. 

best wishes, 

Pete



****************************************************
* Dr Pete Bunting
* Reader in Remote Sensing
* Earth Observation and Ecosystem Dynamics Group
* Department of Geography and Earth Sciences
* Aberystwyth University
* Aberystwyth
* Ceredigion
* SY23 3DB
* UK

* Ph: +44 (0) 1970 622615
* Mob: +44 (0) 7917 842743
* Email: p...@aber.ac.uk
* ORCID: http://orcid.org/0000-0002-7435-0148
****************************************************
"Please consider the environment before printing this email or any documents attached”

02_calc_indices.py

Benjamin Mewes

unread,
Feb 3, 2017, 5:00:45 AM2/3/17
to RSGISLib Support, bennam...@gmail.com, p...@aber.ac.uk
Hi,

same error here: QGIS finds Inf values which leads to a crash. So I decided to use the finite Mask which works now!

Pete Bunting [pfb]

unread,
Feb 3, 2017, 5:04:41 AM2/3/17
to Benjamin Mewes, RSGISLib Support
Hi Benjamin, 

Great that it works, strange inf values cause QGIS to crash - sounds like a bug in QGIS. 

Best wishes, 

Pete

****************************************************
* Dr Pete Bunting
* Reader in Remote Sensing
* Earth Observation and Ecosystem Dynamics Group
* Department of Geography and Earth Sciences
* Aberystwyth University
* Aberystwyth
* Ceredigion
* SY23 3DB
* UK

* Ph: +44 (0) 1970 622615
* Mob: +44 (0) 7917 842743
* Email: p...@aber.ac.uk
* ORCID: http://orcid.org/0000-0002-7435-0148
****************************************************
"Please consider the environment before printing this email or any documents attached”

--
You received this message because you are subscribed to the Google Groups "RSGISLib Support" group.
To unsubscribe from this group and stop receiving emails from it, send an email to rsgislib-suppo...@googlegroups.com.
To post to this group, send email to rsgislib...@googlegroups.com.
Visit this group at https://groups.google.com/group/rsgislib-support.
For more options, visit https://groups.google.com/d/optout.
Reply all
Reply to author
Forward
0 new messages