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

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

fixed subsetting bug when handling times, and also removed dependence on longitude/latitude names

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