source: TI02-CSML/trunk/csml/API/ops_GridSeriesFeature.py @ 1937

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

synching with repository before downloading

Line 
1''' ops_GridSeriesFeature  contains operations for GridSeriesFeatures'''
2
3import csml.parser
4import csml.csmllibs.csmltime
5import csml.csmllibs.csmldocument
6import csml.API.ops_AbstractFeature
7import csml.csmllibs.netCDFWriter
8
9import sys  #remove later
10
11def testmethod(self):
12    #print 'testmethod for gridseries feature'
13    return 'testmethod - gridseries'
14
15def getAllowedSubsettings(self):
16    return ['subsetToGridSeries']  #other operations
17
18
19def getDomain(self):
20    #returns domain as a dictionary of ordinates {name: [values], ...}
21    self.domain={}
22    for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
23        name=gridOrd.coordAxisLabel.CONTENT
24        if isinstance(gridOrd.coordAxisValues, csml.parser.FileExtract):
25            #not tested with file extracts yet: (01/01/07)
26            self.domain[name]=gridOrd.coordAxisValues.getData()
27        else:
28            vals=gridOrd.coordAxisValues.coordinateList.CONTENT
29            valList=[]
30            for val in vals.split(','):  #remove commas
31                valList.append(val)
32            self.domain[name]=valList
33    return self.domain
34
35def subsetToGridSeries(self, csmlpath=None, ncpath=None,**kwargs):
36    #pathToSubsetCSML = container.csmlpath
37    if ncpath is not None:
38        pathToSubsetNetCDF=ncpath
39    else:
40        pathToSubsetNetCDF='temp.nc'
41   
42    #deal with longitude requests
43    #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.
44    dc = self.getDomain()
45    for key in dc.keys():
46        if key == 'longitude': #how do we test if it is longitude properly?
47            for val in dc[key]:
48                if val < 0:
49                    pass
50                else:
51                    if kwargs[key][0] < 0:
52                        kwargs[key]=(kwargs[key][0]+360,kwargs[key][1])
53                    if kwargs[key][1] < 0:
54                        kwargs[key]=(kwargs[key][0],kwargs[key][1]+360)
55   
56    timeAxis=self._identifyTimeAxis(dc.keys())
57    self.times=kwargs[timeAxis]
58    self.files=[]
59    strTimes=''
60    fulldata=[]
61    if len(self.times) == 2:
62        #then this is a range of times (t1, tn)
63        try:
64            tone=ops_AbstractFeature.__getCDtime(self.times[0])
65        except:
66            tone=self.times[0]
67        try:
68            ttwo=ops_AbstractFeature.__getCDtime(self.times[1])
69        except:
70            ttwo=self.times[1]
71        self.times=[]
72        if hasattr(self, 'domain'):
73            pass
74        else:
75            self.getDomain()
76        for time in self.domain[timeAxis][0].split():
77            timeok=csml.API.ops_AbstractFeature.__compareTimes(tone,time,ttwo)
78            if timeok ==1:
79                self.times.append(time)
80    print dir(self.value.rangeSet)
81    if hasattr(self.value.rangeSet, 'aggregatedArray'):
82        #handle aggregatedArray
83        if type(self.value.rangeSet.aggregatedArray.components) is list:
84            numFiles= len(self.value.rangeSet.aggregatedArray.components[0].fileName.CONTENT.split())
85        else:
86            numFiles= len(self.value.rangeSet.aggregatedArray.components.fileName.CONTENT.split())
87        timeToFileRatio=len(self.domain[timeAxis][0].split())/numFiles
88        #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...
89        filesFetched=[]
90        #get data:
91        for time in self.times:
92            listPosition=self.domain[timeAxis][0].split().index(time)
93            strTimes= strTimes + ' ' + time
94            for comp in iter([self.value.rangeSet.aggregatedArray.components]): #iter??
95                filePos=(listPosition)/timeToFileRatio
96                if filePos in filesFetched:
97                    continue #already got data from this file, try next time
98                data=comp.getData(fileposition=filePos, times=self.times, **kwargs)
99                self.files.append(comp.fileName.split()[filePos])
100                if fulldata ==[]:
101                    fulldata = data.tolist()
102                else:
103                    for item in data.tolist():
104                        fulldata.append(item)
105                filesFetched.append(filePos)
106            axisorder = data.getAxisIds()  #will need later!
107            sys.exit()
108            for comp in iter([self.value.rangeSet.aggregatedArray.components]):
109                filePos=(listPosition)/timeToFileRatio
110                if filePos in filesFetched:
111                    continue #already got data from this file, try next time
112                data=comp.getData(fileposition=filePos, times=self.times, **kwargs)
113                self.files.append(comp.fileName.split()[filePos])
114                if fulldata ==[]:
115                    fulldata = data.tolist()
116                else:
117                    for item in data.tolist():
118                        fulldata.append(item)
119                filesFetched.append(filePos)
120            axisorder = data.getAxisIds()  #will need later!
121        try:
122            caltype=self.domain.domainReference.frame.split(':',1)[0]
123            calunits=self.domain.domainReference.frame.split(':',1)[1]
124            csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
125        except:
126            csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar)
127    elif hasattr(self.value.rangeSet, 'datablock'): #not tested
128        pass
129   
130   
131   
132    ### define domain and rangeSet to use for feature in csml document####
133    domain=csml.parser.GridSeriesDomain()
134    domain.domainReference=csml.parser.TimePositionList(timePositions=strTimes)
135    grid=csml.parser.Grid()
136    #dc = self.getDomainComplement()
137    ordinates= []
138    i=0
139    valueStore=[]  # use the values again later to generate netcdf
140    arraySize=0
141    totalArraySize=1
142    for key in dc.keys():
143        arraySize=0
144        i=i+1
145        god=csml.parser.GridOrdinateDescription()
146        god.gridAxesSpanned='dim%s'%i
147        god.sequenceRule='+x+y+z'
148        god.definesAxis=key
149        straxisValues=''
150        #now deal with each argument:
151
152        if key in kwargs:
153            if kwargs[key][0] < kwargs[key][1]:   
154                for val in dc[key]:
155                    if val is not '':
156                        if float(val) >= kwargs[key][0]:
157                            if float(val) <= kwargs[key] [1]:
158                                arraySize=arraySize+1
159                                straxisValues=straxisValues+ str(val) + ', '
160            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.
161                    for val in dc[key]:
162                        if val is not '':
163                            if val >= kwargs[key][0]:
164                                arraySize=arraySize+1
165                                straxisValues=straxisValues+ str(val) + ', '
166                    for val in dc[key]:
167                        if val is not '':
168                            if val <= kwargs[key] [1]:
169                                arraySize=arraySize+1
170                                straxisValues=straxisValues+ str(val) + ', '
171        else: # this dimension has not been subsetted at all
172            for val in dc[key]:
173                if val is not '':
174                    arraySize=arraySize+1
175                    straxisValues=straxisValues+ str(val) + ', '       
176        totalArraySize=totalArraySize*arraySize
177        god.axisValues=straxisValues[:-2]
178        ordinates.append(god)
179    totalArraySize=totalArraySize*len(self.times)
180    grid.ordinates=ordinates
181    domain.domainComplement=grid
182    rangeSet=csml.parser.RangeSet()
183    rangeSet.arrayDescriptor=csml.parser.NetCDFExtract(id=self.id,fileName=pathToSubsetNetCDF,variableName=self.id,arraySize=[arraySize])
184
185    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
186    subsettedFeature=csmlWrap.createGridSeriesFeature(domain,rangeSet,datasetID="A",featureID="B",description="C")
187    #container.appendFeature(subsettedFeature)
188
189    ### write netcdf using NCWriter class (wraps cdms) ###
190    nc=csml.csmllibs.netCDFWriter.NCwriter(pathToSubsetNetCDF)
191    floatTimes=[]
192    for time in self.times:
193        time=csml.API.ops_AbstractFeature.__getCDtime(time).torel(calunits)
194        floatTimes.append(time.value)
195    nc.addAxis('t',floatTimes,isTime=1,units=calunits,calendar=caltype)
196    for ordinate in ordinates:
197        lon,lat=None,None
198        if ordinate.definesAxis=='longitude':
199            lon=1
200        if ordinate.definesAxis=='latitude':
201            lat=1
202        #convert to list
203        vals=[]
204        for val in ordinate.axisValues.split(','):
205            vals.append(float(val))
206        nc.addAxis(ordinate.definesAxis,vals,isLon=lon,isLat=lat,units='')#to do, units attribute for CF compliance
207    if len(ordinates)==3:
208        axes=['t',axisorder[1],axisorder[2],axisorder[3]]
209    elif len(grid.ordinates)==2:
210        axes=['t',axisorder[1],axisorder[2]]
211    nc.addVariable(fulldata,self.id, axes,units='') #to do, units attribute for CF compliance
212    nc.closeFinishedFile()
213    return subsettedFeature, pathToSubsetNetCDF
214    #container.attachNetCDFFile(nc)
215    #MUST be supplied with a CSMLContainer object to store the subsetted feature in
216    return subsettedFeature, pathToSubsetNetCDF
Note: See TracBrowser for help on using the repository browser.