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

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

subsetting seems to be working, but needs more testing and tidying of code

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            for val in vals.split(): 
30                valList.append(val)
31            self.domain[name]=valList
32    return self.domain
33
34def __subsetDomain(self,time=None,times=None,**kwargs):
35    #takes the domain and returns a subset according to the keyword selection criteria
36    strTimes=times
37    subsetDomain={}
38    totalArraySize=1
39    for key in self.domain.keys():
40        straxisValues=''
41        if key in kwargs:
42            if key ==time:
43                straxisValues=strTimes
44                arraySize=len(strTimes.split())
45            elif kwargs[key][0] < kwargs[key][1]:   
46                for val in self.domain[key]:
47                      if val is not '':
48                            if float(val) >= kwargs[key][0]:
49                                if float(val) <= kwargs[key] [1]:
50                                    straxisValues=straxisValues+ str(val) + ','
51                straxisValues=straxisValues[:-1]
52            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.
53                    for val in self.domain[key]:
54                        if val is not '':
55                            if val >= kwargs[key][0]:
56                                straxisValues=straxisValues+ str(val) + ','
57                    for val in self.domain[key]:
58                        if val is not '':
59                            if val <= kwargs[key] [1]:
60                                straxisValues=straxisValues+ str(val) + ','
61                    straxisValues=straxisValues[:-1]                                             
62        else: # this dimension has not been subsetted at all - CHECK  THIS
63            for val in self.domain[key]:
64                if val is not '':
65                    straxisValues=straxisValues+ str(val) + ','       
66            straxisValues=straxisValues[:-1]                       
67                                   
68        subsetDomain[key]=straxisValues
69        if key != time:
70            arraySize=len(subsetDomain[key].split(','))
71            #print subsetDomain[key].split(',')
72            #print key + ' ' + str(arraySize)
73        else:
74            arraySize=len(subsetDomain[key].split())
75            #print key + ' ' + str(arraySize)
76        totalArraySize = totalArraySize * arraySize
77    #print subsetDomain
78    #print totalArraySize
79    #sys.exit()
80    return subsetDomain, totalArraySize
81
82
83def subsetToGridSeries(self, csmlpath=None, ncpath=None,**kwargs):
84    #pathToSubsetCSML = container.csmlpath
85    if ncpath is not None:
86        pathToSubsetNetCDF=ncpath
87    else:
88        pathToSubsetNetCDF='temp.nc'
89       
90    self.getDomain()
91    #Now get the data subset:
92   
93    #TODO - incoporate the crsCatalogue. into this
94    #deal with longitude requests
95    #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.
96    kwargs=csmlutils.fixLongitude(self.domain,kwargs)
97    #deal with times   
98    #need to find the time axis:
99    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
100    crs=cat.getCRS(self.value.gridSeriesDomain.srsName)
101    #get the position of the time axis in crs/ axis labels
102    timeAxisPos=crs.timeAxis
103    #get the name of the time axis in the crs
104    timeName=crs.axes[timeAxisPos]
105    #get the ordinate with that name and find the original time name for the file.
106    for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
107        if gridOrd.coordAxisLabel.CONTENT==timeName:
108            timeAxis = gridOrd.gridAxesSpanned.CONTENT  #only works for regular grids atm
109        else: 
110            print 'crs not supported'
111   
112   
113    #here, timeName is the name given to the time axis by the CRS,
114    #and timeAxis is the name given to the time axis in the file.
115    #currently supporting domain subsetting only by CRS name
116    #(but should be easy to extend later)
117    self.times=kwargs[timeName]
118    self.files=[]
119    strTimes=''
120    fulldata=[]
121   
122    #set the reference system for the time axis
123    calset=False
124    for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
125        if gridOrd.coordAxisLabel.CONTENT==timeName:
126            try:
127                caltype=gridOrd.coordAxisValues.frame.split(':',1)[0]
128                calunits=gridOrd.coordAxisValues.frame.split(':',1)[1]
129                csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
130                calset=True
131            except:pass
132    if calset!=True:
133        csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar)   
134    try:
135        caltype=self.domain.domainReference.frame.split(':',1)[0]
136        calunits=self.domain.domainReference.frame.split(':',1)[1]
137        csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
138    except:
139        csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar)
140 
141       
142    if len(self.times) == 2:
143        #then this is a range of times (t1, tn)
144        try:
145            tone=csml.API.ops_AbstractFeature.__getCDtime(self.times[0])
146        except:
147            tone=self.times[0]
148        try:
149            ttwo=csml.API.ops_AbstractFeature.__getCDtime(self.times[1])
150        except:
151            ttwo=self.times[1]
152        self.times=[]
153           
154        for time in self.domain[timeName]:
155            timeok=csml.API.ops_AbstractFeature.__compareTimes(tone,time,ttwo)
156            if timeok ==1:
157                self.times.append(time)
158   
159    if hasattr(self.value.rangeSet, 'valueArray'):
160        if hasattr(self.value.rangeSet.valueArray.valueComponent.quantityList, '__insertedExtract'):
161            numFiles= len( csmlutils.listify(self.value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components)[0].fileList.fileNames.CONTENT.split())
162            timeToFileRatio=len(self.domain[timeName])/numFiles
163           
164            #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...
165            filesFetched=[]
166            #get data:
167            selection={}
168           
169            #Get names of variables in file
170            cat=csml.csmllibs.csmlcrs.CRSCatalogue()
171            crs=cat.getCRS(self.value.gridSeriesDomain.srsName) 
172           
173            for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
174                selection[gridOrd.gridAxesSpanned.CONTENT]=kwargs[gridOrd.coordAxisLabel.CONTENT]
175           
176                         
177            for time in self.times:
178                listPosition=self.domain[timeName].index(time)
179                strTimes= strTimes + ' ' + time
180                for comp in csmlutils.listify(self.value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components): 
181                    filePos=int(float(listPosition)/timeToFileRatio)
182                    if filePos in filesFetched:
183                        continue #already got data from this file, try next time
184                    print selection
185                   
186                    data=comp.getData(fileposition=filePos, **selection)
187                    print data
188                   
189                     
190                    self.files.append(comp.fileList.fileNames.CONTENT[filePos]) #TODO, get the right file name
191                    if fulldata ==[]:
192                        fulldata = data.tolist()
193                        print fulldata
194                    else:
195                        for item in data.tolist():
196                            fulldata.append(item)
197                    filesFetched.append(filePos)
198                axisorder = data.getAxisIds()  #will need later!
199
200    elif hasattr(self.value.rangeSet, 'datablock'): #not tested
201        pass
202   
203    #Okay, got the data now. Need to write out CSML and NetCDF files...
204 
205
206    #Writing out the CSML file
207    # define domain/coverage  to use in 'value' attribute   
208   
209    domain=csml.parser.GridSeriesDomain()
210    grid=csml.parser.GridCoordinatesTable()
211    domainSubset, totalArraySize=__subsetDomain(self, time=timeName,times=strTimes,**kwargs)
212    cTT=csml.parser.GridCoordinatesTable()
213    ords =[]
214    for key in domainSubset.keys():
215        go=csml.parser.GridOrdinateDescription()
216        go.coordAxisLabel=csml.parser.csString(key)
217        go.gridAxesSpanned=csml.parser.csString(key)
218        go.coordAxisValues = csml.parser.SpatialOrTemporalPositionList()
219        if key==crs.axes[crs.timeAxis]:
220            go.coordAxisValues.timePositionList=csml.parser.csString(domainSubset[key]) #self.domain.key placeholder
221        else:
222            go.coordAxisValues.coordinateList=csml.parser.csString(domainSubset[key]) #self.domain.key placeholder
223        ords.append(go) #go needs a few more properties setting
224    cTT.gridOrdinates=ords
225    grid=self.value.gridSeriesDomain.coordTransformTable
226    domain.coordTransformTable=cTT
227    rangeSet=csml.parser.RangeSet()
228    rangeSet.arrayDescriptor=csml.parser.NetCDFExtract(id=self.id,fileName=csml.parser.csString(pathToSubsetNetCDF),variableName=csml.parser.csString(self.id),arraySize=csml.parser.csString(totalArraySize))
229    cvg=csml.parser.GridSeriesCoverage()
230    cvg.rangeSet=rangeSet
231    cvg.gridSeriesDomain=domain   
232    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
233    subsettedFeature=csmlWrap.createGridSeriesFeature(value=cvg,datasetID="A",featureID="B",description="C")
234   
235   
236    ### write netcdf using NCWriter class (wraps cdms) ###
237    nc=csml.csmllibs.netCDFWriter.NCwriter(pathToSubsetNetCDF)
238    floatTimes=[]
239    for time in self.times:
240        time=csml.API.ops_AbstractFeature.__getCDtime(time).torel(calunits)
241        floatTimes.append(time.value)
242    nc.addAxis('t',floatTimes,isTime=1,units=calunits,calendar=caltype)
243    for ord in ords:
244        print ord.coordAxisLabel.CONTENT
245        vals=[]
246        lon,lat=None,None
247        if ord.coordAxisLabel.CONTENT=='Time':
248            continue
249        else:
250            for val in ord.coordAxisValues.coordinateList.CONTENT.split(','):
251                if val != ' ':
252                    vals.append(float(val)) 
253                    print vals
254        if ord.coordAxisLabel.CONTENT=='Lon':
255            lon=1
256            name='longitude'
257        elif ord.coordAxisLabel.CONTENT=='Lat':
258            lat=1
259            name='latitude'
260        else:
261            name=ord.coordAxisLabel.CONTENT
262 
263        nc.addAxis(name,vals,isLon=lon,isLat=lat,units='')#to do, units attribute for CF compliance
264    if len(ords)==3:
265        axes=['t',axisorder[1],axisorder[2]]
266    elif len(ords)==2:
267        axes=['t',axisorder[1]]
268    #print fulldata
269    print axes
270    nc.addVariable(fulldata,self.id, axes,units='') #to do, units attribute for CF compliance
271    print 'added'
272    nc.closeFinishedFile()
273    return subsettedFeature, pathToSubsetNetCDF
274   
Note: See TracBrowser for help on using the repository browser.