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

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

making selection more robust

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
6
7
8def subsetDomain(timeaxis,times, domain,**kwargs):
9    '''takes the domain and returns a subset according to the keyword selection criteria
10    time = name of time dimension
11    times = string of times
12    '''   
13    strTimes=times
14    subsetDomain={}
15    totalArraySize=1
16    for key in domain.keys():
17        straxisValues=''
18        if key in kwargs:
19            if key ==timeaxis:
20                straxisValues=strTimes
21                arraySize=len(strTimes.split())
22            elif type(kwargs[key]) is not tuple:
23                #only one value not a range.
24                straxisValues=str(kwargs[key]) 
25            elif kwargs[key][0] < kwargs[key][1]:   
26                for val in domain[key]:
27                      if val is not '':
28                            if float(val) >= kwargs[key][0]:
29                                if float(val) <= kwargs[key] [1]:
30                                    straxisValues=straxisValues+ str(val) + ','
31                straxisValues=straxisValues[:-1]
32            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.
33                    for val in domain[key]:
34                        if val is not '':
35                            if val >= kwargs[key][0]:
36                                straxisValues=straxisValues+ str(val) + ','
37                    for val in domain[key]:
38                        if val is not '':
39                            if val <= kwargs[key] [1]:
40                                straxisValues=straxisValues+ str(val) + ','
41                    straxisValues=straxisValues[:-1]                                             
42        else: # this dimension has not been subsetted at all - CHECK  THIS
43            for val in domain[key]:
44                if val is not '':
45                    straxisValues=straxisValues+ str(val) + ','       
46            straxisValues=straxisValues[:-1]                       
47                                           
48        subsetDomain[key]=straxisValues
49        if key != timeaxis:
50            arraySize=len(subsetDomain[key].split(','))
51        else:
52            arraySize=len(subsetDomain[key].split())
53        totalArraySize = totalArraySize * arraySize
54    return subsetDomain, totalArraySize
55   
56def _setCalendar(feature, timeName):
57    '''sets the cdtime calendar needed for conversion from CSML to NetCDF'''
58    calset=False
59    for gridOrd in feature.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
60        if gridOrd.coordAxisLabel.CONTENT==timeName:
61            try:
62                caltype=gridOrd.coordAxisValues.frame.split(':',1)[0]
63                calunits=gridOrd.coordAxisValues.frame.split(':',1)[1]
64                csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
65                calset=True
66            except:pass
67    if calset!=True:
68        csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar)   
69    try:
70        caltype=gridOrd.coordAxisValues.frame.split(':',1)[0]
71        calunits=gridOrd.coordAxisValues.frame.split(':',1)[1]
72        csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
73    except:
74        csml.csmllibs.csmltime.setcdtimeCalendar(csml.csmllibs.csmltime.cdtime.DefaultCalendar)
75    return calunits, caltype
76
77
78def _getTimes(timeSelection, timeName, domain):
79    #currently supporting domain subsetting only by CRS name
80    #(but should be easy to extend later)
81    if type(timeSelection) is str:
82        times=(timeSelection,)
83    else:
84        times=timeSelection
85   
86   
87    #times=timeSelection
88   
89    print times
90    print len(times)
91    if len(times)==0:
92        #no time specified, select all.
93        times= domain[timeName]
94    tone=csml.csmllibs.csmltime.getCDtime(times[0])
95    if len(times) == 2:
96        #then this is a range of times (t1, tn), and not a list
97        try:
98            tone=csml.csmllibs.csmltime.getCDtime(times[0])
99        except:
100            tone=times[0]
101        try:
102            ttwo=csml.csmllibs.csmltime.getCDtime(times[1])
103        except:
104            ttwo=times[1]
105        times=[]         
106        for time in domain[timeName]:
107            timeok=csml.csmllibs.csmltime.compareTimes(tone,time,ttwo)
108            if timeok ==1:
109                times.append(time)
110    return times
111
112def _getTimeToFileRatio(feature,domain, timeName):
113    if hasattr(feature.value.rangeSet, 'valueArray'):
114        if hasattr(feature.value.rangeSet.valueArray.valueComponent.quantityList, '__insertedExtract'):
115            numFiles= len( csmlutils.listify(feature.value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components)[0].fileList.fileNames.CONTENT.split())
116            return len(domain[timeName])/numFiles
117
118def getCoordTransformTable(domainSubset, crs):
119    cTT=csml.parser.GridCoordinatesTable()
120    ords =[]
121    for key in domainSubset.keys():
122        go=csml.parser.GridOrdinateDescription()
123        go.coordAxisLabel=csml.parser.csString(key)
124        go.gridAxesSpanned=csml.parser.csString(key)
125        go.coordAxisValues = csml.parser.SpatialOrTemporalPositionList()
126        if key==crs.axes[crs.timeAxis]:
127            go.coordAxisValues.timePositionList=csml.parser.csString(domainSubset[key]) #self.domain.key placeholder
128        else:
129            go.coordAxisValues.coordinateList=csml.parser.csString(domainSubset[key]) #self.domain.key placeholder
130        ords.append(go) #go needs a few more properties setting
131    cTT.gridOrdinates=ords
132    return cTT
133
134
135
136def getTheData(feature, selection, times,timeName):
137    #SOME OF THIS SHOULD PROBABLY BE IN THE DATA IO LAYER
138    domain = feature.domain
139    value=feature.value
140    files=[]
141    strTimes=''
142    fulldata=[]
143    #get the ratio of times to files
144    timeToFileRatio = _getTimeToFileRatio(feature, domain, timeName)
145       
146           
147    #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...
148    filesFetched=[]
149    for time in times:
150        listPosition=domain[timeName].index(time)
151        strTimes= strTimes + ' ' + time
152        for comp in csmlutils.listify(value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components): 
153            filePos=int(float(listPosition)/timeToFileRatio)
154            if filePos in filesFetched:
155                continue #already got data from this file, try next time
156            data, fillvalue, axisorder, units=comp.getData(fileposition=filePos, **selection)
157            files.append(comp.fileList.fileNames.CONTENT[filePos]) #TODO, get the right file name
158            if fulldata ==[]:
159                fulldata = data
160            else:
161                for item in data:
162                    fulldata.append(item)
163            filesFetched.append(filePos)
164    units.append(value.rangeSet.valueArray.valueComponent.quantityList.uom) # final unit is that of the
165    return strTimes, axisorder, units, fulldata, fillvalue
166
167def genericSubset(feature, csmlpath, ncpath, domain, kwargs):
168    if ncpath is not None:
169        pathToSubsetNetCDF=ncpath
170    else:
171        pathToSubsetNetCDF='temp.nc'
172   
173    #deal with longitude requests
174    #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.
175    kwargs=csmlutils.fixLongitude(domain,kwargs)
176       
177    #get the name of the time axis in the coordinate reference system
178    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
179    crs=cat.getCRS(feature.value.gridSeriesDomain.srsName)
180    timeName=crs.axes[crs.timeAxis]
181     #set the reference system for the time axis
182    calunits, caltype=_setCalendar(feature, timeName)
183   
184    #selection by time - get explcit time values, not just a range mint, maxt
185    try:
186        timeSelection=kwargs[timeName]
187    except KeyError:
188        timeSelection=[]
189    times=_getTimes(timeSelection, timeName,domain)
190   
191   
192    return pathToSubsetNetCDF, kwargs, timeName, calunits, caltype, times
Note: See TracBrowser for help on using the repository browser.