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

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

fixed bug with inline axis values

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