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

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

netcdf writing now done with cdms rather than scientific python

Line 
1''' ops_GridSeriesFeature  contains operations for GridSeriesFeatures'''
2from API import NetCDFWriter
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        pathToSubsetNetCDF='temp.nc'
61   
62    #deal with longitude requests
63    #if the request is in -ve,+ve eg (-30,30) but the data is in (0,360) need to handle this by changing the args.
64    dc = self.getDomainComplement()
65    for key in dc.keys():
66        if key == 'longitude': #how do we test if it is longitude properly?
67            print dc[key]
68            for val in dc[key]:
69                if val < 0:
70                    pass
71                else:
72                    if kwargs[key][0] < 0:
73                        kwargs[key]=(kwargs[key][0]+360,kwargs[key][1])
74                    if kwargs[key][1] < 0:
75                        kwargs[key]=(kwargs[key][0],kwargs[key][1]+360)
76         
77    domainref = getDomainReference(self)
78    self.times=timeSubset
79    self.files=[]
80    strTimes=''
81    fulldata=[]
82    if len(self.times) == 2:
83        try:
84            tone=ops_AbstractFeature.__getCDtime(self.times[0])
85        except:
86            tone=self.times[0]
87        try:
88            ttwo=ops_AbstractFeature.__getCDtime(self.times[1])
89        except:
90            ttwo=self.times[1]
91        dr=getDomainReference(self)
92        self.times=[]
93        for time in dr['t'].split():
94            timeok=ops_AbstractFeature.__compareTimes(tone,time,ttwo)
95            if timeok ==1:
96                self.times.append(time)
97   
98    #to keep track of files that have already been fetched. eg. if multiple times are in a single file only need to get data from that file once...
99    numFiles=len(self.rangeSet.aggregatedArray.component[0].fileName.split())
100    timeToFileRatio=len(domainref['t'].split())/numFiles
101    filesFetched=[]
102    #get data:
103    for time in self.times:
104        listPosition=domainref['t'].split().index(time)
105        strTimes= strTimes + ' ' + time
106        timedim=self.rangeSet.aggregatedArray.component[0].variableName
107        for comp in self.rangeSet.aggregatedArray.component:
108            filePos=(listPosition)/timeToFileRatio
109            if filePos in filesFetched:
110                continue #already got data from this file, try next time
111            data=comp.getData(fileposition=filePos, times=self.times, **kwargs)
112            self.files.append(comp.fileName.split()[filePos])
113            if fulldata ==[]:
114               fulldata = data.tolist()
115            else:
116                for item in data.tolist():
117                    fulldata.append(item)
118            filesFetched.append(filePos)
119        axisorder = data.getAxisIds()  #will need later!
120        print 'filesFetched %s'%filesFetched
121    try:
122        caltype=self.domain.domainReference.frame.split(':',1)[0]
123        calunits=self.domain.domainReference.frame.split(':',1)[1]
124        csmltime.setcdtimeCalendar(caltype)
125    except:
126        csmltime.setcdtimeCalendar(csmltime.cdtime.DefaultCalendar)
127    ### define domain and rangeSet to use for feature in csml document####
128    domain=Parser.GridSeriesDomain()
129    domain.domainReference=Parser.TimePositionList(timePositions=strTimes)
130    grid=Parser.Grid()
131    #dc = self.getDomainComplement()
132    ordinates= []
133    i=0
134    valueStore=[]  # use the values again later to generate netcdf
135    arraySize=0
136    totalArraySize=1
137    for key in dc.keys():
138        arraySize=0
139        i=i+1
140        god=Parser.GridOrdinateDescription()
141        god.gridAxesSpanned='dim%s'%i
142        god.sequenceRule='+x+y+z'
143        god.definesAxis=key
144        straxisValues=''
145       
146        #now deal with each argument:
147        if key in kwargs:
148            if kwargs[key][0] < kwargs[key][1]:   
149                for val in dc[key]:
150                    if val >= kwargs[key][0]:
151                        if val <= kwargs[key] [1]:
152                            arraySize=arraySize+1
153                            straxisValues=straxisValues+ str(val) + ', '
154            else:#this if deals with selections such as longitude (330,30) where the lower limit is 'greater' than the upper limit in a mathematical sense.
155                    for val in dc[key]:
156                        if val >= kwargs[key][0]:
157                            arraySize=arraySize+1
158                            straxisValues=straxisValues+ str(val) + ', '
159                    for val in dc[key]:
160                        if val <= kwargs[key] [1]:
161                            arraySize=arraySize+1
162                            straxisValues=straxisValues+ str(val) + ', '
163        else: # this dimension has not been subsetted at all
164            for val in dc[key]:
165                arraySize=arraySize+1
166                straxisValues=straxisValues+ str(val) + ', '
167        totalArraySize=totalArraySize*arraySize
168        god.axisValues=straxisValues[:-2]
169        ordinates.append(god)
170    totalArraySize=totalArraySize*len(self.times)
171    grid.ordinates=ordinates
172    domain.domainComplement=grid
173    rangeSet=Parser.RangeSet()
174    rangeSet.arrayDescriptor=Parser.NetCDFExtract(id=self.id,fileName=pathToSubsetNetCDF,variableName=self.id,arraySize=[arraySize])
175
176    #### write csml document #####
177    csmldoc=CSMLDocument("mydoc", "mymetadata")
178    csmldoc.addGridSeriesFeature(domain,rangeSet,datasetID="A",featureID="B",description="C")
179    csmldoc=csmldoc.consolidate()
180    output=open(pathToSubsetCSML,'w')
181    output.write(csmldoc)
182    output.close()
183
184    ### write netcdf using NCWriter class (wraps cdms) ###
185    #nc=NetCDFWriter.NCwriter(pathToSubsetNetCDF)
186    nc=NetCDFWriter.NCwriter('temp.nc')
187    floatTimes=[]
188    for time in self.times:
189        time=ops_AbstractFeature.__getCDtime(time).torel(calunits)
190        floatTimes.append(time.value)
191    nc.addAxis('t',floatTimes,isTime=1,units=calunits,calendar=caltype)
192    for ordinate in ordinates:
193        #convert to list
194        vals=[]
195        for val in ordinate.axisValues.split(','):
196            vals.append(float(val))
197        nc.addAxis(ordinate.definesAxis,vals)
198    #nc.addVariable(data,'lala', axes,history='blah blah', units='flakesperminute',long_name='snowfall_at_surface')
199    if len(ordinates)==3:
200        axes=['t',axisorder[1],axisorder[2],axisorder[3]]
201    elif len(grid.ordinates)==2:
202        axes=['t',axisorder[1],axisorder[2]]
203    nc.addVariable(fulldata,self.id, axes)
204    nc.closeFinishedFile()
205    return pathToSubsetCSML,pathToSubsetNetCDF, totalArraySize
Note: See TracBrowser for help on using the repository browser.