Unexpected Results from Formulas.

Errors and unexpected results
Post Reply
User avatar
joleenf
Posts: 1120
Joined: Mon Jan 19, 2009 7:16 pm

Unexpected Results from Formulas.

Post by joleenf »

Hi,

Please see the code below. I have noticed some odd results in NDVI results and square roots of data that has units of "%" or a Universal Unit.

1.) For division, when a result should be unitless, the unit of percent is retained. This leads to an unexpected magnitude.
2.) I can't explain what is happening with the Universal Unit, it does not make much sense to me. Please see code below and if you have questions about the problem, we can discuss that. There seems to be a division by 100 when using he square root function and I don't understand why.
3.) The results seem as expected if the unit is forced to None before calculation. In a whatType(data), the unit is reported as null.

Code: Select all


parmsList=dict(
    #accounting = ('YOU',9999),    #accounting information required
                                   #if ADDE dataset is not registered
                                   #in ADDE Manager
    server='dcarchivex.ssec.wisc.edu',
    dataset='GREX',
    descriptor='FD',
    band=1,
    day='2016277',
    time=('16:04', '16:09'),
    position='ALL'
    )

dateTimes=listADDEImageTimes(**parmsList)

if (len(dateTimes) >= 1):
    #panel=buildWindow(height=900,width=900)[0]

    for dateTime in dateTimes[::2]:
       index=1
       parms=dict(
           #accounting = ('YOU',9999),    #accounting information required
                                   #if ADDE dataset is not registered
                                   #in ADDE Manager
           server=parmsList['server'],
           dataset=parmsList['dataset'],
           descriptor='FD',
           day=dateTime['day'],
           time=dateTime['time'],
           unit='ALB',
           location=(0.1,-98.4),
           coordinateSystem=LATLON,
           place=ULEFT,
           mag=(1,1)
       )
       rdata=loadADDEImage(band=2, size=(8,8), **parms)
   
       gdata=loadADDEImage(band=3, size=(4,4), **parms)

       g_ds=getDomainSet(gdata)
       rdata=resample(rdata,g_ds)
       
       testSubtraction=sub(gdata,rdata)
       testAddition=sub(gdata,-rdata)
       normalizedField=testSubtraction/testAddition

       print (gdata.getRangeUnits()[0][0])
       print 'Green Band Min/Max:  {}/{} {}\n'.format(getMinMax(gdata)[0], getMinMax(gdata)[1], '%')
       
       sqrt_gdata=sqrt(gdata)
       print 'Sqrt of Green Band Min/Max requires unit to make sense {}/{} {}\n'.format(getMinMax(sqrt_gdata)[0], getMinMax(sqrt_gdata)[1],'%')
       print 'The NDVI unit does not make sense, because the result should be unitless.  If this result is unitless, then the magnitude is incorrect.'
       print 'NDVI Min/Max should be between -1/1 but:  {}/{} {}\n'.format(getMinMax(normalizedField)[0],getMinMax(normalizedField)[1] , '%')

       newRed=noUnit(rdata)
       newGreen=noUnit(gdata)
       
       testSubtraction=sub(newGreen,newRed)
       testAddition=sub(newGreen,-newRed)
       normalizedField=testSubtraction/testAddition
       
       print whatType(newGreen)
       print 'Green Band Min/Max:  {}/{} {}'.format(getMinMax(newGreen)[0], getMinMax(newGreen)[1], '1')
       sqrt_rdata=sqrt(newRed)
       print 'If the unit is removed, sqrt should be 1.11/2.21, correct?  Try using noUnit'
       print 'Sqrt of Green Band Min/Max divides by 100 (not the numbers expected)?  {}/{}'.format(getMinMax(sqrt_gdata)[0], getMinMax(sqrt_gdata)[1],'1')
       print 'NDVI Min/Max should be between -1/1 and this seems to be okay:  {}/{} {}'.format(getMinMax(normalizedField)[0],getMinMax(normalizedField)[1] , '1')

       print 'Try making the unit None'
       noneUnitRed=newUnit(rdata,'Band2_ALB',None)
       noneUnitGreen=newUnit(gdata,'Band3_ALB',None)

       testSubtraction=sub(noneUnitGreen,noneUnitRed)
       testAddition=sub(noneUnitGreen,-noneUnitRed)
       normalizedField=testSubtraction/testAddition

       print whatType(noneUnitGreen)
       print 'Green Band Min/Max:  {}/{} {}'.format(getMinMax(noneUnitGreen)[0], getMinMax(noneUnitGreen)[1], 'Null')
       sqrt_rdata=sqrt(noneUnitRed)
       print 'If the unit is removed, sqrt should be 1.11/2.21, correct?  Try using noUnit'
       print 'Sqrt of Green Band Min/Max divides by 100?  {}/{}'.format(getMinMax(sqrt_gdata)[0], getMinMax(sqrt_gdata)[1],'Null')
       print 'NDVI Min/Max should be between -1/1 and this seems to be okay:  {}/{} {}'.format(getMinMax(normalizedField)[0],getMinMax(normalizedField)[1] , 'Null')
       



Thanks,
Joleen
User avatar
bobc
Posts: 938
Joined: Mon Nov 15, 2010 5:57 pm

Re: Unexpected Results from Formulas.

Post by bobc »

Hi Joleen -

Thanks for the script and the thorough description of the problem. I'm replicating what you are seeing. For example, here are some probe values after displaying the individual layers:

red data layer: 3.4 %
green data layer: 3.8 %
green square root layer: .193
sub layer: 0.47 %
add layer: 7.2 %
normalized layer: 578.5 %

There are a few things going on here:
  1. Add and subtract look to be working as expected. The values and units of % match up and are correct.
  2. I agree that the behavior of standard deviation isn't consistent with add/subtract. You are correct that the green layer value is being divided by 100. For example, sqrt(3.8/100) is 0.194.
  3. With normalized layer, there looks to be a bug with the division. What it looks like is going on is that the sub and add data objects are divided as-is and then that result is multiplied by 10,000.

I'll keep investigating this a bit further and then write up any inquiries that are needed.

Thanks again -
Bob Carp
User avatar
joleenf
Posts: 1120
Joined: Mon Jan 19, 2009 7:16 pm

Re: Unexpected Results from Formulas.

Post by joleenf »

Hi Bob,

I am also seeing inconsistent results with rescale. I have attached the code and an image of the result in my jython shell. This unfortunately does not consistently break.

Code: Select all


parms=dict(
        server='geoarc.ssec.wisc.edu',
        dataset='QHIM08',
        #accounting=()    #accounting information is required if this ADDE dataset is not already registered in your ADDE Manager
        unit='ALB',
        day='2016-10-6',
        time='4:30',
        size=(1980,3104),
        coordinateSystem=AREA,
        place=CENTER,
        location=(6820,3740),
        mag=(-9,-9)
       
)
   
rdata=loadADDEImage(band=3, descriptor='FD03', **parms)  #8891 RED

print whatType(rdata)
rng_rdata=(GridUtil.fieldMinMax(rdata))[0]
print '% - Range Min: {}, Range Max: {}, Result of getMinMax:  {}/{}'.format(rng_rdata.min, rng_rdata.max, getMinMax(rdata)[0],getMinMax(rdata)[1])
new_rdata=rescale(rdata,rng_rdata.min,rng_rdata.max,0,100)
print 'AFTER RESCALE: {}/{}'.format(getMinMax(new_rdata)[0],getMinMax(new_rdata)[1])


red=newUnit(rdata,'Band3_ALB',"1")
rng_data=(GridUtil.fieldMinMax(red))[0]
print whatType(red)
print 'Universal Unit - Range Min: {}, Range Max: {}, Result of getMinMax:  {}/{}'.format(rng_data.min, rng_data.max, getMinMax(red)[0], getMinMax(red)[1])
new_red=rescale(red,rng_data.min,rng_data.max,0,100)
print 'AFTER RESCALE: {}/{}'.format(getMinMax(new_red)[0], getMinMax(new_red)[1])

red1=newUnit(rdata,'Band3_ALB',None)
print whatType(red1)
rng_data1=(GridUtil.fieldMinMax(red1))[0]
print 'Null Unit - Range Min: {}, Range Max: {}, Result of getMinMax:  {}/{}'.format(rng_data1.min, rng_data1.max, getMinMax(red1)[0], getMinMax(red1)[1])
new_red1=rescale(red1,rng_data1.min,rng_data1.max,0,100)
print 'AFTER RESCALE: {}/{}'.format(getMinMax(new_red1)[0],getMinMax(new_red1)[0])

Screen Shot 2017-01-31 at 10.29.06 AM.png


And from this code run just a second before the previous code:

Code: Select all

parms=dict(
        server='geoarc.ssec.wisc.edu',
        dataset='QHIM08',
        #accounting=()    #accounting information is required if this ADDE dataset is not already registered in your ADDE Manager
        unit='ALB',
        day='2016-10-6',
        time='4:30',
        size=(1980,3104),
        coordinateSystem=AREA,
        place=CENTER,
        location=(6820,3740),
        mag=(1,1)
       
)
   
rdata=loadADDEImage(band=3, descriptor='FD03', **parms)  #8891 RED

print whatType(rdata)
rng_rdata=(GridUtil.fieldMinMax(rdata))[0]
print '% - Range Min: {}, Range Max: {}, Result of getMinMax:  {}/{}'.format(rng_rdata.min, rng_rdata.max, getMinMax(rdata)[0],getMinMax(rdata)[1])
new_rdata=rescale(rdata,rng_rdata.min,rng_rdata.max,0,100)
print getMinMax(new_rdata)


red=newUnit(rdata,'Band3_ALB',"1")
rng_data=(GridUtil.fieldMinMax(red))[0]
print whatType(red)
print 'Universal Unit - Range Min: {}, Range Max: {}, Result of getMinMax:  {}/{}'.format(rng_data.min, rng_data.max, getMinMax(red)[0], getMinMax(red)[1])
new_red=rescale(red,rng_data.min,rng_data.max,0,100)
print getMinMax(new_red)

red1=newUnit(rdata,'Band3_ALB',None)
print whatType(red1)
rng_data1=(GridUtil.fieldMinMax(red1))[0]
print 'Null Unit - Range Min: {}, Range Max: {}, Result of getMinMax:  {}/{}'.format(rng_data1.min, rng_data1.max, getMinMax(red1)[0], getMinMax(red1)[1])
new_red1=rescale(red1,rng_data1.min,rng_data1.max,0,100)
print getMinMax(new_red1)

Screen Shot 2017-01-31 at 10.29.21 AM.png
User avatar
bobc
Posts: 938
Joined: Mon Nov 15, 2010 5:57 pm

Re: Unexpected Results from Formulas.

Post by bobc »

Hi Joleen -

Thanks again for the scripts. I'm replicating that having a unit of '1' seems to give unexpected results, both before and after the rescale. Before rescaling, getMinMax() isn't returning the correct values (they are being divided by 100). After the rescale, getMinMax() again isn't returning the correct values. I wrote this up, along with the problem involving division from your earlier post, as Inquiry 2490.

Thanks -
Bob
Post Reply