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

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

More docstrings

Line 
1''' ops_RaggedSectionFeature  contains operations for RaggedSectionFeatures. These methods are attached to GridFeature instances at runtime.'''
2import csml
3import MA,MV,Numeric
4import cdms
5
6
7def testmethod(self):
8    ''' test method'''
9    print 'testmethod for ragged section feature'
10    return 'testmethod 2'
11
12def getDomain(self):
13    '''returns domain as a dictionary of ordinates {name: [values], ...}'''
14    self.domain={}
15   
16    #get lat lon station locations
17    locations= self.stationLocations.CONTENT.split()
18    latitude=[]
19    longitude=[]
20    for i in range(0, len(locations)-1,2):
21        latitude.append(locations[i])
22        longitude.append(locations[i+1])
23       
24    self.domain['latitude']=latitude
25    self.domain['longitude']=longitude
26   
27    #get station times
28    times=[]
29    for t in self.stationTimes.CONTENT.split():
30        times.append(t)
31    self.domain['time']=times
32
33    #get vertical dimension
34    gridOrd= self.value.sectionDomain.coordTransformTable.gridOrdinates
35    name=gridOrd.coordAxisLabel.CONTENT
36    if hasattr(gridOrd.coordAxisValues, 'insertedExtract'):
37        self.domain[name], fill, axisorder, units=gridOrd.coordAxisValues.insertedExtract.getData()
38    else:           
39        valList=[]
40        vals=gridOrd.coordAxisValues.coordinateList.CONTENT
41        for val in vals.split(): 
42            valList.append(eval(val))           
43        self.domain[name]=valList   
44    self._depths=self.domain[name]
45    self._depthname=name
46           
47    return self.domain
48       
49
50def subsetByTime(self, outputdir=None, ncname='pointseries.nc' ,time=None):
51    ''' subset a ragged section by time
52    @param outputdir:    name of directory to write data files to
53    @param ncname:     name of netcdf file to write out
54    @param time:      subset selection by time.
55    @return:     subsettedFeature (RaggedSectionFeature)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)
56    '''
57   
58    descriptors=[]
59    if outputdir is not None:
60        self.outputdir=outputdir
61    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
62        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
63     
64    if self.outputdir is not None:
65        pathToSubsetNetCDF=self.outputdir + '/' +ncname
66    else: 
67        pathToSubsetNetCDF=ncname
68    if self.outputdir is not None:
69        csml.csmllibs.csmlextra.checkDirExists(self.outputdir)
70    if type(time) is str:
71        selection = tuple(csml.csmllibs.csmlextra.listify(time))
72    else:
73        selection = time 
74    tlist=[]   
75    for val in self.stationTimes.CONTENT.split():
76        tlist.append(val)
77   
78    indices=[]
79    if  len(selection) ==2:       
80        selectedTimes=[]
81        for domaintime in tlist:
82            if selection[0] == domaintime:
83                minIndex=tlist.index(domaintime)
84                selectedTimes.append(domaintime)
85            if domaintime > selection[0]:
86                if domaintime < selection[1]:
87                    selectedTimes.append(domaintime)         
88            if selection[1] == domaintime:
89                maxIndex=tlist.index(domaintime)
90                selectedTimes.append(domaintime)         
91                break   
92                   
93    elif len(selection)==1:   #point only
94        for domaintime in tlist:
95            if selection[0] == domaintime:             
96                minIndex=tlist.index(domaintime)
97                maxIndex=minIndex
98                selectedTimes=selection
99                break
100           
101    if not hasattr(self, 'domain'):
102        self.domain=self.getDomain()
103    #need to get the lat and lon to correspond with the selected times:
104    selectedLats=[eval(lat) for lat in self.domain['latitude'][minIndex:maxIndex+1]]
105    selectedLons= [eval(lon) for lon in self.domain['longitude'][minIndex:maxIndex+1]]
106   
107    #get maximum profile length:
108    lengths=self.profileLength.CONTENT.split()
109    maxlen=0
110    for l in lengths:
111        testlen=eval(l)
112        if testlen > maxlen:
113            maxlen =testlen
114    #and store selected lengths:
115    selectedLengths=[l for l in lengths[minIndex:maxIndex+1]]       
116   
117
118    #assuming each timestep is in a separate file, minIndex and maxIndex refer to the position of files in the series.       
119    fulldata=[]
120    fileList=[]
121    filenames=self.value.rangeSet.arrayDescriptor.components.fileList.fileNames.CONTENT.split() 
122    for i in range(minIndex, maxIndex+1):
123        fileList.append(filenames[i])
124        data,fillvalue, axisorder, units=self.value.rangeSet.arrayDescriptor.components.getData(fileposition=i)       
125        if len(data)<maxlen:
126            data.resize(maxlen)
127        fulldata.append(data)
128           
129   
130    #Write out CSML feature
131    #create a RaggedSectionCoverage
132    cvg=csml.parser.SectionCoverage()
133    cvg.id=csml.csmllibs.csmlextra.getRandomID()
134   
135     
136    #create a Domain (appropriate domain for a PointSeriesFeature
137    sect=csml.parser.SectionDomain()
138    sect.id=csml.csmllibs.csmlextra.getRandomID()
139   
140    #get min and max profile lengths:
141   
142    minl = maxl=eval(selectedLengths[0])
143    for l in selectedLengths:
144        testl=eval(l)
145        if testl < minl:
146            minl=testl
147        if testl > maxl:
148            maxl=testl
149       
150    lowlim=csml.parser.csString(minl)
151    highlim=csml.parser.csString(maxl)
152    sect.limits=csml.parser.GridEnvelope(low=lowlim, high=highlim)
153    sect.srsName=self.value.sectionDomain.srsName
154    sect.axisLabels=self.value.sectionDomain.axisLabels
155    sect.dimension=self.value.sectionDomain.dimension
156    sect.srsDimension=self.value.sectionDomain.srsDimension
157    sect.aLabels=self.value.sectionDomain.aLabels
158   
159    sect.coordTransformTable=self.value.sectionDomain.coordTransformTable
160    descriptorid=csml.csmllibs.csmlextra.getRandomID()
161    sect.coordTransformTable.gridOrdinates.coordAxisValues.href='#'+descriptorid
162    vname=sect.coordTransformTable.gridOrdinates.coordAxisLabel
163    descriptors.append(csml.parser.NetCDFExtract(id=descriptorid,fileName=csml.parser.csString(ncname),variableName=vname,arraySize=csml.parser.csString(len(self._depths))))
164       
165    sdid=csml.csmllibs.csmlextra.getRandomID() # for the storage descriptor
166    #create a RangeSet
167    rs=csml.parser.RangeSet()
168    va=csml.parser.ValueArray()
169    vc=csml.parser.MeasureOrNullList()
170    vc.href='#%s'%sdid
171    vc.arcrole="http://ndg.nerc.ac.uk/xlinkUsage/insert#QuantityList"
172    vc.role="http://ndg.nerc.ac.uk/fileFormat/csmlStorageDescriptor"
173    vc.show='embed'
174   
175   
176    try:
177        vc.uom=self.value.rangeSet.valueArray.valueComponent.uom
178    except:
179        try:
180            vc.uom=self.value.rangeSet.arrayDescriptor.uom.CONTENT
181        except:
182            vc.uom = 'unknown' 
183    va.valueComponent=vc
184    va.id=csml.csmllibs.csmlextra.getRandomID()
185    rs.valueArray=va
186   
187   
188    #Add the domain and rangeSet as attributes of the coverage
189    cvg.sectionDomain=sect
190    cvg.rangeSet=rs
191   
192    descriptors.append(csml.parser.NetCDFExtract(id=sdid,fileName=csml.parser.csString(ncname),variableName=self.name,arraySize=csml.parser.csString(len(selectedTimes))))
193   
194    #the parameter of the feature is of type Phenomenon, here href creates "xlink:href=..."
195    param=self.parameter
196   
197     
198    #create a stand alone raggedsection feature containing this coverage
199    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
200    subsettedFeature=csmlWrap.createRaggedSectionFeature(value=cvg,parameter=param, stationtimes=selectedTimes, stationlats=selectedLats, stationlons=selectedLons, profilelengths=selectedLengths,featureID=csml.csmllibs.csmlextra.getRandomID(),name=self.name,description=self.description)
201   
202   
203       
204    #Write out NetCDF file
205        #Now write out the NetCDF File:   
206   
207    nc=cdms.open(pathToSubsetNetCDF,'w')
208   
209    #create profile axis
210    profAx=cdms.createAxis(range(0,len(selectedTimes)))
211    profAx.id='profile'
212   
213    #create depths/pressure/height axis and variable
214    depthAx=cdms.createAxis(self._depths)
215    depthAx.id=self._depthname
216   
217    #create variable for feature
218    var=MV.array(fulldata)
219    var.id=self.name.CONTENT
220    var.name=self.name.CONTENT
221    var.units=units[0] # hopefully there is just one unit in this list..
222    setattr(var, 'missing_value' ,fillvalue)
223    var.setAxis(0,profAx)       
224    var.setAxis(1,depthAx)   
225    setattr(var, 'standard_name',param.getStandardName())
226    nc.write(var)
227   
228    #create variables for latitude/longitude
229    varlat = MV.array(selectedLats)
230    varlat.id='latitude'
231    varlat.name='latitude'
232    varlat.setAxis(0,profAx)
233    nc.write(varlat)
234   
235    varlon = MV.array(selectedLons)
236    varlon.id='longitude'
237    varlon.name='longitude'
238    varlon.setAxis(0,profAx)
239    nc.write(varlon)
240   
241    #create variable for times
242    floatTimes=[]
243   
244    tOne=csml.csmllibs.csmltime.getCDtime(selection[0])
245    tbase=csml.csmllibs.csmltime.getBaseUnits(tOne)
246    for time in selectedTimes:
247        time=csml.csmllibs.csmltime.getCDtime(time).torel(tbase)
248        floatTimes.append(time.value)
249   
250    vart=MV.array(floatTimes)
251    vart.id='time'
252    vart.units=tbase
253    vart.setAxis(0,profAx)
254    nc.write(vart)
255   
256    nc.close()   
257    print 'NetCDF file written to %s'%pathToSubsetNetCDF
258    return subsettedFeature, pathToSubsetNetCDF, descriptors
Note: See TracBrowser for help on using the repository browser.