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

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

better handling of arrays to speed up subsetting

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