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

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

subsetting ProfileSeries? to Profile working more or less
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, gridnames, **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[gridnames[key]]) is not list:
16                axeslist=domain[gridnames[key]].tolist()
17            else:
18                axeslist = domain[gridnames[key]]
19            nearestNeighbour=csml.API.csmlutils.nudgeSingleValuesToAxisValues(kwargs[key],axeslist)
20            if nearestNeighbour is not None:
21                kwargs[key]=nearestNeighbour   
22        else:
23            #handle ranges
24            tmpkey=[]
25            for val in kwargs[key]:
26                nearestNeighbour=csml.API.csmlutils.nudgeSingleValuesToAxisValues(val, domain[gridnames[key]])
27                if nearestNeighbour is not None:
28                    tmpkey.append(nearestNeighbour)   
29                kwargs[key]=tuple(tmpkey)
30    return kwargs
31
32def subsetDomain(timeaxis,times, domain,**kwargs):
33    '''takes the domain and returns a subset according to the keyword selection criteria
34    time = name of time dimension
35    times = string of times
36    '''   
37    strTimes=times
38    subsetDomain={}
39    totalArraySize=1
40    for key in domain.keys():
41        print key
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    return subsetDomain, totalArraySize
80   
81def _setCalendar(feature, timeName, ords):
82    '''sets the cdtime calendar needed for conversion from CSML to NetCDF'''
83    calset=False
84    for gridOrd in ords:
85        if gridOrd.coordAxisLabel.CONTENT==timeName:
86            try:
87                caltype=gridOrd.coordAxisValues.frame.split(':',1)[0]
88                calunits=gridOrd.coordAxisValues.frame.split(':',1)[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.frame.split(':',1)[0]
96        calunits=gridOrd.coordAxisValues.frame.split(':',1)[1]
97        csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
98    except:
99        csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar)
100    return calunits, 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.quantityList, '__insertedExtract'):
135            numFiles= len( csmlutils.listify(feature.value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components)[0].fileList.fileNames.CONTENT.split())
136            return len(domain[timeName])/numFiles
137
138def getCoordTransformTable(domainSubset, crs):
139    cTT=csml.parser.GridCoordinatesTable()
140    ords =[]
141    for key in domainSubset.keys():
142        go=csml.parser.GridOrdinateDescription()
143        go.coordAxisLabel=csml.parser.csString(key)
144        go.gridAxesSpanned=csml.parser.csString(key)
145        go.coordAxisValues = csml.parser.SpatialOrTemporalPositionList()
146        if key==crs.axes[crs.timeAxis]:
147            go.coordAxisValues.timePositionList=csml.parser.csString(domainSubset[key]) #self.domain.key placeholder
148        else:
149            go.coordAxisValues.coordinateList=csml.parser.csString(domainSubset[key]) #self.domain.key placeholder
150        ords.append(go) #go needs a few more properties setting
151    cTT.gridOrdinates=ords
152    return cTT
153
154
155
156def getTheData(feature, selection, times,timeName):
157    #SOME OF THIS SHOULD PROBABLY BE IN THE DATA IO LAYER
158    domain = feature.domain
159    value=feature.value
160    files=[]
161    strTimes=''
162    fulldata=None
163    #get the ratio of times to files
164    timeToFileRatio = _getTimeToFileRatio(feature, domain, timeName)
165       
166           
167    #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...
168    filesFetched=[]
169    for time in times:
170        listPosition=domain[timeName].index(time)
171        strTimes= strTimes + ' ' + time
172        for comp in csmlutils.listify(value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components): 
173            filePos=int(float(listPosition)/timeToFileRatio)
174            if filePos in filesFetched:
175                continue #already got data from this file, try next time
176            data, fillvalue, axisorder, units=comp.getData(fileposition=filePos, **selection)
177            files.append(comp.fileList.fileNames.CONTENT[filePos]) #TODO, get the right file name
178            if fulldata is None:
179                fulldata=data
180            else:
181                newfulldata=Numeric.concatenate([fulldata,data], axis=0)   
182                fulldata=newfulldata
183            filesFetched.append(filePos)
184    units.append(value.rangeSet.valueArray.valueComponent.quantityList.uom) # final unit is that of the
185    return strTimes, axisorder, units, fulldata, fillvalue
186
187def genericSubset(feature, csmlpath, ncpath, domain, kwargs):
188    if ncpath is not None:
189        pathToSubsetNetCDF=ncpath
190    else:
191        pathToSubsetNetCDF='temp.nc'
192   
193    #deal with longitude requests
194    #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.
195    kwargs=csmlutils.fixLongitude(domain,kwargs)
196       
197    #get the name of the time axis in the coordinate reference system
198    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
199    if type(feature) == csml.parser.GridSeriesFeature:
200        crs=cat.getCRS(feature.value.gridSeriesDomain.srsName)
201        gridordinates = feature.value.gridSeriesDomain.coordTransformTable.gridOrdinates
202    elif type(feature) == csml.parser.ProfileSeriesFeature:
203        crs=cat.getCRS(feature.value.profileSeriesDomain.srsName)
204        gridordinates = feature.value.profileSeriesDomain.coordTransformTable.gridOrdinates
205    timeName=crs.axes[crs.timeAxis]
206    #Find the original time name (timeAxis) for the file.
207    for gridOrd in gridordinates:
208        if gridOrd.coordAxisLabel.CONTENT==timeName:
209            timeAxis = gridOrd.gridAxesSpanned.CONTENT  #only works for regular grids atm
210            break
211    #set the reference system for the time axis
212    calunits, caltype=_setCalendar(feature, timeName, gridordinates)
213   
214    #selection by time - get explcit time values, not just a range mint, maxt
215   
216    try:
217        timeSelection=kwargs[timeAxis]
218    except KeyError:
219        timeSelection=[]
220       
221    times=_getTimes(timeSelection, timeName,domain)   
222   
223    return pathToSubsetNetCDF, kwargs, timeAxis,timeName, calunits, caltype, times
Note: See TracBrowser for help on using the repository browser.