source: TI02-CSML/trunk/csml/API/ops_PointSeriesFeature.py @ 3530

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

fixed two small bugs in pointseries subsetting code

Line 
1''' ops_PointSeriesFeature  contains operations for PointSeriesFeatures. These methods are attached to PointSeriesFeature instances at runtime.'''
2import csml,sys
3import cdms, MV, MA
4
5def testmethod(self):
6    ''' testmethod for gridseries feature'''
7    print 'testmethod for pointSeries feature'
8    return 'testmethod pointSeries'
9
10def getAllowedSubsettings(self):
11    '''get the allowed operations
12    @return:    list of operations'''
13    return ['subsetToPointSeries']
14   
15def getDomain(self):
16    ''' returns domain as a dictionary of ordinates {name: [values], ...} '''
17    self.domain={}
18    valList=[]
19    for val in self.value.pointSeriesDomain.timePositionList.CONTENT.split():
20        valList.append(val)
21    self.domain['time']=valList
22    [lat,lon]=self.location.CONTENT.split()
23    self.domain['longitude']=lon
24    self.domain['latitude']=lat
25
26    return self.domain   
27   
28def getLongitudeAxis(self):
29    ''' get the index of the longitude axis'''
30    return 'longitude' #dummy value not used, subsetting is by time only
31
32def getLatitudeAxis(self):
33    ''' get the index of the latitude axis'''
34    return 'latitude' #dummy value not used,  subsetting is by time only
35
36def getTimeAxis(self):
37    ''' get the index of the time axis'''
38    return 'times'
39   
40def subsetToPointSeries(self, outputdir=None, ncname='pointseries.nc' ,times=None):
41    '''Subset a PointSeriesFeature to a PointSeriesFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs.
42   
43    @param outputdir:    name of directory to write data files to
44    @param ncname:     name of netcdf file to write out
45    @param times:       selection by times
46    @return:     subsettedFeature (PointSeriesFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)'''
47    if outputdir is not None:
48        self.outputdir=outputdir
49    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
50        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
51     
52    if self.outputdir is not None:
53        pathToSubsetNetCDF=self.outputdir + '/' +ncname
54    else: 
55        pathToSubsetNetCDF=ncname
56    if self.outputdir is not None:
57        csml.csmllibs.csmlextra.checkDirExists(self.outputdir)
58 
59    if type(times) is str:
60        selection = tuple(csml.csmllibs.csmlextra.listify(times))
61    else:
62        selection = times 
63    tlist=[]   
64    for val in self.value.pointSeriesDomain.timePositionList.CONTENT.split():
65        tlist.append(val)
66   
67    indices=[]
68    if  len(selection) ==2:       
69        selectedTimes=[]
70        for domaintime in tlist:
71            if selection[0] == domaintime:
72                minIndex=tlist.index(domaintime)
73                selectedTimes.append(domaintime)
74            if domaintime > selection[0]:
75                if domaintime < selection[1]:
76                    selectedTimes.append(domaintime)         
77            if selection[1] == domaintime:
78                maxIndex=tlist.index(domaintime)
79                selectedTimes.append(domaintime)         
80                break   
81                   
82    elif len(selection)==1:   #point only
83        for domaintime in tlist:
84            if selection[0] == domaintime:
85               
86                minIndex=tlist.index(domaintime)
87                maxIndex=minIndex
88                selectedTimes=selection
89                break
90   
91   
92   
93    #fileList=[]
94       
95    if hasattr(self.value.rangeSet, 'quantityList'):       
96        #inline
97        data=self.value.rangeSet.quantityList.CONTENT.split()[minIndex:maxIndex+1]
98        fillvalue=None
99        units=self.value.rangeSet.quantityList.uom
100    else:
101        #file extract
102        if hasattr(self.value.rangeSet, 'valueArray'):
103            fextract=self.value.rangeSet.valueArray.valueComponent.insertedExtract.components
104            uom = self.value.rangeSet.valueArray.valueComponent.insertedExtract.uom
105        else:
106            fextract=self.value.rangeSet.arrayDescriptor.components
107            uom = self.value.rangeSet.arrayDescriptor.uom
108        #for f in fextract.fileList.fileNames.CONTENT.split():
109            #fileList.append(f)
110       
111        data, fillvalue=fextract.getDataFromChunks(minIndex, maxIndex)
112        units=[uom.CONTENT]
113   
114    #Now write out the CSML Feature:
115    #create a PointSeriesCoverage
116    cvg=csml.parser.PointSeriesCoverage()
117    cvg.id=csml.csmllibs.csmlextra.getRandomID()
118   
119    #create a TimeSeriesDomain (appropriate domain for a PointSeriesFeature
120    ptsd=csml.parser.TimeSeries()
121    ptsd.id=csml.csmllibs.csmlextra.getRandomID()
122    #create (and populate) a TimePositionList. Using keyword arguements for conciseness.
123    ptsd.timePositionList=csml.parser.TimePositionList(CONTENT=csml.csmllibs.csmlextra.stringify(selectedTimes))
124       
125    sdid=csml.csmllibs.csmlextra.getRandomID() # for the storage descriptor
126    #create a RangeSet
127    rs=csml.parser.RangeSet()
128    va=csml.parser.ValueArray()
129    vc=csml.parser.MeasureOrNullList()
130    vc.href='#%s'%sdid
131    vc.arcrole="http://ndg.nerc.ac.uk/xlinkUsage/insert#QuantityList"
132    vc.role="http://ndg.nerc.ac.uk/fileFormat/csmlStorageDescriptor"
133    vc.show='embed'
134    try:
135        vc.uom=self.value.rangeSet.valueArray.valueComponent.uom
136    except:
137        vc.uom = 'unknown'  #TODO,
138    va.valueComponent=vc
139    va.id=csml.csmllibs.csmlextra.getRandomID()
140    rs.valueArray=va
141   
142   
143    #Add the domain and rangeSet as attributes of the coverage
144    cvg.pointSeriesDomain=ptsd
145    cvg.rangeSet=rs
146   
147    totalArraySize=10
148    if hasattr(self, 'name'):
149        varname=self.name.CONTENT
150    else:
151        varname=self.id
152    descriptor=csml.parser.NetCDFExtract(id=sdid,fileName=csml.parser.csString(ncname),variableName=csml.parser.csString(varname),arraySize=csml.parser.csString(totalArraySize))
153
154    #the parameter of the feature is of type Phenomenon, here href creates "xlink:href=..."
155    param=self.parameter
156   
157     
158    #create a stand alone pointseries feature containing this coverage
159    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
160    subsettedFeature=csmlWrap.createPointSeriesFeature(value=cvg,parameter=param,location=self.location,featureID=self.id,description=self.description)
161     
162   
163    #Now write out the NetCDF File:   
164    nc=cdms.open(pathToSubsetNetCDF,'w')
165    var=MV.array(data)
166    var.id=varname
167    var.name=varname
168    var.units=units[0] # hopefully there is just one unit in this list..
169    setattr(var, 'missing_value' ,fillvalue)
170   
171    floatTimes=[]
172   
173    tOne=csml.csmllibs.csmltime.getCDtime(selection[0])
174    tbase=csml.csmllibs.csmltime.getBaseUnits(tOne)
175    for time in selectedTimes:
176        time=csml.csmllibs.csmltime.getCDtime(time).torel(tbase)
177        floatTimes.append(time.value)
178    timeAx=cdms.createAxis(floatTimes)
179    timeAx.designateTime()
180    timeAx.id='time'
181    timeAx.units=tbase 
182    var.setAxis(0,timeAx)     
183    nc.write(var)
184    nc.close()   
185    print 'NetCDF file written to %s'%pathToSubsetNetCDF
186    return subsettedFeature, pathToSubsetNetCDF, descriptor
187   
188   
189def subsetToPoint(self, outputdir=None, ncname='point.nc',time=None):
190    '''Subset a PointSeriesFeature to a PointFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs.
191       
192    @param outputdir:    name of directory to write data files to
193    @param ncname:     name of netcdf file to write out
194    @param time:       selection by single time
195    @return:     subsettedFeature (PointFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)'''
196    if outputdir is not None:
197        self.outputdir=outputdir
198    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
199        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
200    pseriesFeature, pathToSubsetNetCDF, descriptor=self.subsetToPointSeries(self.outputdir, ncname,time)         
201    #now need to take this point series feature containing one time and rewrite it as a point feature.
202    pointFeature=csml.parser.PointFeature()
203    pointFeature.location=pseriesFeature.location
204    pointFeature.parameter=pseriesFeature.parameter
205    pointFeature.description=pseriesFeature.description
206    cvg=csml.parser.PointCoverage()
207    domain=csml.parser.PointDomain()
208    domain.id=csml.csmllibs.csmlextra.getRandomID()
209    domain.pointMember=csml.parser.csString(pseriesFeature.location.CONTENT)
210    #create (and populate) a TimePositionList. Using keyword arguements for conciseness.
211    cvg.pointDomain=domain
212    pointFeature.value=cvg
213    pointFeature.time=csml.parser.csString(time)
214    return pointFeature, pathToSubsetNetCDF, descriptor
Note: See TracBrowser for help on using the repository browser.