Open durack1 opened 10 years ago
You folks might want to assign this as a 2.1 milestone to get if off the to-do list
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 ;)
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?
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]])
@doutriaux1 another 2.3 candidate..
@dnadeau4 @doutriaux1 just pinging you both on this..
@durack1 Can you check issue_402 and tell me if it works for you? I added ESMF for non uniform grids as you suggest.
@dnadeau4 sure, will do this week.. What have you changed in the branch?
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())
@dnadeau4 ok, so changes are in genutil
and specifically the area_weights
function as outlined above?
@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?
This example outlines the issue: