source: CedaCdoScripts/branch/c1/ncExtras.py @ 126

Subversion URL: http://proj.badc.rl.ac.uk/svn/exarch/CedaCdoScripts/branch/c1/ncExtras.py@126
Revision 126, 3.8 KB checked in by mjuckes, 6 years ago (diff)

added files

Line 
1
2## script to create a netcdef file with variables to be inserted in another by ncks.
3##
4## classes of variables:
5## time: time axis to overwrite existing axis messed up by CDO
6## climatology:  climatological bounds, typically derived from time dimension of a reference netcdf file
7## scalar: a scalar variable, e.g to specify a threshold.
8##
9
10## usage [--cb name] [--sc name,value,a1=val;a2=val;a3=val] [--time calendar,units,len,inc,incUnits] <ref netcdf> <out netcdf>
11
12import cdms2 as cdms
13import os, string
14import numpy
15
16def leny( y ):
17     if (y/4)*4 != y:
18       return 365
19     if (y/400)*400 == y:
20       return 366
21     if (y/100)*100 == y:
22       return 365
23     return 366
24
25
26def run(argsIn):
27
28  args = argsIn[:]
29  print args
30  climBounds=False
31  scalarVar = False
32  timeAxis = False
33  while len(args) > 2:
34    if args[0] == '--cb':
35      climBounds=True
36      cbName = args[1]
37      args.pop(0)
38    elif args[0] == '--sc':
39      scalarVar = True
40      bits = string.split( args[1], ',' )
41      svName = bits[0]
42      svVal = float(bits[1])
43      svAts = {}
44      if len(bits) > 2:
45        b2 = string.split( bits[2], ';' )
46        for b in b2:
47          a,v = string.split(b, '=' )
48          svAts[a] = v
49      args.pop(0)
50
51    elif args[0] == '--time':
52      timeAxis = True
53      bits = string.split( args[1], ',' )
54      cal = bits[0]
55      units = string.strip(bits[1], '"' )
56      ll = int( bits[2] )
57      inc = int( bits[3] )
58      incUnits = bits[4]
59      print cal,units,ll,inc,incUnits
60      assert ll < 10000, 'to many bits, %s' % ll
61      args.pop(0)
62    args.pop(0)
63 
64  if not (scalarVar or climBounds or timeAxis):
65    print 'Nothing to do', argsIn
66    sys.exit(0)
67 
68  if climBounds and not timeAxis:
69    nc = cdms.open( args[0] )
70    t = nc.getAxis('time')
71    tvals = t.getValue()
72    tub = tvals[-1] + [tvals[-1]-tvals[0]]/(len(tvals)-1)
73    ll = len(tvals)
74    nc.close()
75  elif timeAxis and ll == -1:
76    nc = cdms.open( args[0] )
77    t = nc.getAxis('time')
78    tv = t.getValue()
79    ll = len(tv)
80    nc.close()
81 
82  if timeAxis:
83    tu = string.split( units )[0]
84    print tu, incUnits
85    if tu == incUnits:
86      tvals = range(ll)
87      tub = ll
88    elif tu == 'days' and incUnits == 'years':
89      if cal in ['365-day','noleap','366-day','all-leap','360-day']:
90        if cal in ['365-day','noleap']:
91          ylen = 365
92        elif cal in ['366-day','all-leap']:
93          ylen = 366
94        else:
95          ylen = 360
96        tvals = map( lambda x: x*ylen, range(ll) )
97        tub = tvals[-1] + ylen
98      else:
99        bu = string.split( units )
100        ystart = int(bu[2][:4] )
101        tvals = [0,]
102        for i in range(ll-1):
103          y = ystart + i
104          tvals.append( tvals[-1] + leny( y ) )
105        tub = tvals[-1] + leny( ystart + ll -1 )
106    else:
107      print 'options not supported'
108      sys.exit(0)
109 
110  if timeAxis or climBounds:
111    tvalsd = numpy.array( tvals, dtype = 'double' )
112 
113  nc = cdms.open( args[1], 'w' )
114  if timeAxis:
115    t = nc.createAxis( 'time', tvalsd, unlimited=True )
116    t.standard_name = 'time'
117    t.long_name = 'time'
118    t.units = units
119    if climBounds:
120      t.climatology = cbName
121 
122  if climBounds:
123    if not timeAxis:
124      t = nc.createAxis( 'time', tvalsd, unlimited=True )
125    ba = numpy.array( [0,1], dtype = 'int' )
126    b = nc.createAxis( 'nb2', ba )
127    cb = nc.createVariable( cbName, 'double', [t,b] )
128   
129    cbd1 = tvals[:]
130    cbd2 = tvals[1:]
131    cbd2.append( tub )
132    cbda = map( lambda x,y: [x,y], cbd1, cbd2 )
133    cbd = numpy.array(  cbda, 'double' )
134    cb.assignValue( cbd )
135    nc.write( cb )
136
137  if scalarVar:
138    #sv = nc.createVariable( svName, 'double', [] )
139    #sv.assignValue( svVal )
140    print svAts
141    nc.write( svVal, id=svName, attributes=svAts )
142   
143  nc.close()
144
145if __name__ == "__main__":
146  import sys
147  args = sys.argv[1:]
148  run(args)
Note: See TracBrowser for help on using the repository browser.