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

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

added fix to check ordering of output axes, fixes latitude reversal problem

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