source: TI02-CSML/trunk/csml/API/ops_TrajectoryFeature.py @ 2767

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

subset trajectory by times working, added some standard names support

Line 
1''' ops_TrajectoryFeature  contains operations for TrajectoryFeature'''
2
3import csml, cdms, MV
4import sys
5import copy
6
7def testmethod(self):
8    print 'testmethod for TrajectoryFeature'
9    return 'testmethod - TrajectoryFeature'
10
11def getAllowedSubsettings(self):
12    return ['subsetToTrajectory'] 
13
14
15def getDomain(self):
16    #returns domain as a dictionary of ordinates {name: [values], ...}
17    self.domain={}
18    self.gridnames={}
19    for gridOrd in self.value.trajectoryDomain.coordTransformTable.gridOrdinates:
20        name=gridOrd.coordAxisLabel.CONTENT
21        self.gridnames[gridOrd.gridAxesSpanned.CONTENT]=name
22        if hasattr(gridOrd.coordAxisValues, 'insertedExtract'):
23            self.domain[name], fill, axisorder, units=gridOrd.coordAxisValues.insertedExtract.getData()
24            #add to cache
25        else:           
26            valList=[]
27            try:
28                vals=gridOrd.coordAxisValues.coordinateList.CONTENT
29                for val in vals.split(): 
30                    valList.append(eval(val))
31            except:
32                vals=gridOrd.coordAxisValues.timePositionList.CONTENT
33                for val in vals.split(): 
34                    valList.append(val)           
35            self.domain[name]=valList   
36    return self.domain   
37
38   
39def subsetToTrajectory(self, outputdir=None, ncname='trajectory.nc' ,times=None, **kwargs):
40    ''' subsetting a trajectory feature by time '''
41   
42    #work out where to store the resulting files:
43    if outputdir is not None:
44        self.outputdir=outputdir
45    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
46        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
47     
48    if self.outputdir is not None:
49        pathToSubsetNetCDF=self.outputdir + '/' +ncname
50    else: 
51        pathToSubsetNetCDF=ncname
52    if self.outputdir is not None:
53        csml.csmllibs.csmlextra.checkDirExists(self.outputdir)
54 
55    #make sure single time string is converted to a tuple
56    if type(times) is str:
57        selection = tuple(csml.csmllibs.csmlextra.listify(times))
58    else:
59        selection = times 
60       
61       
62    #repurpose the existing trajectoryDomain for the output csml.
63    td=copy.deepcopy(self.value.trajectoryDomain)
64
65     
66    #get the requested times, and also max/min index values, needed to get the corresponding spatial axes 
67    tlist=[]
68    for gridOrd in td.coordTransformTable.gridOrdinates:
69        if hasattr(gridOrd.coordAxisValues, 'timePositionList'):
70            for val in gridOrd.coordAxisValues.timePositionList.CONTENT.split():
71                tlist.append(val)       
72            indices=[]
73            if  len(selection) ==2:       
74                selectedTimes=[]
75                for domaintime in tlist:
76                    if selection[0] == domaintime:
77                        minIndex=tlist.index(domaintime)
78                        selectedTimes.append(domaintime)
79                    if domaintime > selection[0]:
80                        if domaintime < selection[1]:
81                            selectedTimes.append(domaintime)         
82                    if selection[1] == domaintime:
83                        maxIndex=tlist.index(domaintime)
84                        selectedTimes.append(domaintime)         
85                        break   
86                           
87            elif len(selection)==1:   #point only
88                for domaintime in tlist:
89                    if selection[0] == domaintime:               
90                        minIndex=tlist.index(domaintime)
91                        maxIndex=minIndex
92                        selectedTimes=selection
93                        break
94            if len(selectedTimes) ==1:
95                timestr=selectedTimes[0]
96            else:
97                timestr=''
98                for t in selectedTimes:
99                    timestr = timestr + ' ' + t
100            #now substitute these times into the feature as we will use it later to write out a new feature.
101            gridOrd.coordAxisValues.timePositionList=csml.parser.csString(timestr)
102            gridOrd.coordAxisValues.id = csml.csmllibs.csmlextra.getRandomID()
103            break
104   
105    #change the high limit  (assumes low limit is 0)
106    highlim=''
107    for item in td.limits.high.CONTENT.split():
108        highlim=highlim+ ' ' + str(len(selectedTimes)-1)
109    td.limits.high.CONTENT=highlim   
110   
111    #Now get the data from the spatial axes using the indices
112    spatialAxes={}
113    for gridOrd in td.coordTransformTable.gridOrdinates:
114        if hasattr(gridOrd.coordAxisValues, 'coordinateList'):
115            data, fillvalue=gridOrd.coordAxisValues.insertedExtract.getDataFromChunks(minIndex, maxIndex)
116            if len(data)==1:
117                datastr=data[0]
118            else:
119                datastr=''
120                for d in data:
121                    datastr=datastr + ' ' + str(d)
122            gridOrd.coordAxisValues.coordinateList=csml.parser.csString(datastr)
123            gridOrd.coordAxisValues.id = csml.csmllibs.csmlextra.getRandomID()
124            spatialAxes[gridOrd.coordAxisLabel.CONTENT]=data       
125           
126            #get rid of the old xlink attributes
127            delattr(gridOrd.coordAxisValues, 'href')
128            delattr(gridOrd.coordAxisValues, 'role')
129            delattr(gridOrd.coordAxisValues, 'arcrole')
130            delattr(gridOrd.coordAxisValues, 'show')
131            delattr(gridOrd.coordAxisValues, '{http://www.w3.org/1999/xlink}href')
132            delattr(gridOrd.coordAxisValues, '{http://www.w3.org/1999/xlink}role')
133            delattr(gridOrd.coordAxisValues, '{http://www.w3.org/1999/xlink}arcrole')
134            delattr(gridOrd.coordAxisValues, '{http://www.w3.org/1999/xlink}show')
135
136           
137           
138    #Now create the file extract
139    fileList=[]
140    try:
141        fileList.append(self.value.rangeSet.valueArray.valueComponent.insertedExtract.fileName.CONTENT)
142        fextract=self.value.rangeSet.valueArray.valueComponent.insertedExtract
143        units=self.value.rangeSet.valueArray.valueComponent.uom
144    except:
145        for f in self.value.rangeSet.arrayDescriptor.components.fileList.fileNames.CONTENT.split():
146            fileList.append(f)
147            fextract=self.value.rangeSet.arrayDescriptor.components
148            units= self.value.rangeSet.arrayDescriptor.uom.CONTENT
149   
150    #Now write out the CSML Feature:
151    #create a TrajectoryCoverage
152    cvg=csml.parser.TrajectoryCoverage()
153    cvg.id=csml.csmllibs.csmlextra.getRandomID()
154   
155    td.id=csml.csmllibs.csmlextra.getRandomID()
156               
157    sdid=csml.csmllibs.csmlextra.getRandomID() # for the storage descriptor
158    #create a RangeSet
159    rs=csml.parser.RangeSet()
160    va=csml.parser.ValueArray()
161    vc=csml.parser.MeasureOrNullList()
162    vc.href='#%s'%sdid
163    vc.arcrole="http://ndg.nerc.ac.uk/xlinkUsage/insert#QuantityList"
164    vc.role="http://ndg.nerc.ac.uk/fileFormat/csmlStorageDescriptor"
165    vc.show='embed'
166    try:
167        vc.uom=self.value.rangeSet.valueArray.valueComponent.uom
168    except:
169        vc.uom = 'unknown'  #TODO,
170    va.valueComponent=vc
171    va.id=csml.csmllibs.csmlextra.getRandomID()
172    rs.valueArray=va
173   
174    #Add the domain and rangeSet as attributes of the coverage
175    cvg.trajectoryDomain=td
176    cvg.rangeSet=rs
177   
178    totalArraySize=10
179    descriptor=csml.parser.NetCDFExtract(id=sdid,fileName=csml.parser.csString(ncname),variableName=self.name,arraySize=csml.parser.csString(totalArraySize))
180
181   
182    #the parameter of the feature is of type Phenomenon, here href creates "xlink:href=..."
183    param=self.parameter
184   
185   
186    #create a stand alone pointseries feature containing this coverage
187    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
188    subsettedFeature=csmlWrap.createTrajectoryFeature(value=cvg,parameter=param,featureID=csml.csmllibs.csmlextra.getRandomID(),description=self.description)
189   
190     
191   
192    #Now open a new NetCDF File:   
193    nc=cdms.open(pathToSubsetNetCDF,'w')
194   
195    #create the variable
196    var=MV.array(data)
197    var.id=self.name.CONTENT
198    var.name=self.name.CONTENT
199    var.units=units[0] # hopefully there is just one unit in this list..
200
201   
202    #create the time axis
203    floatTimes=[]
204    tOne=csml.csmllibs.csmltime.getCDtime(selection[0])
205    tbase=csml.csmllibs.csmltime.getBaseUnits(tOne)
206    for time in selectedTimes:
207        time=csml.csmllibs.csmltime.getCDtime(time).torel(tbase)
208        floatTimes.append(time.value)
209    timeAx=cdms.createAxis(floatTimes)
210    timeAx.designateTime()
211    timeAx.id='time'
212    timeAx.standard_name='time'
213    timeAx.units=tbase 
214    var.setAxis(0,timeAx)     
215    nc.write(var)
216   
217   
218    #get standard_name table
219    sntable=csml.csmllibs.standardnames.__path__[0]+'/cf-standard-name-table.xml'
220    snt=csml.csmllibs.standardnames.snames.StandardNameTable(tablelocation=sntable)
221   
222    #add auxiliary coordinate variables to hold the locations
223    for ax in spatialAxes:
224        var=MV.array(spatialAxes[ax])
225        var.id=ax
226        #if id is standard_name, use it as a standard_name       
227        if snt.lookup(ax) is 1:
228            var.standard_name=ax
229        var.setAxis(0,timeAx)
230        nc.write(var)   
231               
232   
233    nc.close()   
234    print 'NetCDF file written to %s'%pathToSubsetNetCDF
235    return subsettedFeature, pathToSubsetNetCDF, descriptor
236
237   
238   
Note: See TracBrowser for help on using the repository browser.