diff --git a/da/ct/tools.py b/da/ct/tools.py index 44f6fda4d28855970ecf48c622d2957b83a06970..dab211de71d4db521c226792e47030f1ddf7f868 100755 --- a/da/ct/tools.py +++ b/da/ct/tools.py @@ -86,23 +86,63 @@ class DaInfo(object): return None ################### End Class DaInfo ################### -def StateToGrid(StateVector): +def StateToGrid(DaInfo,values,reverse=False,avg=False): """ - This method converts parameters in a CarbonTracker StateVector object to a gridded map of linear multiplication values. These - can subsequently be used in the transport model code to multiply/manipulare fluxes + This method converts parameters from a CarbonTracker StateVector object to a gridded map of linear multiplication values. These + can subsequently be used in the transport model code to multiply/manipulate fluxes """ + import Nio + mapfile = os.path.join(DaInfo.da_settings['datadir'],DaInfo.da_settings['regionsfile']) + ncf = Nio.open_file(mapfile,'r') + regionmap = ncf.variables['budget_region'].get_value() + dummy = ncf.close() + nregions = regionmap.max() + + # dictionary for region <-> map conversions + + regs={} + for r in arange(1,nregions+1): + sel=(regionmap.flat == r).nonzero() + if len(sel[0])>0: regs[r]=sel + + regionselect=regs + + if reverse: + + """ project 1x1 degree map onto ecoregions """ + + result=zeros(nregions,float) + for k,v in regionselect.iteritems(): + if avg: + result[k-1]=values.ravel().take(v).mean() + else : + result[k-1]=values.ravel().take(v).sum() + return result + + else: + + """ project ecoregion properties onto 1x1 degree map """ + + result=zeros((180,360,),float) + for k,v in regionselect.iteritems(): + result.put(v,values[k-1]) + + return result if __name__ == "__main__": + sys.path.append('../../') + import os import sys from da.tools.general import StartLogger from da.tools.initexit import CycleControl import numpy as np import da.tools.rc as rc + from pylab import * opts = ['-v'] args = {'rc':'da.rc','logfile':'da_initexit.log','jobrcfilename':'test.rc'} @@ -112,6 +152,28 @@ if __name__ == "__main__": DaCycle.Initialize() print DaCycle + + a=arange(240)+100 + + b = StateToGrid(DaCycle.DaSystem,a) + + figure() + imshow(b) + colorbar() + print b.max() + + c = StateToGrid(DaCycle.DaSystem,b,reverse=True,avg=True) + + figure() + print c.max() + + plot(a,label='original') + plot(c,label='reconstructed') + legend(loc=0) + + + show() +