genutil.area_weights doesn't work with non-cartesian grids (most CMIP5 ocean grids)

[durack1]

This example outlines the issue:

import cdms2 as cdm
import genutil as gutil
infile = '/work/cmip5/historical/ocn/mo/thetao/cmip5.ACCESS1-0.historical.r1i1p1.mo.ocn.Omon.thetao.ver-1.latestX.xml'
f_h = cdm.open(infile)
mask = f_h('thetao',time=slice(0,1),lev=slice(0,1))(squeeze=1)
ocean_area_frac = gutil.area_weights(mask)

>>> ocean_area_frac = gutil.area_weights(mask)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/uvcdat/latest/lib/python2.7/site-packages/genutil/averager.py", line 188, in area_weights
    wt = __myGetAxisWeights(ds, 0, axisoptions)
  File "/usr/local/uvcdat/latest/lib/python2.7/site-packages/genutil/averager.py", line 54, in __myGetAxisWeights
    raise AveragerError, 'Bounds not available to compute weights for dimension: '+ax.id
genutil.averager.AveragerError: ('B', 'o', 'u', 'n', 'd', 's', ' ', 'n', 'o', 't', ' ', 'a', 'v', 'a', 'i', 'l', 'a', 'b', 'l', 'e', ' ', 't', 'o', ' ', 'c', 'o', 'm', 'p', 'u', 't', 'e', ' ', 'w', 'e', 'i', 'g', 'h', 't', 's', ' ', 'f', 'o', 'r', ' ', 'd', 'i', 'm', 'e', 'n', 's', 'i', 'o', 'n', ':', ' ', 'j')

[durack1]

You folks might want to assign this as a 2.1 milestone to get if off the to-do list

[doutriaux1]

is is a regular lat/lon grid? area_weight only works on this type of grids, I guess we could add a check and exit in a more civilized maner if that's the case ;)

[durack1]

Nope the issue was with the non lat/lon grids.. Most of the CMIP5 oceans.. A better way would be to use the ESMF tools to deal with non lat/lon grids instead.. Is this technically possible?

[durack1]

I do believe it would be possible to generate the weights using a call to ESMF, so like this example:

import cdms2 as cdm
infile = '/work/durack1/Shared/cmip5/historical/ocn/an/tos/cmip5.GFDL-CM3.historical.r1i1p1.an.ocn.Omon.tos.ver-v1.latestX.1860-2005.nc'
tos = cdm.open(infile)('tos')[0,...]
diag = {'srcAreaFractions':None, 'srcAreas':None, 'dstAreas':None} ; # Create dictionary, with additional keys for function to return
None = tos.regrid(tos.getGrid(),regridTool='esmf',regridMethod='conserve',diag=diag)
srcAreas = diag.get('srcAreas')
srcAreas
array([[  4.50225393e-05,   4.50225393e-05,   4.50225393e-05, ...,
          4.50225393e-05,   4.50225393e-05,   4.50225393e-05],
       [  5.02732898e-05,   5.02732898e-05,   5.02732898e-05, ...,
          5.02732898e-05,   5.02732898e-05,   5.02732898e-05],
       [  5.55087334e-05,   5.55087334e-05,   5.55087334e-05, ...,
          5.55087334e-05,   5.55087334e-05,   5.55087334e-05],
       ..., 
       [  2.07087362e-06,   6.05800431e-06,   9.80637603e-06, ...,
          9.80561961e-06,   6.05782234e-06,   2.07127099e-06],
       [  2.08306596e-06,   6.08933391e-06,   9.85335671e-06, ...,
          9.85375251e-06,   6.08972623e-06,   2.08286389e-06],
       [  2.08856097e-06,   6.10596819e-06,   9.87812170e-06, ...,
          9.87831337e-06,   6.10577112e-06,   2.08856097e-06]])

[durack1]

@doutriaux1 another 2.3 candidate..

[durack1]

@dnadeau4 @doutriaux1 just pinging you both on this..

[dnadeau4]

@durack1 Can you check issue_402 and tell me if it works for you? I added ESMF for non uniform grids as you suggest.

[durack1]

@dnadeau4 sure, will do this week.. What have you changed in the branch?

[dnadeau4]

I actually added your code to call ESMF regridder in the averager and return the srcAreas weight when the grid is non uniform. No use of Numpy for non uniform grid.

The branch is called issue_402.

. diag = {'srcAreaFractions':None, 'srcAreas':None, 'dstAreas':None} ; # Create dictionary, with additional keys for function to return
. dummy=ds.regrid(ds.getGrid(),regridTool='esmf',regridMethod='conserve',diag=diag)
. wt = numpy.array(diag.get('srcAreas'))

. return cdms2.createVariable(numpy.ma.masked_array(wt, numpy.ma.getmask(ds)), axes=ds.getAxisList())

[durack1]

@dnadeau4 ok, so changes are in genutil and specifically the area_weights function as outlined above?

[durack1]

@dnadeau4 ok so it seems to work and output reasonable fields - though I agree with you, we need to return srcAreaFractions rather than the current srcArea ESMF output:

[durack1@oceanonly ~]$ ipython
Python 2.7.10 (default, Sep 16 2015, 12:49:34)
Type "copyright", "credits" or "license" for more information.

IPython 3.0.0 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

In [1]: import cdms2 as cdm
In [2]: import genutil as gutil
In [3]: infile = '/work/cmip5/historical/ocn/mo/thetao/cmip5.ACCESS1-0.historical.r1i1p1.mo.ocn.Omon.thetao.ver-1.latestX.xml'
In [4]: f_h = cdm.open(infile)
In [5]: mask = f_h('thetao',time=slice(0,1),lev=slice(0,1))(squeeze=1)
In [6]: ocean_area_frac = gutil.area_weights(mask)
/export/durack1/built/lib/python2.7/site-packages/ESMP/src/ESMP_API.py:1241: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  if (srcMaskValues != None):
/export/durack1/built/lib/python2.7/site-packages/ESMP/src/ESMP_API.py:1248: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  if (dstMaskValues != None):
In [7]: import vcs
In [8]: x = vcs.init()
In [9]: x.plot(ocean_area_frac)
Out[9]: <vcs.displayplot.Dp at 0x7f3bca331750>

I think we need to go one further, and tweak the output field to include bounds so the genutil.averager function can then also work.. Currently, it's not:

In [13]: ocean_area_frac.shape
Out[13]: (300, 360)
In [14]: gutil.averager(ocean_area_frac,axis=0)
---------------------------------------------------------------------------
AveragerError                             Traceback (most recent call last)
<ipython-input-14-cd98dcb1f59e> in <module>()
----> 1 gutil.averager(ocean_area_frac,axis=0)

/export/durack1/built/lib/python2.7/site-packages/genutil/averager.pyc in averager(V, axis, weights, action, returned, weight, combinewts)
   1079     if __DEBUG__: print 'Checking weights= options:',weights
   1080     #
-> 1081     filled_wtoptions = __check_weightoptions(V, axis, weights)
   1082     if __DEBUG__: print 'The weights options are ', filled_wtoptions
   1083     #

/export/durack1/built/lib/python2.7/site-packages/genutil/averager.pyc in __check_weightoptions(x, axisoptions, weightoptions)
    354             #
    355             if __DEBUG__: print 'I have only 1 axis and 1 option'
--> 356             weightoptions = __check_each_weight_option(x, axislist[0], axisindex[0], weightoptions)
    357         else:
    358             #

/export/durack1/built/lib/python2.7/site-packages/genutil/averager.pyc in __check_each_weight_option(x, ax, index, wtopt)
    288             # NOTE: THIS FUNCTION CAN BE CHANGED WHEN BOB PUTS THE FUNCTION getAxisWeights() INTO cdms2
    289             if __DEBUG__: print 'GENERATE weights specified'
--> 290             wtopt = __myGetAxisWeights(x, index)
    291             if __DEBUG__: print wtopt
    292         else:

/export/durack1/built/lib/python2.7/site-packages/genutil/averager.pyc in __myGetAxisWeights(x, i, axisoptions)
     52                     axis_wts=numpy.ones(x.getAxis(i)[:].shape)
     53             else:
---> 54                 raise AveragerError, 'Bounds not available to compute weights for dimension: '+ax.id
     55         else:
     56             axis_wts =  abs(axis_bounds[..., 1] - axis_bounds[..., 0])

AveragerError: ('B', 'o', 'u', 'n', 'd', 's', ' ', 'n', 'o', 't', ' ', 'a', 'v', 'a', 'i', 'l', 'a', 'b', 'l', 'e', ' ', 't', 'o', ' ', 'c', 'o', 'm', 'p', 'u', 't', 'e', ' ', 'w', 'e', 'i', 'g', 'h', 't', 's', ' ', 'f', 'o', 'r', ' ', 'd', 'i', 'm', 'e', 'n', 's', 'i', 'o', 'n', ':', ' ', 'j')

It'd also be nice to clean up the error messages too..

@doutriaux1 what are your thoughts on this?