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

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

more on subsetting trajectory

Line 
1''' ops_TrajectoryFeature  contains operations for TrajectoryFeature'''
2
3import csml, cdms, MV
4import sys
5
6def testmethod(self):
7    print 'testmethod for TrajectoryFeature'
8    return 'testmethod - TrajectoryFeature'
9
10def getAllowedSubsettings(self):
11    return ['subsetToTrajectory'] 
12
13
14def getDomain(self):
15    #returns domain as a dictionary of ordinates {name: [values], ...}
16    self.domain={}
17    self.gridnames={}
18    for gridOrd in self.value.trajectoryDomain.coordTransformTable.gridOrdinates:
19        name=gridOrd.coordAxisLabel.CONTENT
20        self.gridnames[gridOrd.gridAxesSpanned.CONTENT]=name
21        if hasattr(gridOrd.coordAxisValues, 'insertedExtract'):
22            self.domain[name], fill, axisorder, units=gridOrd.coordAxisValues.insertedExtract.getData()
23            #add to cache
24        else:           
25            valList=[]
26            try:
27                vals=gridOrd.coordAxisValues.coordinateList.CONTENT
28                for val in vals.split(): 
29                    valList.append(eval(val))
30            except:
31                vals=gridOrd.coordAxisValues.timePositionList.CONTENT
32                for val in vals.split(): 
33                    valList.append(val)           
34            self.domain[name]=valList   
35    return self.domain   
36               
37   
38               
39   
40def subsetToTrajectory(self, outputdir=None, ncname='trajectory.nc' ,times=None):
41    ''' subsetting a trajectory feature by time '''
42   
43    #work out where to store the resulting files:
44    if outputdir is not None:
45        self.outputdir=outputdir
46    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
47        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
48     
49    if self.outputdir is not None:
50        pathToSubsetNetCDF=self.outputdir + '/' +ncname
51    else: 
52        pathToSubsetNetCDF=ncname
53    if self.outputdir is not None:
54        csml.csmllibs.csmlextra.checkDirExists(self.outputdir)
55 
56    #make sure single time string is converted to a tuple
57    if type(times) is str:
58        selection = tuple(csml.csmllibs.csmlextra.listify(times))
59    else:
60        selection = times 
61       
62       
63    #repurpose the existing trajectoryDomain for the output csml.
64    td=self.value.trajectoryDomain                         
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    setattr(var, 'missing_value' ,fillvalue)
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.units=tbase 
213    var.setAxis(0,timeAx)     
214    nc.write(var)
215   
216   
217    #add auxiliary coordinate variables to hold the locations
218    for ax in spatialAxes:
219        var=MV.array(spatialAxes[ax])
220        var.id=ax
221        var.setAxis(0,timeAx)
222        nc.write(var)   
223               
224   
225    nc.close()   
226    print 'NetCDF file written to %s'%pathToSubsetNetCDF
227    return subsettedFeature, pathToSubsetNetCDF, descriptor
228
229   
230   
Note: See TracBrowser for help on using the repository browser.