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

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

subsetting working, but not complete

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