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

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

unravelled grid/axis name confusion

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