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

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

changes to !subsetGridSeries

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