source: TI02-CSML/trunk/csml/API/genSubset.py @ 3049

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

2 minor bugs relating to cf output

Line 
1'''#getSubset.py contains generic subsetting code that can be reused for multiple features
2Dominic Lowe, CCLRC, 31/01/2007'''
3
4import csml
5import csmlutils
6import sys
7import Numeric
8
9def checkNeighbours(domain, **kwargs):
10    #for any non-range requests, get nearest neighbour
11    #e.g if 'latitude': (11) requested, latitude=12 may be the nearest value
12    for key in kwargs:
13        #handle single values
14        if type(kwargs[key]) is not tuple:
15            if type(domain[key]) is not list:
16                axeslist=domain[key].tolist()
17            else:
18                axeslist = domain[key]
19            nearestNeighbour=csml.API.csmlutils.nudgeSingleValuesToAxisValues(kwargs[key],axeslist)
20            if nearestNeighbour is not None:
21                kwargs[key]=nearestNeighbour   
22        else:                                 
23            tmpkey=[]
24            for val in kwargs[key]:
25                nearestNeighbour=csml.API.csmlutils.nudgeSingleValuesToAxisValues(val, domain[key], minBound=kwargs[key][0], maxBound=kwargs[key][1])               
26                if nearestNeighbour is not None:
27                    tmpkey.append(nearestNeighbour)   
28            kwargs[key]=tuple(tmpkey)           
29    return kwargs
30
31def subsetDomain(timeaxis,times, domain,**kwargs):
32    '''takes the domain and returns a subset according to the keyword selection criteria
33    time = name of time dimension
34    times = string of times
35    '''       
36    #reduce requests with same max and min eg (45,45) down to single value.
37    for kw in kwargs:
38        if type (kwargs[kw]) is tuple:
39            if kwargs[kw][0] == kwargs[kw][1]:
40                kwargs[kw]=kwargs[kw][0]
41   
42   
43    strTimes=times
44    subsetDomain={}
45    totalArraySize=1
46    for key in domain.keys():
47        straxisValues=''
48        if key in kwargs:
49            if key ==timeaxis:
50                straxisValues=strTimes
51                arraySize=len(strTimes.split())
52            elif type(kwargs[key]) is not tuple:
53                #only one value not a range
54                straxisValues=str(kwargs[key]) 
55            elif kwargs[key][0] < kwargs[key][1]:   
56                for val in domain[key]:
57                      if val is not '':
58                            if float(val) >= kwargs[key][0]:
59                                if float(val) <= kwargs[key] [1]:
60                                    straxisValues=straxisValues+ str(val) + ','
61                straxisValues=straxisValues[:-1]
62            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.
63                    for val in domain[key]:
64                        if val is not '':
65                            if val >= kwargs[key][0]:
66                                straxisValues=straxisValues+ str(val) + ','
67                    for val in domain[key]:
68                        if val is not '':
69                            if val <= kwargs[key] [1]:
70                                straxisValues=straxisValues+ str(val) + ','
71                    straxisValues=straxisValues[:-1]                                             
72        else: # this dimension has not been subsetted at all - CHECK  THIS
73            for val in domain[key]:
74                if val is not '':
75                    straxisValues=straxisValues+ str(val) + ','       
76            straxisValues=straxisValues[:-1]                       
77                                           
78        subsetDomain[key]=straxisValues
79        if key != timeaxis:
80            arraySize=len(subsetDomain[key].split(','))
81        else:
82            arraySize=len(subsetDomain[key].split())
83        totalArraySize = totalArraySize * arraySize
84   
85   
86    #Finally, check the order of the subsetDomain is the same as the order of the kwargs:
87    #eg. if the kwargs = (-90,90), then the subsetDomain  should not be 90 to -90.
88    KWordering={}   
89    for kw in kwargs:
90        if type (kwargs[kw]) is tuple:
91            if kw in ['longitude', 'time']:  #TODO should get these names from the CRS...
92                continue
93            if len(kwargs[kw]) >1: 
94                if kwargs[kw][0] > kwargs[kw][-1]:
95                    KWordering[kw]= 'highlow'
96                else:
97                    KWordering[kw]= 'lowhigh'
98           
99    for kw in KWordering:
100        if eval(subsetDomain[kw].split(',')[0]) > eval(subsetDomain[kw].split(',')[-1]):
101            subOrder='highlow'
102        else:
103            subOrder='lowhigh'
104        #if order needs reversing, reverse the values in the subset domain string           
105        if subOrder != KWordering[kw]:
106            valList=[]
107            valList= subsetDomain[kw].split(',')
108            valStr=''
109            valList.reverse()
110            for val in valList:
111                valStr=valStr + ',' + val   
112            subsetDomain[kw] = valStr[1:]
113    return subsetDomain, totalArraySize
114   
115def _setCalendar(feature, timeName, ords):
116    '''sets the cdtime calendar needed for conversion from CSML to NetCDF'''
117    calset=False
118    for gridOrd in ords:
119        if gridOrd.coordAxisLabel.CONTENT==timeName:     
120            try:
121                caltype=gridOrd.coordAxisValues.timePositionList.frame.split(':',2)[1]
122                csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
123                calset=True
124            except:pass
125    if calset!=True:
126        csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar)       
127    try:
128        caltype=gridOrd.coordAxisValues.timePositionList.frame.split(':',2)[1]
129        csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
130    except:
131        csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar)     
132        caltype='standard' 
133    return caltype
134
135
136def _getTimes(timeSelection, timeName, domain):
137    #currently supporting domain subsetting only by CRS name
138    #(but should be easy to extend later)
139    if type(timeSelection) in [str,unicode]:
140        times=(timeSelection,)
141    else:
142        times=timeSelection
143    if len(times)==0:
144        #no time specified, select all.
145        times= domain[timeName]
146    tone=csml.csmllibs.csmltime.getCDtime(times[0])
147    if len(times) == 2:
148        #then this is a range of times (t1, tn), and not a list
149        try:
150            tone=csml.csmllibs.csmltime.getCDtime(times[0])
151        except:
152            tone=times[0]
153        try:
154            ttwo=csml.csmllibs.csmltime.getCDtime(times[1])
155        except:
156            ttwo=times[1]
157        times=[]         
158        for time in domain[timeName]:
159            timeok=csml.csmllibs.csmltime.compareTimes(tone,time,ttwo)
160            if timeok ==1:
161                times.append(time)
162    return times
163
164def _getTimeToFileRatio(feature,domain, timeName):
165    if hasattr(feature.value.rangeSet, 'valueArray'):
166        if hasattr(feature.value.rangeSet.valueArray.valueComponent, 'insertedExtract'):           
167            try:
168                numFiles= len( csmlutils.listify(feature.value.rangeSet.valueArray.valueComponent.insertedExtract.components)[0].fileList.fileNames.CONTENT.split())
169            except:
170                try:
171                    numFiles= len(csmlutils.listify(feature.value.rangeSet.valueArray.valueComponent.insertedExtract.components)[0].fileName.CONTENT)
172                except:
173                    numFiles= len(csmlutils.listify(feature.value.rangeSet.valueArray.valueComponent.insertedExtract.fileName.CONTENT))
174            timeToFileRatio= len(domain[timeName])/numFiles
175            if timeToFileRatio==0:
176                timeToFileRatio=None           
177            return timeToFileRatio
178
179def getCoordTransformTable(domainSubset, crs, frame):
180    cTT=csml.parser.GridCoordinatesTable()
181    ords =[]
182    for key in domainSubset.keys():
183        go=csml.parser.GridOrdinateDescription()
184        go.coordAxisLabel=csml.parser.csString(key)
185        go.gridAxesSpanned=csml.parser.csString(key)
186        go.coordAxisValues = csml.parser.SpatialOrTemporalPositionList()
187        go.coordAxisValues.id=csml.csmllibs.csmlextra.getRandomID()
188        vals=str(domainSubset[key]).replace(',',' ')
189       
190        if key==crs.axes[crs.timeAxis]:
191            go.coordAxisValues.timePositionList=csml.parser.csString(vals)
192            go.coordAxisValues.frame = frame
193            #self.domain.key placeholder
194        else:
195            go.coordAxisValues.coordinateList=csml.parser.csString(vals) #self.domain.key placeholder
196        seqRule= csml.parser.SequenceRule()
197        seqRule.CONTENT='Linear'
198        seqRule.axisOrder='+1'  #TO DO. Work this out.
199        go.sequenceRule=seqRule
200        ords.append(go) #go needs a few more properties setting
201    cTT.gridOrdinates=ords
202    return cTT
203
204
205
206def getTheData(feature, selection, times,timeName):
207    #SOME OF THIS SHOULD PROBABLY BE IN THE DATA IO LAYER
208    #add time range back into the selection
209    if len(times)==1:
210        selection[timeName]=times[0]
211    else:
212         selection[timeName]=(times[0], times[len(times)-1])
213
214    domain = feature.domain
215    value=feature.value
216    files=[]
217    strTimes=''
218    fulldata=None
219    #get the ratio of times to files
220    timeToFileRatio = _getTimeToFileRatio(feature, domain, timeName)
221           
222    #list 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...
223    filesFetched=[]
224     
225    #get list of file extracts
226    try:
227        componentlist = csmlutils.listify(value.rangeSet.valueArray.valueComponent.insertedExtract.components)
228    except:
229        componentlist = csmlutils.listify(value.rangeSet.valueArray.valueComponent.insertedExtract)
230    for time in times:
231        listPosition=domain[timeName].index(time)
232        strTimes= strTimes + ' ' + time                   
233        for comp in componentlist:
234            if timeToFileRatio == None:
235                filePos=None
236            else:               
237                filePos=int(float(listPosition)/timeToFileRatio)
238            if filePos in filesFetched:
239                continue #already got data from this file, try next time         
240            data, fillvalue, axisorder, units=comp.getData(fileposition=filePos, **selection)         
241            if filePos is None:
242                files.append(comp.fileName.CONTENT)
243            else:
244                try:
245                    files.append(comp.fileList.fileNames.CONTENT[filePos]) #TODO, get the right file name
246                except:
247                    files.append(comp.fileName.CONTENT)
248               
249    #open the file
250           
251           
252            if fulldata is None:
253                fulldata=data
254            else:
255                newfulldata=Numeric.concatenate([fulldata,data], axis=0)   
256                fulldata=newfulldata
257            filesFetched.append(filePos)
258    if hasattr(value.rangeSet, 'valueArray'):
259        units.append(value.rangeSet.valueArray.valueComponent.uom) # final unit is that of the parameter
260    else:
261        units.append ('unitsTBA')
262    return strTimes, axisorder, units, fulldata, fillvalue
263
264def genericSubset(feature, outputdir, ncname, domain, kwargs):
265    if outputdir is not None:
266        pathToSubsetNetCDF=outputdir+'/' +ncname
267    else:
268        pathToSubsetNetCDF=ncname
269   
270    #deal with longitude requests
271    #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.
272    kwargs=csmlutils.fixLongitude(domain,kwargs)
273   
274   
275    #if request doesn't match domain points find nearest neighbours
276    kwargs=csml.API.genSubset.checkNeighbours(domain, **kwargs)
277   
278    #get the name of the time axis in the coordinate reference system
279    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
280    if type(feature) == csml.parser.GridSeriesFeature:
281        crs=cat.getCRS(feature.value.gridSeriesDomain.srsName,labels=feature.value.gridSeriesDomain.axisLabels)
282        gridordinates = feature.value.gridSeriesDomain.coordTransformTable.gridOrdinates
283    elif type(feature) == csml.parser.ProfileSeriesFeature:
284        crs=cat.getCRS(feature.value.profileSeriesDomain.srsName,labels = feature.value.profileSeriesDomain.axisLabels)
285        gridordinates = feature.value.profileSeriesDomain.coordTransformTable.gridOrdinates
286    elif type(feature) == csml.parser.TrajectoryFeature:
287        crs=cat.getCRS(feature.value.trajectoryDomain.srsName,labels = feature.value.trajectoryDomain.axisLabels)
288        gridordinates = feature.value.trajectoryDomain.coordTransformTable.gridOrdinates
289    timeName=crs.axes[crs.timeAxis]
290    #Find the original time name (timeAxis) for the file.
291    for gridOrd in gridordinates:
292        if gridOrd.coordAxisLabel.CONTENT==timeName:
293            timeAxis = gridOrd.gridAxesSpanned.CONTENT  #only works for regular grids atm
294            break
295    #set the reference system for the time axis
296    caltype=_setCalendar(feature, timeName, gridordinates)
297    try:
298        timeSelection=kwargs[timeName]
299    except KeyError:
300        timeSelection=[]
301   
302    times=_getTimes(timeSelection, timeName,domain)       
303    return pathToSubsetNetCDF, kwargs, timeAxis,timeName, caltype, times
Note: See TracBrowser for help on using the repository browser.