How do I extract units consistently?

Post any questions, ideas, or topics related to Jython and Python scripting.
Post Reply
User avatar
joleenf
Posts: 1123
Joined: Mon Jan 19, 2009 7:16 pm

How do I extract units consistently?

Post by joleenf »

Hi,

I can't find a way to extract the unit consistently from the data object in McIDAS-V.

I have used:
1.) flatField.getRangeUnits()[0][0]
I have used:
1.) flatField.getRangeUnits()[0][0]
- <class '__main__._MappedAreaImageFlatField'> original unit, correct result
- noUnit(<class '__main__._MappedAreaImageFlatField'>) universal unit, incorrect result
- <class '__main__._MappedGeoGridFlatField'> correct result
- <type 'ucar.visad.data.GeoGridFlatField'> correct result (for ADDE or Grids loaded through choosers)
- <type 'visad.FieldImpl'> doesn't work, but did not expect this to work since the Impl can be an array of data fields.

2.) getRangeType(gdata)
- <class '__main__._MappedAreaImageFlatField'> incorrect result
- noUnit(<class '__main__._MappedAreaImageFlatField'>) correct result
- <class '__main__._MappedGeoGridFlatField'> different path needed
- <type 'ucar.visad.data.GeoGridFlatField'> unit is in result, but result is not a unit
(for ADDE or Grids loaded through choosers)
- <type 'visad.FieldImpl'> unit is in result, but result is not a unit

3.) GridUtil.getParamUnits(flatField)[0] does not always report the correct unit
- <class '__main__._MappedAreaImageFlatField'> correct result
- noUnit(<class '__main__._MappedAreaImageFlatField'>) incorrect result
- <class '__main__._MappedGeoGridFlatField'> correct result
- <type 'ucar.visad.data.GeoGridFlatField'> correct result (for ADDE or Grids loaded through choosers)
- <type 'visad.FieldImpl'> correct

4.) whatType always reports corrects units as far as I can tell, but how can the units be extracted?

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)
       )
   
       gdata=loadADDEImage(band=3, size=(4,4), **parms)     

       print 'gdata is {}'.format(type(gdata))
       print '\n(CORRECT) gdata.getRangeUnits()[0][0]:  {}'.format((gdata.getRangeUnits()[0][0]))
       print '(WRONG) getRangeType(gdata) will not work in this case {}'.format(getRangeType(gdata))
       print '(CORRECT) GridUtil.getParamUnits: {}\n'.format(GridUtil.getParamUnits(gdata)[0])

       print 'Apply noUnit\n\n'
       newGreen=noUnit(gdata)
       print 'newGreen is {}'.format(type(newGreen))
       print '(CORRECT) getRangeType(newGreen):  {}\n'.format(getRangeType(newGreen))
       print '(WRONG) As reported above, this fails on the original object, try getRangeType(gdata): {}'.format(getRangeType(gdata))
       
       thisNewUnit=newGreen.getRangeUnits()[0][0]
       print '(WRONG) newGreen.getRangeUnits()[0][0]:  noUnit has been applied and the unit is still being reported as {}'.format(thisNewUnit)
       print '(CORRECT) Try using whatType(newGreen)....This should report the Units correctly as UniversalUnit'
       print whatType(newGreen)
       print '(WRONG) GridUtil.getParamUnits(newGreen)[0]:  {}\n'.format(GridUtil.getParamUnits(newGreen)[0])

       print 'The result for flatFields loaded via loadGrid is different'
       filename='/Users/joleenf/data/goesr/DOE-4/CONUS/OR_ABI-L2-CMIPC-M3C07_G16_s20162231442163_e20162231444535_c20162231444592.nc'
       thisData=loadGrid(filename=filename,field='CMI',stride=50)
       
       print 'thisData is {}'.format(type(thisData))
       print '(CORRECT) thisData.getRangeUnits()[0][0]:  {}'.format(thisData.getRangeUnits()[0][0])
       print '(DIFFERENT PATH) getRangeType(thisData):  {}\n'.format(getRangeType(thisData))
       print 'Getting the unit is via a different method getRangeType(thisData).range.defaultUnit.toString():  {}\n'.format(getRangeType(thisData).defaultUnit.toString())
       print '(CORRECT):  GridUtil.getParamUnits(thisData): {}\n'.format(GridUtil.getParamUnits(thisData)[0])

       a=selectData('Stage same grid as before via grid chooser')
       print 'a is {}'.format(type(a))
       print 'For netCDF grid data loaded via the gui selection and then accessed via select data, getRangeType produces a different answer\n'
       
       print '(NOT QUITE A UNIT) Try getRangeType(a):  {}'.format(getRangeType(a))
       print '(CORRECT BUT DIFFERENT PATH) longer path provides the unit getRangeType(a).range.defaultUnit: {}'.format(getRangeType(a).range.defaultUnit)
       print '(CORRECT) GridUtil.getParamUnits(a): {}\n'.format(GridUtil.getParamUnits(a)[0])

       
       print 'a[0] is {}'.format(type(a[0]))
       print '(CORRECT) a[0].getRangeUnits()[0][0]:  {}\n'.format(a[0].getRangeUnits()[0][0])
       print '(NOT QUITE A UNIT) Try getRangeType(a[0]):  {}.  Is there a path to get the unit from this?'.format(getRangeType(a[0]))
       print '(CORRECT) GridUtil.getParamUnits(a[0]): {}\n'.format(GridUtil.getParamUnits(a[0])[0])


       print 'Try an ADDE image'
       b=selectData('Use ADDE Chooser')
       print 'b[0] is {}'.format(type(a[0]))
       print '(CORRECT) b[0].getRangeUnits()[0][0]:  {}\n'.format(b[0].getRangeUnits()[0][0])
       print '(NOT QUITE A UNIT) Try getRangeType(b[0]):  {}.  Is there a path to get the unit from this?'.format(getRangeType(b[0]))
       print '(CORRECT) GridUtil.getParamUnits(b[0]): {}\n'.format(GridUtil.getParamUnits(b[0])[0])

      print 'whatType is consistent with expected units, but need to extract units'
      print whatType(gdata)
      print whatType(newGreen)
      print whatType(thisData)
      print whatType(a)
      print whatType(b)
      print whatType(a[0])
      print whatType(b[0])


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

Re: How do I extract units consistently?

Post by bobc »

Hi Joleen -

As for extracting the units from whatType(), I found a solution (although I'm not sure if it is the most elegant way of doing it):

Code: Select all

print whatType(gdata).toString().rstrip().splitlines()[-1].split('Unit: ',1)[1]
print whatType(newGreen).toString().rstrip().splitlines()[-1].split('Unit: ',1)[1]
print whatType(thisData).toString().rstrip().splitlines()[-1].split('Unit: ',1)[1]
print whatType(a).toString().rstrip().splitlines()[-1].split('Unit: ',1)[1]
print whatType(b).toString().rstrip().splitlines()[-1].split('Unit: ',1)[1]
print whatType(a[0]).toString().rstrip().splitlines()[-1].split('Unit: ',1)[1]
print whatType(b[0]).toString().rstrip().splitlines()[-1].split('Unit: ',1)[1]

Here's a description of the above lines:
  • toString - Change what is returned from whatType to a string.
  • rstrip - Remove any blank lines or characters at the end of the string.
  • splitlines - Pull the last line from the string. In some cases, whatType returns multiple lines that includes the string 'Unit:', and the last line contains the unit we are interested in.
  • split - Pull the text that comes after the string 'Unit: '.

Does this get you what you are looking for?

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

Re: How do I extract units consistently?

Post by joleenf »

Hi Bob,

Yes and no. I don't know how to look into where whatType() is actually getting the unit. There is probably a more direct path that would provide the same consistency. I also wonder if this inconsistency is to blame for some of the odd behavior noted in the formulas, namely whatever is deciding to magnitudes to use in the formulas is extracting the wrong unit and in some cases maintaining a unit when the result should become unit-less:

http://mcidasv.ssec.wisc.edu/forums/vie ... d753e996e1

When the formula inconsistencies are taken into account, this becomes a more critical problem.

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

Re: How do I extract units consistently?

Post by bobc »

Hi Joleen -

I just wanted to let you know that I sent this around to a few programmers. I tried looking up documentation of how units are returned from whatType(), but I haven't been able to find anything. I agree that this looks to be the one consistent way to return the units of a data object. When I hear back from the programmers, I'll follow up with you.

Thanks -
Bob
User avatar
bobc
Posts: 988
Joined: Mon Nov 15, 2010 5:57 pm

Re: How do I extract units consistently?

Post by bobc »

Hi Joleen -

I heard back from a programmer about how whatType returns the unit, and here is what he said:
I think the trouble is stemming from the amount of flexibility VisAD can accommodate; dumpMT looks like it basically recurses through things until it hits a RealType that has a non-null return value for RealType.getDefaultUnit().

FYI, dumpMT is referring to the the code found in visad.jar (/visad/jmet/DumpType.java) that controls the whatType output:

Code: Select all

  private static void dumpMT(MathType t, String prefix) {
    PrintStream out = os instanceof PrintStream ?
      (PrintStream) os : new PrintStream(os);
    if (init) {
      out.println("VisAD MathType analysis");
    }
    init = false;

    try {

      if (t instanceof FunctionType) {
        out.println(prefix + " FunctionType: ");
        RealTupleType domain = ((FunctionType) t).getDomain();
        int           num    = domain.getDimension();
        out.println(prefix + " Domain has " + num + " components:");
        for (int i = 0; i < num; i++) {
          MathType comp = domain.getComponent(i);
          dumpMT(comp, prefix + "  " + i + ".");
        }

        out.println(prefix + " Range:");
        MathType range = ((FunctionType) t).getRange();
        dumpMT(range, prefix + "  ");


      } else if (t instanceof SetType) {
        out.println(prefix + " SetType: " + t);

      } else if (t instanceof RealTupleType) {
        int num = ((RealTupleType) t).getDimension();
        out.println(prefix + " RealTupleType has " + num
                           + " components:");
        for (int i = 0; i < num; i++) {
          MathType comp = ((RealTupleType) t).getComponent(i);
          dumpMT(comp, prefix + "  " + i + ".");
        }


      } else if (t instanceof TupleType) {
        int num = ((TupleType) t).getDimension();
        out.println(prefix + " TupleType has " + num + " components:");
        for (int i = 0; i < num; i++) {
          MathType comp = ((TupleType) t).getComponent(i);
          dumpMT(comp, prefix + "  " + i + ".");
        }

      } else if (t instanceof TextType) {
        out.println(prefix + " TextType: " + t);

      } else if (t instanceof RealType) {
        out.println(prefix + " RealType: " + t);
        prefix = prefix + "  ";
        out.println(prefix + " Name = " + ((RealType) t).toString());
        Unit   du = ((RealType) t).getDefaultUnit();
        String s  = null;
        if (du != null) {
          s = du.toString();
        }
        if (s != null) {
          out.println(prefix + " Unit: " + s);

Our programmer attempted replicating this with the following Jython code that you can put in your Jython Library:

Code: Select all

from visad import Data
from visad import FunctionType
from visad import RealType
from visad import TupleType

def _onlyrange(mathtype, units, inrange):
    #print("Type: {0} str: {1} units: {2} inrange: {3}".format(type(mathtype), mathtype, units, inrange))
    if isinstance(mathtype, FunctionType):
        domain = mathtype.getDomain()
        dims = domain.getDimension()
        for dim in range(dims):
            comp = domain.getComponent(dim)
            _onlyrange(comp, units, False)
        rng = mathtype.getRange()
        _onlyrange(rng, units, True)
    elif isinstance(mathtype, TupleType):
        dims = mathtype.getDimension()
        for dim in range(dims):
            comp = mathtype.getComponent(dim)
            _onlyrange(comp, units, inrange)
    elif isinstance(mathtype, RealType):
        unit = mathtype.getDefaultUnit()
        #print("SHOULD WORK: Type: {0} str: {1} units: {2} inrange: {3}".format(type(mathtype), mathtype, units, inrange))
        if unit is not None and inrange:
            units.append(unit)
           
def _getunit(mathtype, units):
    if isinstance(mathtype, FunctionType):
        domain = mathtype.getDomain()
        dims = domain.getDimension()
        for dim in range(dims):
            comp = domain.getComponent(dim)
            _getunit(comp, units)
        rng = mathtype.getRange()
        _getunit(rng, units)
    elif isinstance(mathtype, TupleType):
        dims = mathtype.getDimension()
        for dim in range(dims):
            comp = mathtype.getComponent(dim)
            _getunit(comp, units)
    elif isinstance(mathtype, RealType):
        unit = mathtype.getDefaultUnit()
        if unit is not None:
            units.append(unit)
           
def getUnitsOrElse(d, rangeOnly=False):
    if not isinstance(d, Data):
        raise TypeError("Unknown type: {0}".format(type(d)))
    units = []
    if rangeOnly:
        _onlyrange(d.getType(), units, False)
    else:
        _getunit(d.getType(), units)
    return units

You can call this getUnitsOrElse function within your script as such:

Code: Select all

print getUnitsOrElse(gdata)
# output = [%]
print getUnitsOrElse(newGreen)
# output = [UniversalUnit]
print getUnitsOrElse(a[0])
# output = [rad, rad, W m-2 sr-1 um-1]
print getUnitsOrElse(a[0], rangeOnly=True)
# output = [W m-2 sr-1 um-1]

Note that for the last example, where the data was loaded in via loadGrid, I had to use the rangeOnly=True keyword to get only the 'W m-2 sr-1 um-1' range units to print out. Without this keyword, the domain unit [rad, rad] is printed out as well.

I haven't done a ton of testing with this function, but it seems to work for most of the cases in your script earlier in this post.

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

Re: How do I extract units consistently?

Post by joleenf »

Thanks Bob (and programmer). The name made me laugh. Thanks for adding a little joy into the process.
Joleen
User avatar
bobc
Posts: 988
Joined: Mon Nov 15, 2010 5:57 pm

Re: How do I extract units consistently?

Post by bobc »

Hi Joleen -

I meant to get this to you sooner, but the programmer came up with an improved method of extracting units from different data objects. From the programmer:
It basically simplifies things, and the function that users would use is now named "findUnits". The default behavior will now only find "range" units, but you can find every unit via "findUnits(data, rangeOnly=False)".

Here's the code for the findUnits() function:

Code: Select all

from visad import Data
from visad import FunctionType
from visad import RealType
from visad import TupleType

def _findUnits(mathtype, foundUnits, rangeOnly, withinRange):
    if isinstance(mathtype, FunctionType):
        domain = mathtype.getDomain()
        dims = domain.getDimension()
        for dim in range(dims):
            comp = domain.getComponent(dim)
            _findUnits(comp, foundUnits, rangeOnly, False)
        rng = mathtype.getRange()
        _findUnits(rng, foundUnits, rangeOnly, True)
    elif isinstance(mathtype, TupleType):
        dims = mathtype.getDimension()
        for dim in range(dims):
            comp = mathtype.getComponent(dim)
            _findUnits(comp, foundUnits, rangeOnly, withinRange)
    elif isinstance(mathtype, RealType):
        unit = mathtype.getDefaultUnit()
        if unit is not None:
            if rangeOnly and withinRange:
                foundUnits.append(unit)
            elif not rangeOnly:
                foundUnits.append(unit)

def findUnits(d, rangeOnly=True):
    if not isinstance(d, Data):
        raise TypeError("Unknown type: {0}".format(type(d)))
    units = []
    if rangeOnly:
        _findUnits(d.getType(), units, True, False)
    else:
        _findUnits(d.getType(), units, False, False)
    return units

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

Re: How do I extract units consistently?

Post by joleenf »

Hi,

Will findUnits be incorporated into the 1.7 release, or will I need to add this to my jython library for the 1.7 official release?

Joleen
User avatar
Jon
Posts: 192
Joined: Fri Jan 09, 2009 8:44 pm
Location: Madison, WI

Re: How do I extract units consistently?

Post by Jon »

I will be sure to make that infernal programmer get this stuff into the 1.7 release!
User avatar
joleenf
Posts: 1123
Joined: Mon Jan 19, 2009 7:16 pm

Re: How do I extract units consistently?

Post by joleenf »

@Jon - your programming, including this little function, has made my work much easier.

Joleen
Post Reply