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

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

added !getData method to GridSeriesFeature (works for inline data only at the moment)

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