source: TI02-CSML/trunk/parser/API/ops_GridSeriesFeature.py @ 1093

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/TI02-CSML/trunk/parser/API/ops_GridSeriesFeature.py@1093
Revision 1093, 6.4 KB checked in by domlowe, 14 years ago (diff)

Can set file paths of output xml and netcdf. Added example to apicalls.py

Line 
1''' ops_GridSeriesFeature  contains operations for GridSeriesFeatures'''
2from API import *
3from CSMLDocument import *
4from Scientific.IO.NetCDF import *  #use this instead of cdms for now for it's simple write interface..
5from Numeric import *
6
7
8def testmethod(self):
9    print 'testmethod for gridseries feature'
10    return 'testmethod - gridseries'
11
12def getAllowedSubsettings(self):
13    return ['subsetToGridSeries']
14
15def getDomainReference(self):
16    #Inspects a time position list for the domain reference.
17    #TODO, does not handle a file extract in place of a list.
18    if isinstance(self.domain.domainReference,Parser.TimePositionList):
19        time = {}
20        time['t'] = self.domain.domainReference.timePositions
21        domainref  = time
22    return domainref
23   
24def getDomainComplement(self):
25    #This will return a list containing one or more ordinates:
26    #currently in form [Name, values]
27    #assumes ordinate.axisValues is a file extract
28    #TODO axisValues may be inline -see also domainReference for similar problem
29    domaincomp ={}
30    dc = self.domain.domainComplement
31    #dc should be a grid!
32    if isinstance(dc, Parser.Grid):
33        for ordinate in dc.ordinates:
34            domaincomp[ordinate.definesAxis]=ordinate.axisValues.getData()
35    return domaincomp
36
37def getDomain(self):
38    #returns both the domain reference axes and domain compliment axes in a single domain dictionary
39    #axes are in no particular order
40    domain = {}
41    dr=getDomainReference(self)
42    dc=getDomainComplement(self)
43    for key in dc.keys():
44        domain[key]=dc[key]
45    for key in dr.keys():
46        domain[key]=dr[key]
47    return domain
48
49
50       
51
52def subsetToGridSeries(self, timeSubset,  csmlpath=None, ncpath=None,**kwargs):
53    if csmlpath is not None:
54        pathToSubsetCSML = csmlpath
55    else:
56        pathToSubsetCSML='temp.xml'
57    if ncpath is not None:
58        pathToSubsetNetCDF=ncpath
59    else:
60        pathToSubsetCSML='temp.nc'
61
62    domainref = getDomainReference(self) 
63    self.times=timeSubset
64    self.files=[]
65    strTimes=''
66    fulldata=[]
67    if len(self.times) == 2:
68        tone=ops_AbstractFeature.__getCDtime(self.times[0])
69        ttwo=ops_AbstractFeature.__getCDtime(self.times[1])
70        dr=getDomainReference(self)
71        self.times=[]
72        for time in dr['t'].split():
73            timeok=ops_AbstractFeature.__compareTimes(tone,time,ttwo)
74            if timeok ==1:
75                self.times.append(time)
76    for time in self.times:
77        listPosition=domainref['t'].split().index(time)
78        strTimes= strTimes + ' ' + time
79        for comp in self.rangeSet.aggregatedArray.component:
80            data=comp.getData(fileposition=listPosition, **kwargs)
81            self.files.append(comp.fileName.split()[listPosition])
82            if fulldata ==[]:
83                fulldata = data.tolist()
84            else:
85                for item in data.tolist():
86                    fulldata.append(item)
87     #get the calendar type
88    try:
89        caltype, calunits = ops_AbstractFileExtract.__calendar(self.rangeSet.aggregatedArray.component[0].fileName.split()[0], 't')    #TODO should accept any time dim!!
90        csmltime.setcdtimeCalendar(caltype)
91    except:
92        caltype=cdtime.DefaultCalendar   
93    ### define domain and rangeSet to use for feature in csml document####
94    domain=Parser.GridSeriesDomain()
95    domain.domainReference=Parser.TimePositionList(timePositions=strTimes) 
96    grid=Parser.Grid()
97    dc = self.getDomainComplement()
98    ordinates= []
99    i=0
100    valueStore=[]  # use the values again later to generate netcdf
101    arraySize=0
102    totalArraySize=1
103    for key in dc.keys():
104        arraySize=0
105        i=i+1
106        god=Parser.GridOrdinateDescription()
107        god.gridAxesSpanned='dim%s'%i
108        god.sequenceRule='+x+y+z'
109        god.definesAxis=key
110        straxisValues=''
111        if key in kwargs:
112            for val in dc[key]:
113                if val >= kwargs[key][0]:
114                    if val <= kwargs[key] [1]:
115                        arraySize=arraySize+1
116                        straxisValues=straxisValues+ str(val) + ', '
117        else: # this dimension has not been subsetted
118            for val in dc[key]:
119                arraySize=arraySize+1
120                straxisValues=straxisValues+ str(val) + ', '
121        totalArraySize=totalArraySize*arraySize
122        god.axisValues=straxisValues[:-2]
123        ordinates.append(god)
124    totalArraySize=totalArraySize*len(self.times)
125    grid.ordinates=ordinates
126    domain.domainComplement=grid
127    rangeSet=Parser.RangeSet()
128    rangeSet.arrayDescriptor=Parser.NetCDFExtract(id=self.id,fileName='temp.nc',variableName=self.id,arraySize=[arraySize])
129
130    #### write csml document ##### -move this to the csmldocument module?
131    subsetCSML=CSMLDocument()
132    subsetCSML=subsetCSML.makeGridSeries(domain,rangeSet)
133    output=open(pathToSubsetCSML,'w')
134    output.write(subsetCSML)
135    output.close()
136
137    ### create and write netcdf - uses scientific python####
138    ncfile=NetCDFFile(pathToSubsetNetCDF,'w')
139    # create the dimensions       
140    ncfile.createDimension ( 'time', len(self.times))
141    time_var = ncfile.createVariable ( 'time', Float, ('time',) )
142    time_var.longname = 'time'
143    floatTimes=[]
144    print len(fulldata[0])
145    for time in self.times:
146        time=ops_AbstractFeature.__getCDtime(time).torel(calunits)
147        floatTimes.append(time.value)
148    time_var[:] =floatTimes[:]
149    time_var.units=calunits
150    time_var.calendar=caltype
151
152    for ordinate in ordinates:
153        ncfile.createDimension(ordinate.definesAxis, len(ordinate.axisValues.split())) 
154        item_var = ncfile.createVariable (ordinate.definesAxis, Float, (ordinate.definesAxis,) )
155        #convert to list
156        vals=[]
157        for val in ordinate.axisValues.split(','):
158            vals.append(float(val))
159        ordinate.axisValues=vals
160        item_var[:]=vals[:]
161        print ordinate.definesAxis
162    #this needs reconsidering - do the shapes always match up??
163    if len(ordinates)==3:
164        feature_var = ncfile.createVariable (self.id, Float, ('time',ordinates[0].definesAxis,ordinates[1].definesAxis,ordinates[2].definesAxis))
165    elif len(grid.ordinates)==2:
166        feature_var = ncfile.createVariable (self.id, Float, ('time',ordinates[0].definesAxis,ordinates[1].definesAxis))
167    print shape(feature_var)
168    print shape(fulldata)
169    feature_var[:]=fulldata[:]
170    ncfile.close()
171
172    return pathToSubsetCSML, pathToSubsetNetCDF, totalArraySize
Note: See TracBrowser for help on using the repository browser.