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

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

fixed problems with longitude subsetting

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