2

I'm new to the Iris library and need to compute a heat index, which is a multivariate non-linear function of temperature and relative humidity that goes something like $HI = temp +rh + temp*rh + temp^2 *rh + rh^2*temp $ . Here, temp has units of Fahrenheit and rh has units of 1.

However, Iris cubes won't add differing units:

In [147]: HI = temp + rh + temp*rh + temp**2*rh + rh**2*temp
---------------------------------------------------------------------------
NotYetImplementedError                    Traceback (most recent call last)
<ipython-input-147-675ea72a5d06> in <module>()
----> 1 HI = temp + rh + temp*rh + temp**2*rh + rh**2*temp

/Users/me/anaconda/lib/python2.7/site-packages/iris/cube.pyc in __add__(self, other)
  2596 
   2597     def __add__(self, other):
-> 2598         return iris.analysis.maths.add(self, other, ignore=True)
   2599     __radd__ = __add__
   2600 

/Users/me/anaconda/lib/python2.7/site-packages/iris/analysis/maths.pyc in add(cube, other, dim, ignore, in_place)
    166     op = operator.iadd if in_place else operator.add
167     return _add_subtract_common(op, 'addition', 'added', cube, other, dim=dim,
--> 168                                 ignore=ignore, in_place=in_place)
169 
170 

/Users/me/anaconda/lib/python2.7/site-packages/iris/analysis/maths.pyc     in _add_subtract_common(operation_function, operation_noun,     operation_past_tense, cube, other, dim, ignore, in_place)
    216     """
    217     _assert_is_cube(cube)
--> 218     _assert_matching_units(cube, other, operation_noun)
219 
220     if isinstance(other, iris.cube.Cube):

/Users/me/anaconda/lib/python2.7/site-packages/iris/analysis/maths.pyc in _assert_matching_units(cube, other, operation_noun)
132         raise iris.exceptions.NotYetImplementedError(
133             'Differing units (%s & %s) %s not implemented' %
--> 134             (cube.units, other.units, operation_noun))
135 
136 

NotYetImplementedError: Differing units (Fahrenheit & 1) addition not implemented

If I call the data instead as a numpy array then this is a workaround, e.g.: heatIndex = -42.379 + temp.data + rh.data + temp.data2 + rh.data2 but this seems to defeat the purpose of using Iris in the first place, and requires re-writing the meta-data.

Is this possible to do with iris cubes? Is there a unit-less udunit that I'm missing that would allow this to happen?

NB: if it wasn't clear from the error, I'm running Python 2.7 (and Iris 1.7).

Anna S.
  • 153
  • 2
  • 7

2 Answers2

0

If I call the data instead as a numpy array ... this seems to defeat the purpose of using Iris in the first place, and requires re-writing the meta-data.

Firstly, making use of numpy and all of the incredible scipy ecosystem is a powerful capability of Iris. Where functionality exists in Iris, it is better to use it rather than having to manage the metadata yourself, but if something doesn't exist then hopefully it shouldn't be too difficult to implement.

In your case, you will probably be wanting to update some of the metadata yourself (i.e. Iris wont know that this mathematical transformation produces data which represents "heat index").

It appears the heat index should be in Fahrenheit, and the addition of the Fahrenheit cube to a scalar cube is causing this error. The simplest solution might be to just add the relative humidity scalar data (i.e. not the cube) to the temperature:

rh = rh.data
heat_index = temp + rh + temp*rh + temp**2*rh + rh**2*temp
heat_index.rename('heat_index')
assert heat_index.units == 'Fahrenheit'
pelson
  • 21,252
  • 4
  • 92
  • 99
  • This is putting me on the right track, though it doesn't really solve the issue of adding temp + temp**2 (heat index units are degrees F + degrees F^2). I got around this by factoring and multiplying as: temp*(a+b*temp) + c. It's taking a horrendously long time to calculate, so hopefully this will be a workable situation in the future. – Anna S. Feb 03 '15 at 22:49
  • Edit: this solution causes python to run out of memory. I'm just sticking with a Numpy array for this problem. – Anna S. Feb 04 '15 at 02:54
0

I think that you were on the right lines when you started working with the data arrays, but it can be made a bit cleaner. The key to this is to start by copying the temperature cube and work from there

C1 = -42.379
C2 = 2.04901523
# etc

def get_heat_index(temp, rh):
    '''Calculate the heat index as defined by George Winterling'''

    # make sure we are using the right temperature units
    # other checks could be implemented as assertions
    temp.convert_units('Fahrenheit')

    # start by copying the temperature cube
    heat_index = temp

    # fix the name of the cube to be returned
    heat_index.rename('heat_index')

    # C2 - do this first as already have temperature data
    heat_index.data = (C1 + C2*temp.data +
                      C3*rh.data +
    # etc
                      C9*temp.data**2*rh.data**2)

    return heat_index