4

I have raster function for computing the EVI index from a Sentinel-2 raster located in a mosaic dataset. The code is similar to https://github.com/Esri/raster-functions/blob/master/functions/NDVI.py I added the block for the blue band:

            {
                'name': 'blue',
                'dataType': 'numeric',
                'value': 1,
                'required': True,
                'displayName': "Blue Band Index",
                'description': "Index for Blue band"
            },

The EVI formula is:

outBlock = 2.5 * (nir - red) / ((nir + 6*red - 7.5*blue) + 1)

Obviously this does not raise a DivisionError because it's impossible due the the '+1' in the formula.

When i run the function on the ArcGIS dataset, i get the error:

 outBlock = 2.5 * (nir - red) / ((nir + 6*red - 7.5*blue) + 1)

FloatingPointError: divide by zero encountered in divide

How can I solve this?

migrated from stackoverflow.com Sep 7 '16 at 11:46

This question came from our site for professional and enthusiast programmers.

  • as mentioned by @T. Desmarez, your denominator could mathematically be equal to zero. If you work with reflectances, this is however quite unlikely, so make sure that your input data is a reflectance. – radouxju Sep 7 '16 at 12:10
  • I managed a workaround with a 'with' statement. I will post an answer – Litwos Sep 7 '16 at 12:15
4

Perhaps "(nir + 6*red - 7.5*blue)" can be equal to "-1".

  • Yes. That can be possible. I haven't thought about that.. – Litwos Aug 19 '16 at 9:18
2

I managed to do it like this:

with np.errstate(divide = 'ignore', invalid = 'ignore'): 
    outBlock[outBlock == np.inf] = 1.0
    outBlock[outBlock == -float('Inf')] = -1.0 

Whenever an infinity value is obatined, is changed to +1 or -1, so it can fit in the EVI range.

1

As pointed out by T. Desmarez, your can indeed be zero. What you can do about this is, is to calculate it separately and handle it according to your needs like this:

divisor = (nir + 6 * red - 7.5 * blue) + 1
if (divisor == o):
    outBlock = -9999
else:
    outBlock = 2.5 * (nir - red) / divisor

In this example, instead of trying to calculate the value, you assign a fixed number to outBlock.

  • I tried that, but I cannot test with 'if/else' because it needs to check on multiple arrays. 'np.all' and np.any' won't work either. I just posted my resolution. – Litwos Sep 7 '16 at 12:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.