source: CMIP6dreq/trunk/src/scope.py @ 390

Subversion URL: http://proj.badc.rl.ac.uk/svn/exarch/CMIP6dreq/trunk/src/scope.py@390
Revision 390, 8.6 KB checked in by mjuckes, 5 years ago (diff)

For 2nd distribution; example volume estimate added

Line 
1import dreq
2import collections, string
3
4class baseException(Exception):
5  """Basic exception for general use in code."""
6
7  def __init__(self,msg):
8    self.msg = 'scope:: %s' % msg
9
10  def __str__(self):
11    return repr( self.msg )
12
13  def __repr__(self):
14    return self.msg
15
16nt_mcfg = collections.namedtuple( 'mcfg', ['nho','nlo','nha','nla','nlas','nls','nh1'] )
17
18npy = {'daily':365, u'Annual':1, u'fx':0.01, u'1hr':24*365, u'3hr':8*365, u'monClim':12, u'Timestep':100, u'6hr':4*365, u'day':365, u'1day':365, u'mon':12, u'yr':1, u'1mon':12, 'month':12, 'year':1, 'monthly':12, 'hr':24*365, 'other':24*365, 'subhr':24*365, 'Day':365, '6h':4*365,
19'3 hourly':8*365, '':1 }
20## There are 4 cmor variables with blank frequency ....
21
22def vol01( sz, v, npy, freq, inx ):
23  n1 = npy[freq]
24  s = sz[inx.uid[v][1].stid]
25  assert type(s) == type(1), 'Non-integer size found for %s' % v
26  assert type(n1) in (type(1),type(0.)), 'Non-number "npy" found for %s, [%s]' % (v,freq)
27  return s*n1
28
29class col_list(object):
30  def __init__(self):
31    self.a = collections.defaultdict(list)
32
33class col_count(object):
34  def __init__(self):
35    self.a = collections.defaultdict(int)
36
37class dreqQuery(object):
38  def __init__(self,dq=None):
39    if dq == None:
40      self.dq = dreq.loadDreq()
41    else:
42      self.dq=dq
43    self.rlu = {}
44    for i in self.dq.coll['objective'].items:
45      k = '%s.%s' % (i.mip,i.label)
46      assert not self.rlu.has_key(k), 'Duplicate label in objectives: %s' % k
47      self.rlu[k] = i.uid
48
49    self.mips = { i.mip for i in  self.dq.coll['requestItem'].items}
50    self.mipls = sorted( list( self.mips ) )
51
52    self.default_mcfg = nt_mcfg._make( [259200,60,64800,40,20,5,100] )
53    self.mcfg = {}
54    for k in self.default_mcfg.__dict__.keys():
55      self.mcfg[k] = self.default_mcfg.__dict__[k]
56    self.szcfg()
57
58  def szcfg(self):
59    self.szss = {}
60    self.sz = {}
61    for i in self.dq.coll['spatialShape'].items:
62      type = 'a'
63      if i.levelFlag == 'false':
64        ds =  string.split( i.dimensions, '|' )
65        if ds[-1] in ['site', 'basin']:
66          vd = ds[-2]
67        else:
68          vd = ds[-1]
69 
70        if vd[:4] == 'olev' or vd == 'rho':
71          type = 'o'
72          nz = self.mcfg['nlo']
73        elif vd[:4] == 'alev':
74          nz = self.mcfg['nla']
75        elif vd in ['slevel','sdepth']:
76          nz = self.mcfg['nls']
77        elif vd == 'aslevel':
78          nz = self.mcfg['nlas']
79        else:
80          print 'Failed to parse dimensions %s' % i.dimensions
81          raise
82      else:
83        nz = i.levels
84
85      dims = set( string.split( i.dimensions, '|' ) )
86      if 'latitude' in dims and 'longitude' in dims:
87        if type == 'o':
88          nh = self.mcfg['nho']
89        else:
90          nh = self.mcfg['nha']
91      else:
92        nh = 10
93
94      self.szss[i.uid] = nh*nz
95    for i in self.dq.coll['structure'].items:
96      s = self.szss[i.spid]
97      if i.odims != '':
98        s = s*5
99      self.sz[i.uid] = s
100
101  def getRequestLinkByMip( self, mipSel ):
102    """Return the set of request links which are associated with specified MIP"""
103    if type(mipSel) == type(''):
104      t1 = lambda x: x == mipSel
105    elif type(mipSel) == type({1,2}):
106      t1 = lambda x: x in mipSel
107    self.rqs = list({self.dq.inx.uid[i.rid][1] for i in self.dq.coll['objectiveLink'].items if t1(i.label) })
108    return self.rqs
109
110  def getRequestLinkByObjective( self, objSel ):
111    """Return the set of request links which are associated with specified objectives"""
112    if type(objSel) == type(''):
113      t1 = lambda x: x == self.rlu[objSel]
114    elif type(objSel) == type({1,2}):
115      t1 = lambda x: x in {self.rlu[i] for i in objSel}
116
117    self.rqs = list({self.dq.inx.uid[i.rid][1] for i in self.dq.coll['objectiveLink'].items if t1(i.oid) })
118    return self.rqs
119
120  def varGroupXexpt(self, rqList ):
121    """For a list of request links, return a list of variable group IDs for each experiment"""
122    self.cc = collections.defaultdict( list )
123    dummy = {self.cc[i.expt].append(i.rlid) for i in self.dq.coll['requestItem'].items if i.rlid in {j.uid for j in rqList} }
124    return self.cc
125
126  def yearsInRequest(self, rql ):
127    self.ntot = sum( [i.ny for i in self.dq.coll['requestItem'].items if i.rlid == rql.uid] )
128    return self.ntot
129
130  def volByExpt( self, l1, ex, pmax=2, cc=None ):
131    """volByExpt: calculates the total data volume associated with an experiment/experiment group and a list of request items.
132          The calculation has some approximations concerning the number of years in each experiment group."""
133##
134## cc: an optional collector, to accumulate indexed volumes
135##
136    inx = self.dq.inx
137    rql = {i.rlid for i in l1 if i.expt == ex}
138
139## The complete set of variables associated with these requests:
140    rqvg = list({inx.uid[i][1].refid for i in rql})
141
142###To obtain a set of variables associated with this collection of variable groups:
143    col1 = set()
144    x = {tuple( {col1.add(i) for i in inx.iref_by_sect[vg].a['requestVar']} ) for vg in rqvg}
145###The collector col1 here accumulates all the record uids, resulting in a single collection. These are request variables, to get a set of CMOR variables at priority <= pmax:
146    vars = {inx.uid[l][1].vid for l in list(col1) if inx.uid[l][1].priority <= pmax}
147
148### filter out cases where the request does not point to a CMOR variable.
149    vars = {vid for vid in vars if inx.uid[vid][0] == u'CMORvar'}
150
151    e = {}
152    for u in rql:
153### for request variables which reference the variable group attached to the link, add the associate CMOR variables, subject to priority
154      i = inx.uid[u][1]
155      e[i.uid] = { inx.uid[x][1].vid for x in inx.iref_by_sect[i.refid].a['requestVar'] if inx.uid[x][1].priority <= pmax}
156
157#
158# for each variable, calculate the maximum number of years across all the request links which reference that variable.
159#
160    nym = {}
161    for v in vars:
162      ### for each request item, check if v is in the set of variables and then add the number of years.
163      nym[v] = max( {i.ny for i in l1 if i.expt == ex and v in e[i.rlid]} )
164
165    szv = {}
166    ov = []
167    for v in vars:
168      szv[v] = self.sz[inx.uid[v][1].stid]*npy[inx.uid[v][1].frequency]
169      ov.append( self.dq.inx.uid[v][1] )
170    ee = self.listIndexDual( ov, 'frequency', 'mipTable', acount=None, alist=None, cdict=szv, cc=cc )
171    self.ngptot = sum( [  self.sz[inx.uid[v][1].stid]* npy[inx.uid[v][1].frequency] *nym[v]  for v in vars] )
172    return (self.ngptot, ee )
173
174  def summaryByMip( self, pmax=1 ):
175    bytesPerFloat = 2.
176    for m in self.mipls:
177      v = self.volByMip( m, pmax=pmax )
178      print '%12.12s: %6.2fTb' % (m,v*bytesPerFloat*1.e-12)
179
180  def volByMip( self, mip, pmax=2):
181
182    if type(mip) in {type( '' ),type( u'') }:
183      if mip not in self.mips:
184        print self.mips
185        raise baseException( 'volByMip: Name of mip not recognised: %s' % mip )
186      l1 = [i for i in  self.dq.coll['requestItem'].items if i.mip == mip]
187    elif type(mip) == type( set()):
188      nf = [ m for m in mip if m not in self.mips]
189      if len(nf) > 0:
190        raise baseException( 'volByMip: Name of mip(s) not recognised: %s' % str(nf) )
191      l1 = [i for i in  self.dq.coll['requestItem'].items if i.mip in mip]
192    else:
193      raise baseException( 'volByMip: "mip" (1st explicit argument) should be type string or set: %s -- %s' % (mip, type(mip))   )
194     
195    #### The set of experiments/experiment groups:
196    exps = {i.expt for i in l1}
197    self.volByE = {}
198    vtot = 0
199    cc = collections.defaultdict( col_count )
200    for e in exps:
201      self.volByE[e] = self.volByExpt( l1, e, pmax=pmax, cc=cc )
202      vtot += self.volByE[e][0]
203    self.indexedVol = cc
204
205    return vtot
206
207  def listIndexDual(self, ll, a1, a2, acount=None, alist=None, cdict=None, cc=None ):
208    do_count = acount != None
209    do_list = alist != None
210    assert not (do_count and do_list), 'It is an error to request both list and count'
211    if not (do_count or do_list):
212      acount = '__number__'
213      do_count = True
214
215    if cc == None:
216      if do_count:
217        cc = collections.defaultdict( col_count )
218      elif do_list:
219        cc = collections.defaultdict( col_list )
220
221    if do_count:
222      for l in ll:
223        if cdict != None:
224          v = cdict[l.uid]
225        elif acount == '__number__':
226          v = 1
227        else:
228          v = l.__dict__[acount]
229
230        cc[ l.__dict__[a1] ].a[ l.__dict__[a2] ] += v
231    elif do_list:
232      for l in ll:
233        if cdict != None:
234          v = cdict[l.uid]
235        elif alist == '__item__':
236          v = l
237        else:
238          v = l.__dict__[alist]
239        cc[ l.__dict__[a1] ].a[ l.__dict__[a2] ].append( v )
240
241    od = {}
242    for k in cc.keys():
243      d2 = {}
244      for k2 in cc[k].a.keys():
245        d2[k2] = cc[k].a[k2]
246      od[k] = d2
247    return od
Note: See TracBrowser for help on using the repository browser.