McIDAS-V Support Forums

Support and discussion of the McIDAS-V software package
Forum HomeMcIDAS-V Home
It is currently 04 Aug 2020, 02:54

All times are UTC




Post new topic Reply to topic  [ 4 posts ] 
Author Message
 Post subject: Unexpected Results from Formulas.
PostPosted: 30 Jan 2017, 16:31 
User avatar

Joined: 19 Jan 2009, 19:16
Posts: 1119
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:

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


Top
 Profile  
 
 Post subject: Re: Unexpected Results from Formulas.
PostPosted: 30 Jan 2017, 20:57 
User avatar

Joined: 15 Nov 2010, 17:57
Posts: 917
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


Top
 Profile  
 
 Post subject: Re: Unexpected Results from Formulas.
PostPosted: 31 Jan 2017, 16:31 
User avatar

Joined: 19 Jan 2009, 19:16
Posts: 1119
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:

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])

Attachment:
Screen Shot 2017-01-31 at 10.29.06 AM.png
Screen Shot 2017-01-31 at 10.29.06 AM.png [ 107.93 KiB | Viewed 313 times ]


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

Code:
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)

Attachment:
Screen Shot 2017-01-31 at 10.29.21 AM.png
Screen Shot 2017-01-31 at 10.29.21 AM.png [ 109.25 KiB | Viewed 318 times ]


Top
 Profile  
 
 Post subject: Re: Unexpected Results from Formulas.
PostPosted: 02 Feb 2017, 21:34 
User avatar

Joined: 15 Nov 2010, 17:57
Posts: 917
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


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 4 posts ] 

All times are UTC


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group