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

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

Any trajectory data outside selected bounding box is set to be missing value using masked array

Line 
1''' ops_TrajectoryFeature  contains operations for TrajectoryFeature'''
2
3import csml, cdms, MV, MA
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,mask=self.m)
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                            dataToRemove.append(item[0])
204                            continue
205                        elif item[1] >maxval:
206                            dataToRemove.append(item[0])
207            self.spatialAxes[gridOrd.coordAxisLabel.CONTENT]=data           
208            #get rid of the old xlink attributes
209            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']:
210                try:
211                    delattr(gridOrd.coordAxisValues, att)
212                except:
213                    pass
214   
215    #create mask to mask out values outside any bounding box selection
216    maskIndices= list(set(dataToRemove))
217    self.m=MA.create_mask((len(data),))
218    for i in maskIndices:
219        self.m[i]=1   
220    #now get the actual rangeset data
221    if hasattr(self.value.rangeSet.valueArray.valueComponent, 'insertedExtract'):
222        self.data,self.fillvalue=self.value.rangeSet.valueArray.valueComponent.insertedExtract.getDataFromChunks(minIndex, maxIndex)
223    else:
224        print 'warning: inline rangeset not implemented'
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.