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

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

subsetting of profileseries, pointseries retested and fixed

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    if not hasattr(self, 'domain'):
94        self.domain=self.getDomain()
95    #need to get the lat and lon to correspond with the selected times:
96    selectedLats=[eval(lat) for lat in self.domain['latitude'][minIndex:maxIndex+1]]
97    selectedLons= [eval(lon) for lon in self.domain['longitude'][minIndex:maxIndex+1]]
98   
99    #get maximum profile length:
100    lengths=self.profileLength.CONTENT.split()
101    maxlen=0
102    for l in lengths:
103        testlen=eval(l)
104        if testlen > maxlen:
105            maxlen =testlen
106    #and store selected lengths:
107    selectedLengths=[l for l in lengths[minIndex:maxIndex+1]]       
108   
109
110    #assuming each timestep is in a separate file, minIndex and maxIndex refer to the position of files in the series.       
111    fulldata=[]
112    fileList=[]
113    filenames=self.value.rangeSet.arrayDescriptor.components.fileList.fileNames.CONTENT.split() 
114    for i in range(minIndex, maxIndex+1):
115        fileList.append(filenames[i])
116        data,fillvalue, axisorder, units=self.value.rangeSet.arrayDescriptor.components.getData(fileposition=i)       
117        if len(data)<maxlen:
118            data.resize(maxlen)
119        fulldata.append(data)
120           
121   
122    #Write out CSML feature
123    #create a RaggedSectionCoverage
124    cvg=csml.parser.SectionCoverage()
125    cvg.id=csml.csmllibs.csmlextra.getRandomID()
126   
127     
128    #create a Domain (appropriate domain for a PointSeriesFeature
129    sect=csml.parser.SectionDomain()
130    sect.id=csml.csmllibs.csmlextra.getRandomID()
131   
132    #get min and max profile lengths:
133   
134    minl = maxl=eval(selectedLengths[0])
135    for l in selectedLengths:
136        testl=eval(l)
137        if testl < minl:
138            minl=testl
139        if testl > maxl:
140            maxl=testl
141       
142    lowlim=csml.parser.csString(minl)
143    highlim=csml.parser.csString(maxl)
144    sect.limits=csml.parser.GridEnvelope(low=lowlim, high=highlim)
145    sect.srsName=self.value.sectionDomain.srsName
146    sect.axisLabels=self.value.sectionDomain.axisLabels
147    sect.dimension=self.value.sectionDomain.dimension
148    sect.srsDimension=self.value.sectionDomain.srsDimension
149    sect.aLabels=self.value.sectionDomain.aLabels
150   
151    sect.coordTransformTable=self.value.sectionDomain.coordTransformTable
152    descriptorid=csml.csmllibs.csmlextra.getRandomID()
153    sect.coordTransformTable.gridOrdinates.coordAxisValues.href='#'+descriptorid
154    vname=sect.coordTransformTable.gridOrdinates.coordAxisLabel
155    descriptors.append(csml.parser.NetCDFExtract(id=descriptorid,fileName=csml.parser.csString(ncname),variableName=vname,arraySize=csml.parser.csString(len(self._depths))))
156       
157    sdid=csml.csmllibs.csmlextra.getRandomID() # for the storage descriptor
158    #create a RangeSet
159    rs=csml.parser.RangeSet()
160    va=csml.parser.ValueArray()
161    vc=csml.parser.MeasureOrNullList()
162    vc.href='#%s'%sdid
163    vc.arcrole="http://ndg.nerc.ac.uk/xlinkUsage/insert#QuantityList"
164    vc.role="http://ndg.nerc.ac.uk/fileFormat/csmlStorageDescriptor"
165    vc.show='embed'
166   
167   
168    try:
169        vc.uom=self.value.rangeSet.valueArray.valueComponent.uom
170    except:
171        try:
172            vc.uom=self.value.rangeSet.arrayDescriptor.uom.CONTENT
173        except:
174            vc.uom = 'unknown' 
175    va.valueComponent=vc
176    va.id=csml.csmllibs.csmlextra.getRandomID()
177    rs.valueArray=va
178   
179   
180    #Add the domain and rangeSet as attributes of the coverage
181    cvg.sectionDomain=sect
182    cvg.rangeSet=rs
183   
184    descriptors.append(csml.parser.NetCDFExtract(id=sdid,fileName=csml.parser.csString(ncname),variableName=self.name,arraySize=csml.parser.csString(len(selectedTimes))))
185   
186    #the parameter of the feature is of type Phenomenon, here href creates "xlink:href=..."
187    param=self.parameter
188   
189     
190    #create a stand alone raggedsection feature containing this coverage
191    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
192    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)
193   
194   
195       
196    #Write out NetCDF file
197        #Now write out the NetCDF File:   
198   
199    nc=cdms.open(pathToSubsetNetCDF,'w')
200   
201    #create profile axis
202    profAx=cdms.createAxis(range(0,len(selectedTimes)))
203    profAx.id='profile'
204   
205    #create depths/pressure/height axis and variable
206    depthAx=cdms.createAxis(self._depths)
207    depthAx.id=self._depthname
208   
209    #create variable for feature
210    var=MV.array(fulldata)
211    var.id=self.name.CONTENT
212    var.name=self.name.CONTENT
213    var.units=units[0] # hopefully there is just one unit in this list..
214    setattr(var, 'missing_value' ,fillvalue)
215    var.setAxis(0,profAx)       
216    var.setAxis(1,depthAx)   
217    setattr(var, 'standard_name',param.getStandardName())
218    nc.write(var)
219   
220    #create variables for latitude/longitude
221    varlat = MV.array(selectedLats)
222    varlat.id='latitude'
223    varlat.name='latitude'
224    varlat.setAxis(0,profAx)
225    nc.write(varlat)
226   
227    varlon = MV.array(selectedLons)
228    varlon.id='longitude'
229    varlon.name='longitude'
230    varlon.setAxis(0,profAx)
231    nc.write(varlon)
232   
233    #create variable for times
234    floatTimes=[]
235   
236    tOne=csml.csmllibs.csmltime.getCDtime(selection[0])
237    tbase=csml.csmllibs.csmltime.getBaseUnits(tOne)
238    for time in selectedTimes:
239        time=csml.csmllibs.csmltime.getCDtime(time).torel(tbase)
240        floatTimes.append(time.value)
241   
242    vart=MV.array(floatTimes)
243    vart.id='time'
244    vart.units=tbase
245    vart.setAxis(0,profAx)
246    nc.write(vart)
247   
248    nc.close()   
249    print 'NetCDF file written to %s'%pathToSubsetNetCDF
250    return subsettedFeature, pathToSubsetNetCDF, descriptors
Note: See TracBrowser for help on using the repository browser.