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

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

More docstrings

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