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

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

more on Trajectory

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