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

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

added EnvelopeWithTimePeriod?, also, more on subsetting

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
8import csmlutils
9
10import sys  #remove later
11
12def testmethod(self):
13    #print 'testmethod for gridseries feature'
14    return 'testmethod - gridseries'
15
16def getAllowedSubsettings(self):
17    return ['subsetToGridSeries']  #other operations
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 hasattr(gridOrd.coordAxisValues, '__insertedExtract'):
25            self.domain[name]=gridOrd.coordAxisValues.__insertedExtract.getData()
26        else:
27            vals=gridOrd.coordAxisValues.coordinateList.CONTENT
28            valList=[]
29            print 'splitting'
30            for val in vals.split(): 
31                valList.append(val)
32                print val
33            self.domain[name]=valList
34    return self.domain
35
36def subsetToGridSeries(self, csmlpath=None, ncpath=None,**kwargs):
37    #pathToSubsetCSML = container.csmlpath
38    if ncpath is not None:
39        pathToSubsetNetCDF=ncpath
40    else:
41        pathToSubsetNetCDF='temp.nc'
42   
43    self.getDomain()
44   
45    #TODO - incoporate the crsCatalogue. into this
46    #deal with longitude requests
47    #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.
48    kwargs=csmlutils.fixLongitude(self.domain,kwargs)
49   
50    #deal with times   
51    #need to find the time axis:
52    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
53    crs=cat.getCRS(self.value.gridSeriesDomain.srsName)
54    #get the position of the time axis in crs/ axis labels
55    timeAxisPos=crs.timeAxis
56    #get the name of the time axis in the crs
57    timeName=crs.axes[timeAxisPos]
58    #get the ordinate with that name and find the original time name for the file.
59    for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
60        if gridOrd.coordAxisLabel.CONTENT==timeName:
61            timeAxis = gridOrd.gridAxesSpanned.CONTENT  #only works for regular grids atm
62        else: 
63            print 'crs not supported'
64   
65   
66    #here, timeName is the name given to the time axis by the CRS,
67    #and timeAxis is the name given to the time axis in the file.
68    #currently supporting domain subsetting only by CRS name
69    #(but should be easy to extend later)
70    self.times=kwargs[timeName]
71    self.files=[]
72    strTimes=''
73    fulldata=[]
74   
75    #set the reference system for the time axis
76    calset=False
77    for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
78        if gridOrd.coordAxisLabel.CONTENT==timeName:
79            try:
80                caltype=gridOrd.coordAxisValues.frame.split(':',1)[0]
81                calunits=gridOrd.coordAxisValues.frame.split(':',1)[1]
82                csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
83                calset=True
84            except:pass
85    if calset!=True:
86        csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar)   
87    try:
88        caltype=self.domain.domainReference.frame.split(':',1)[0]
89        calunits=self.domain.domainReference.frame.split(':',1)[1]
90        csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
91    except:
92        csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar)
93 
94       
95    if len(self.times) == 2:
96        #then this is a range of times (t1, tn)
97        try:
98            tone=csml.API.ops_AbstractFeature.__getCDtime(self.times[0])
99        except:
100            tone=self.times[0]
101        try:
102            ttwo=csml.API.ops_AbstractFeature.__getCDtime(self.times[1])
103        except:
104            ttwo=self.times[1]
105        self.times=[]
106           
107        for time in self.domain[timeName]:
108            timeok=csml.API.ops_AbstractFeature.__compareTimes(tone,time,ttwo)
109            if timeok ==1:
110                self.times.append(time)
111     #okay up to here!
112   
113    if hasattr(self.value.rangeSet, 'valueArray'):
114        if hasattr(self.value.rangeSet.valueArray.valueComponent.quantityList, '__insertedExtract'):
115            numFiles= len( csmlutils.listify(self.value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components)[0].fileList.fileNames.CONTENT.split())
116           
117            timeToFileRatio=len(self.domain[timeName])/numFiles
118            #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...
119            filesFetched=[]
120            #get data:
121            selection={}
122            for kw in kwargs:
123                        if kw != timeAxis:
124                            selection[kw]=kwargs[kw]
125            for time in self.times:
126                print time
127                listPosition=self.domain[timeName].index(time)
128                strTimes= strTimes + ' ' + time
129                for comp in csmlutils.listify(self.value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components): 
130                    filePos=(listPosition)/timeToFileRatio
131                    if filePos in filesFetched:
132                        continue #already got data from this file, try next time
133                   
134                    data=comp.getData(fileposition=filePos, **selection)
135                    print comp.fileList.fileNames.CONTENT
136                    self.files.append(comp.fileList.fileNames.CONTENT[filePos]) #TODO, get the right file name
137                    if fulldata ==[]:
138                        fulldata = data.tolist()
139                    else:
140                        for item in data.tolist():
141                            fulldata.append(item)
142                    filesFetched.append(filePos)
143                axisorder = data.getAxisIds()  #will need later!
144                for comp in csmlutils.listify(self.value.rangeSet.aggregatedArray.components):
145                    filePos=(listPosition)/timeToFileRatio
146                    if filePos in filesFetched:
147                        continue #already got data from this file, try next time
148                    data=comp.getData(fileposition=filePos, **selection)
149                    self.files.append(comp.fileName.CONTENT.split()[filePos])
150                    if fulldata ==[]:
151                        fulldata = data.tolist()
152                    else:
153                        for item in data.tolist():
154                            fulldata.append(item)
155                    filesFetched.append(filePos)
156                axisorder = data.getAxisIds()  #will need later!
157
158    elif hasattr(self.value.rangeSet, 'datablock'): #not tested
159        pass
160
161 
162
163    ### define domain and rangeSet to use for feature in csml document####
164    #TODO - rewrite this section for version 2 parser ###
165   
166    #domain=csml.parser.GridSeriesDomain()
167    #domain.domainReference=csml.parser.TimePositionList(timePositions=strTimes)
168    #grid=csml.parser.Grid()
169    ##dc = self.getDomainComplement()
170    ordinates= []
171    i=0
172    valueStore=[]  # use the values again later to generate netcdf
173    arraySize=0
174    totalArraySize=1
175    for key in self.domain.keys():
176        arraySize=0
177        i=i+1
178        god=csml.parser.GridOrdinateDescription()
179        god.gridAxesSpanned='dim%s'%i
180        god.sequenceRule='+x+y+z'
181        god.definesAxis=key
182        straxisValues=''
183        #now deal with each argument:
184   
185        if key in kwargs:
186            if key ==timeName:
187                continue #dealt with time earlier           
188            elif kwargs[key][0] < kwargs[key][1]:   
189                for val in self.domain[key]:
190                    if val is not '':
191                            if float(val) >= kwargs[key][0]:
192                                if float(val) <= kwargs[key] [1]:
193                                    arraySize=arraySize+1
194                                    straxisValues=straxisValues+ str(val) + ', '
195            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.
196                    for val in self.domain[key]:
197                        if val is not '':
198                            if val >= kwargs[key][0]:
199                                arraySize=arraySize+1
200                                straxisValues=straxisValues+ str(val) + ', '
201                    for val in self.domain[key]:
202                        if val is not '':
203                            if val <= kwargs[key] [1]:
204                                arraySize=arraySize+1
205                                straxisValues=straxisValues+ str(val) + ', '
206        else: # this dimension has not been subsetted at all
207            for val in self.domain[key]:
208                if val is not '':
209                    arraySize=arraySize+1
210                    straxisValues=straxisValues+ str(val) + ', '       
211        totalArraySize=totalArraySize*arraySize
212        god.axisValues=straxisValues[:-2]
213        ordinates.append(god)
214    #totalArraySize=totalArraySize*len(self.times)
215    #grid.ordinates=ordinates
216    #domain.domainComplement=grid
217    #rangeSet=csml.parser.RangeSet()
218    #rangeSet.arrayDescriptor=csml.parser.NetCDFExtract(id=self.id,fileName=pathToSubsetNetCDF,variableName=self.id,arraySize=[arraySize])
219
220    #csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
221    #subsettedFeature=csmlWrap.createGridSeriesFeature(domain,rangeSet,datasetID="A",featureID="B",description="C")
222    ##container.appendFeature(subsettedFeature)
223
224    ### write netcdf using NCWriter class (wraps cdms) ###
225    nc=csml.csmllibs.netCDFWriter.NCwriter(pathToSubsetNetCDF)
226    floatTimes=[]
227    for time in self.times:
228        time=csml.API.ops_AbstractFeature.__getCDtime(time).torel(calunits)
229        floatTimes.append(time.value)
230    nc.addAxis('t',floatTimes,isTime=1,units=calunits,calendar=caltype)
231    #USE CRS!   
232    for ordinate in ordinates:
233        lon,lat=None,None
234        if ordinate.definesAxis=='longitude':
235            lon=1
236        if ordinate.definesAxis=='latitude':
237            lat=1
238        #convert to list
239        vals=[]
240        for val in ordinate.axisValues.split(','):
241            vals.append(float(val))
242        nc.addAxis(ordinate.definesAxis,vals,isLon=lon,isLat=lat,units='')#to do, units attribute for CF compliance
243    if len(ordinates)==3:
244        axes=['t',axisorder[1],axisorder[2],axisorder[3]]
245    elif len(ordinates)==2:
246        axes=['t',axisorder[1],axisorder[2]]
247    nc.addVariable(fulldata,self.id, axes,units='') #to do, units attribute for CF compliance
248    nc.closeFinishedFile()
249    return subsettedFeature, pathToSubsetNetCDF
250    #container.attachNetCDFFile(nc)
251    #MUST be supplied with a CSMLContainer object to store the subsetted feature in
252    return subsettedFeature, pathToSubsetNetCDF
Note: See TracBrowser for help on using the repository browser.