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

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

Improved csmlContainer mechanism to allow subsetting of multiple features and delivery in a single csml document with multiple CSML files

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