how to indentify infinite values in array

[joleenf]

Hi,

I recently have had some problems with and NDVI calculation in which I am seeing min/max bounds from -inf/inf. I decided to try to diagnose the problem by using the "find" command to identify points in the array which are a problem. However, I found that "find" cloud not identify infinite values properly. When the product is displayed, values are listing as missing to missing (actually, I can reset the range to -1,1 and see an NDVI result, so the array is not full of missing or even infinite values). I tried the following tests to identify infinite values in an array with no success.

Code: Select all

import java.lang.Float

infinity=Float.POSITIVE_INFINITY

afield=field([infinity,infinity,1])

print afield

l= find(afield,'==',infinity)

print len(l)

# this fails to identify locations of infinite value as well (find nothing)

l=find(afield,'==', float('inf'))
print len(l)

# can't get the expected answer when looking for the opposite either (would expect an array with dimension of 1)

l=find(afield,'!=', infinity)

print l

Screen Shot 2017-02-06 at 8.44.45 AM.png (9.64 KiB) Viewed 5080 times

Joleen

[bobc]

Hi Joleen -

I talked this over with a programmer, and it is likely a limitation of find that it can't make use of 'infinity' since infinity isn't a real number. NaN has the same issue.

I found that there are a couple of ways you can return the max value of your 'afield' example:

Code: Select all

print getMinMax(afield)
# returns:
# array('d', [1.0, inf])

describe(afield)
# returns a variety of parameters including:
# Max: Infinity

These two methods allow you to determine if 'infinity' is present in your data object.

Does this give you what you are looking for?

Thanks -
Bob Carp

[joleenf]

Hi Bob,

I know infinite values are in my data and they should not be there. I need to know why there are infinite values and I can't identify the problem without knowing where it is happening. I have tried removing denominators that are incredibly small. This did not remove the -inf/inf problem either.

Joleen

[bobc]

Hi Joleen -

Can you please send me the NDVI code/formula you are using? Also, what data are you using with this formula, is it AHI/ABI? With the code/data I'll see if I can replicate the problem on my end.

Thanks -
Bob

[joleenf]

Hi Bob,

I found there was a problem when I used the fixed grid files of AHI produced by Ray (please see email with Re: Newest AHI from Ray) and the latest NavUtility.normalizedDifference formula. I then tried to use the "traditional" method of calculating AHI with

redData=loadGrid()
grnData=loadGrid()

ds=getDomainSet(grnData)
newRed=resample(redData,ds) # I don't use resampleGrid because I can never remember the order

ndvi=sub(grnData,newRed)/sub(grnData,-newRed)

A work-around was found by a colleague, so I can at least set the infinite values to missing using

ndvi=ndvi/abs(ndvi)

Joleen

[bobc]

Hi Joleen -

The first thing I did was create a Grid Table display of redData and grnData and I saw that there are a lot of non-navigated points that still have albedo values. To set these non-navigated points to missing, I used the setMissingNoNavigation JPythonMethod on both redData and grnData (note that all of the print statements below aren't necessary, I just used them to let me know when the formulas were complete):

Code: Select all

redMiss = setMissingNoNavigation(redData)
print 'redMiss completed'
grnMiss = setMissingNoNavigation(grnData)
print 'grnMiss completed'

I used these redMiss and grnMiss objects later on in your code:

Code: Select all

ds = getDomainSet(grnMiss)
print 'getDomainSet completed'
newRed = resample(redMiss, ds)
print 'resample completed'

From your NDVI formula:

Code: Select all

ndvi=sub(grnData,newRed)/sub(grnData,-newRed)

I pulled out the numerator and denominator to make them their own objects (using grnMiss instead of grnData):

Code: Select all

denom = sub(grnMiss,-newRed)
numer = sub(grnMiss,newRed)
print 'done denom and numer'

After creating the numerator and denominator objects, I ran the division and replicated the -inf to +inf data value range:

Code: Select all

ndvi = numer/denom
getMinMax(ndvi)
# result: array('d', [-inf, inf])

Next, I created a Grid Table display of 'denom':

Code: Select all

layer = activeDisplay().createLayer('Grid Table',denom)

After sorting the right column in the Grid Table (which does take some time with how many values there are) I can see that there are 1,344 points with a value of 0 (zero) which is causing the infinity values:

- When you have a positive value in the numerator divided by a value of 0 the result is positive infinity.
- When you have a negative value in the numerator divided by a value of 0 the result is negative infinity.

I believe this explains the -inf to +inf range for NDVI. Do you agree?

Thanks -
Bob

[bobc]

Hi again Joleen -

Quick followup. If you want to, you can set the denominator values equal to 0 to be missing using the setToMissing JPythonMethod:

Code: Select all

denomMiss = setToMissing(denom,0)
ndviMiss = numer/denomMiss
getMinMax(ndviMiss)
# returns: array('d', [-7.000000953674316, 4.999999523162842])

I created a Grid Table display of this ndviMiss object and by far the majority of the data values are within the expected range of -1 to 1. In my case, there are 1,328 data points not in this range while there are a total of 7,562,500 data points in the object (per the describe() function). Perhaps there are some bad/unexpected albedo values in the file that are causing these values outside of the expected range?

It looks like the original redMiss data has albedo values ranging from -0.002 to 1.192 and grnMiss has albedo values ranging from -0.001 to 1.192.

- Bob

[joleenf]

Hi Bob,

Oh oops! I should have tried that.

Thank-you for the solution. I think I noticed something else about the range issue outlined in
http://mcidasv.ssec.wisc.edu/forums/vie ... f4abfd774a

This range of this ndvi should be -1 to 1. However, you see that the range is off by an order of magnitude. However, if the image is displayed and the data values are probed, it looks like the data values in the display are correct.

Joleen

[bobc]

Hi Joleen -

I was going to look into this, but I wasn't sure exactly how you are computing the NDVI. From the forum post you linked, there were a few different ways that the NDVI was being calculated:

• Use the red and green flat fields as-is without modifying the unit
• Pass the red and green flat fields through noUnit()
• Pass the red and green flat fields through newUnit()
• There are references to rescale(). Are you rescaling the red and green flat fields before doing the subtraction and division?

Which method are you using?

Thanks -
Bob

[joleenf]

Hi Bob,

I am using newUnit to remove the unit...

Code: Select all

rfield=getRangeType(rdata)
    bfield=getRangeType(bdata)
    nIRfield=getRangeType(nIR)
    
    #rUnit=(rfield.defaultUnits)[0]
    #bUnit=(bfield.defaultUnits)[0]
    #nIRUnit=(nIRfield.defaultUnits)[0]
    rUnit=whatType(rdata).toString().rstrip().splitlines()[-1].split('Unit: ',1)[1]
    nIRUnit=whatType(nIR).toString().rstrip().splitlines()[-1].split('Unit: ',1)[1]
    bUnit=whatType(bdata).toString().rstrip().splitlines()[-1].split('Unit: ',1)[1]

    rdata=newUnit(rdata,rfield.toString(),None)
    bdata=newUnit(bdata,bfield.toString(),None)
    nIR=newUnit(nIR,nIR.toString(),None)

I am rescaling after the NDVI calculation. I thought that rescale would just set values that are not within the range to outLo and outHi.

Joleen