source: TI02-CSML/branches/csml-cdms2/API/ops_TrajectoryFeature.py @ 3627

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

This branch contains CSML converted to use cdat_lite-5.

  • convertcdms was run on the source
  • the MA.set_print_limit call was changed to the numpy equivilent
  • The tests were changed to account for the existence of numpy scalar types.

All tests, except the two known to fail, pass on i686 ubuntu.

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