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

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

fixed bug where nearest neighbour search conflicted with longitude transformation

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