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

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

correct url for cf converions

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