source: TI02-CSML/trunk/csml/API/ops_GridSeriesFeature.py @ 4019

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

more changes to allow preservation of netcdf attributes. Not complete, but not broken

Line 
1''' ops_GridSeriesFeature  contains operations for GridSeriesFeatures. These methods are attached to GridFeature instances at runtime.'''
2
3import csml.parser
4import csml.csmllibs.csmltime
5import csml.csmllibs.csmldocument
6import csml.API.ops_AbstractFeature
7import csml.API.genSubset
8import csml.csmllibs.netCDFWriter
9import csml.csmllibs.csmlcrs
10import csmlutils
11import os
12
13def testmethod(self):
14    ''' testmethod for gridseries feature'''
15    return 'testmethod - gridseries'
16
17def getAllowedSubsettings(self):
18    '''get the allowed operations
19    @return:    list of operations'''
20    return ['subsetToGridSeries', 'subsetToProfileSeries']  #other operations
21
22
23def getBoundingBox(self):
24    ''' get the simple boundingBox '''
25    try:
26        lower1=eval(bbox.boundedBy.lowerCorner.CONTENT.split()[0])
27        lower2=eval(bbox.boundedBy.lowerCorner.CONTENT.split()[1])
28        upper1=eval(bbox.boundedBy.upperCorner.CONTENT.split()[0])
29        upper2=eval(bbox.boundedBy.upperCorner.CONTENT.split()[1])
30        return [lower1,lower2,upper1,upper2]
31    except:
32        return None
33
34def _preserveAttributes(self):
35    '''  attempts to preserve information from the original file if information exist, such as CF attributes, projection info etc
36    returns?? tbd'''   
37    filename = self.value.rangeSet.valueArray.valueComponent.insertedExtract.components.fileList.fileNames.CONTENT.split(' ')[0]
38    di = csml.csmllibs.csmldataiface.DataInterface()
39    di = di.getUnknownInterfaceType(filename)
40    di.openFile(filename)
41    preservedData=csml.csmllibs.csmlpreservation.PreservedData(self.name.CONTENT, di)
42    return preservedData
43   
44def getDomain(self):
45    ''' get domain of the feature
46    @return:   domain as a dictionary of ordinates {name: [values], ...}'''
47    self.domain={}
48    self.gridnames={}
49    for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
50        name=gridOrd.coordAxisLabel.CONTENT
51        self.gridnames[gridOrd.gridAxesSpanned.CONTENT]=name
52        if hasattr(gridOrd.coordAxisValues, 'insertedExtract'):
53            self.domain[name], fill, axisorder, units=gridOrd.coordAxisValues.insertedExtract.getData()
54            #add to cache
55        else:           
56            valList=[]
57            try:
58                vals=gridOrd.coordAxisValues.coordinateList.CONTENT
59                for val in vals.split(): 
60                    valList.append(eval(val))
61            except:
62                vals=gridOrd.coordAxisValues.timePositionList.CONTENT
63                for val in vals.split(): 
64                    valList.append(val)           
65            self.domain[name]=valList   
66    return self.domain
67
68def getNativeCRS(self):
69    ''' get the the coordinate reference system of the domain of the feature
70    return:   a string, the srsName'''
71    return self.value.gridSeriesDomain.srsName
72
73def getUom(self):
74    ''' get the uom of the phenomenon.
75    @return: string - uom attribute'''
76    uom=None
77    try: 
78        uom=self.value.rangeSet.valueArray.valueComponent.quantityList.insertedExtract.uom.CONTENT
79    except AttributeError:
80        uom =None
81    return uom
82   
83def getDomainUnits(self):
84    ''' get the units used in the domain
85    @return:     list of units'''   
86    #srsname=self.value.gridSeriesDomain.srsName
87    #cat=csml.csmllibs.csmlcrs.CRSCatalogue()
88    #crs=cat.getCRS(srsname)
89    #return crs.units
90    return self.value.gridSeriesDomain.uomLabels.split()
91   
92def __getAxis(self, name):
93    '''called by getLongitudeAxis, getTimeAxis and getLatitudeAxis'''   
94    srsname=self.value.gridSeriesDomain.srsName
95    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
96    crs=cat.getCRS(srsname) 
97    axID=None
98    if name == 'lon':
99        axID = crs.lonAxis
100    elif name == 'lat':
101        axID = crs.latAxis
102    elif name == 'time':
103        axID = crs.timeAxis     
104    if axID is not None:
105        return crs.getAxisLabels()[axID]
106    else:
107        return None
108   
109def getLongitudeAxis(self):
110    ''' get the index of the longitude axis'''
111    return self.__getAxis('lon')
112
113def getLatitudeAxis(self):
114    ''' get the index of the latitude axis'''
115    return self.__getAxis('lat')
116
117def getTimeAxis(self):
118    ''' get the index of the time axis'''
119    return self.__getAxis('time')
120
121def getAxisLabels(self):
122    '''get the names of the axes'''
123    labels=[]
124    for label in self.value.gridSeriesDomain.axisLabels.split():
125        labels.append(label)
126    return labels
127
128def _subsetGrid(self, **kwargs):
129    '''this takes a selection from a gridseries object (be it a profile, a grid, whatever, it is still just a selection).  It is called internally  by the subsetToGridSeries method.
130    @return: a load of things which can then be used to rebuild a new csml feature: pathToSubsetNetCDF, crs, frame, timeName, times, strTimes, caltype, axisorder,units, fulldata, fillvalue, kwargs '''
131           
132    #set self.domain:
133    self.getDomain()     
134   
135    #get the CRS from a  the  catalogue
136    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
137    crs=cat.getCRS(self.value.gridSeriesDomain.srsName) 
138   
139    #non-feature specific setup code, mainly handles the time dimension/calendar
140   
141    pathToSubsetNetCDF, kwargs, timeAxis, timeName, caltype, times=csml.API.genSubset.genericSubset(self, self.outputdir, self.ncname, self.domain, kwargs)
142    ##Get names of variables in file and relate them to the subset selection
143    selection={}
144    frame=''
145   
146    for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
147        try:
148            selection[gridOrd.gridAxesSpanned.CONTENT]=kwargs[gridOrd.coordAxisLabel.CONTENT]           
149        except KeyError:
150            allValues=tuple(self.domain[gridOrd.coordAxisLabel.CONTENT])
151        if hasattr(gridOrd.coordAxisValues, 'timePositionList'):
152            if hasattr(gridOrd.coordAxisValues, 'frame'):
153                frame=gridOrd.coordAxisValues.frame
154    try:
155        del selection[timeAxis] #NOTE: Haven't resolved the single CRS issue so cannot subset by time within an individual file for now
156    except KeyError:
157        pass
158   
159    #get slice indices
160    selection=self.getSliceIndices(selection)                 #commented out while fixing subsetting
161
162    strTimes, axisorder, units, fulldata, fillvalue =csml.API.genSubset.getTheData(self, selection, times, timeName, timeAxis)
163    return pathToSubsetNetCDF, crs, frame, timeName, times, strTimes, caltype, axisorder,units, fulldata, fillvalue, kwargs
164   
165def subsetToGridSeries(self, outputdir=None, ncname='gridseries.nc',**kwargs):
166    '''Subset a GridSeriesFeature to another smaller  GridSeriesFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs
167    @param outputdir:    name of directory to write data files to
168    @param ncname:     name of netcdf file to write out
169    @param kwargs:      subset selection by axis name.
170    @return:     subsettedFeature (GridSeriesFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)
171   
172    '''
173    if outputdir is not None:
174        self.outputdir=outputdir
175    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
176        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
177    csml.csmllibs.csmlextra.checkDirExists(self.outputdir)
178    self.ncname=ncname
179    pathToSubsetNetCDF, crs, frame, timeName, times, strTimes, caltype, axisorder,units, fulldata, fillvalue, kwargs=self._subsetGrid(**kwargs) 
180    #Okay, got the data now. Need to write CSML feature and NetCDF files.
181    #Writing out the CSML feature   
182   
183    # define domain/coverage  to use in 'value' attribute   
184    newdomain=csml.parser.GridSeriesDomain()
185    domainSubset, totalArraySize=csml.API.genSubset.subsetDomain(timeName,strTimes,self.domain, **kwargs)
186    cTT=csml.API.genSubset.getCoordTransformTable(domainSubset, crs, frame)
187    newdomain.id=csml.csmllibs.csmlextra.getRandomID()
188    newdomain.coordTransformTable=cTT
189    newdomain.srsName=self.value.gridSeriesDomain.srsName 
190    newdomain.axisLabels=self.value.gridSeriesDomain.axisLabels
191    newdomain.srsDimension=self.value.gridSeriesDomain.srsDimension
192    newdomain.dimension=self.value.gridSeriesDomain.dimension
193    env=csml.parser.GridEnvelope()
194    env.low=csml.parser.csString('0 0 0') #TODO
195    env.high=csml.parser.csString('0 0 0')
196    newdomain.limits=env
197    newdomain.aLabels=self.value.gridSeriesDomain.aLabels
198    rs=csml.parser.RangeSet()
199    sdid=csml.csmllibs.csmlextra.getRandomID()
200    descriptor=csml.parser.NetCDFExtract(id=sdid,fileName=csml.parser.csString(self.ncname),variableName=self.name,arraySize=csml.parser.csString(totalArraySize))
201    rs=csml.parser.RangeSet()
202    va=csml.parser.ValueArray()
203    vc=csml.parser.MeasureOrNullList()
204    vc.href='#%s'%sdid
205    vc.arcrole="http://ndg.nerc.ac.uk/xlinkUsage/insert#QuantityList"
206    vc.role="http://ndg.nerc.ac.uk/fileFormat/csmlStorageDescriptor"
207    vc.show='embed'
208    try:
209        vc.uom=self.value.rangeSet.valueArray.valueComponent.uom
210    except:
211        vc.uom = 'unknown'  #TODO,
212    va.valueComponent=vc
213    va.id=csml.csmllibs.csmlextra.getRandomID()
214    rs.valueArray=va
215    #gridseries coverage
216    cvg=csml.parser.GridSeriesCoverage()
217    cvg.id=csml.csmllibs.csmlextra.getRandomID()
218    cvg.rangeSet=rs
219    cvg.gridSeriesDomain=newdomain   
220   
221    #parameter, as before subsetting.
222    param = self.parameter
223       
224    #create a stand alone gridseries feature containing this coverage
225    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
226    subsettedFeature=csmlWrap.createGridSeriesFeature(value=cvg,parameter=param,featureID=csml.csmllibs.csmlextra.getRandomID(),name=self.name,description=self.description)
227 
228    #Attempt to preserve any attributes from original data
229    preservedData=self._preserveAttributes()
230   
231 
232    ### write netcdf using NCWriter class (wraps cdms) ###
233    nc=csml.csmllibs.netCDFWriter.NCwriter(pathToSubsetNetCDF)
234    ords=cTT.gridOrdinates
235    axislist=[]
236    for a in axisorder:
237        axislist.append(self.gridnames[a])
238    stdname=param.getStandardName()
239   
240    #get preserved attributes for the variable and assoicated axes
241    preservedAttributes=preservedData.keyVariable.attributes
242    axisattdict=preservedData.keyVariable.axesattribs
243    #remove attributes that will clash with other values already defined
244    for kw in ['units', 'fillvalue']:
245        if kw in preservedAttributes:
246            del preservedAttributes[kw]
247        for item in axisattdict:
248            if kw in axisattdict[item]:
249                del axisattdict[item][kw]
250   
251    nc.genWriteVar(self.name.CONTENT,ords, times, caltype, axislist, units, stdname, fulldata, fillvalue,axisattribs=axisattdict, **preservedAttributes)
252   
253    #now add any other data we preserved earlier, e.g. projection variables:
254    for pv in preservedData.projectionVariables:
255        for pvaxis in pv.axes:
256            nc.addAxis(pvaxis,pv.axes[pvaxis], pv.attributes) #of form: nc.addAxis(axisname, data)
257        print pv.data
258        #TODO slice the extra variables
259        #nc.addVariable(pv.data, pv.name, pv.axes.keys(),**pv.attributes)   
260    nc.closeFinishedFile()
261    print 'NetCDF file written to %s'%pathToSubsetNetCDF
262    return subsettedFeature, pathToSubsetNetCDF, descriptor
263   
264
265def subsetToProfileSeries(self, outputdir=None, ncname='profileseries.nc',**kwargs):
266    '''Subset a GridSeriesFeature to a ProfileSeriesFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs
267    @param outputdir:    name of directory to write data files to
268    @param ncname:     name of netcdf file to write out
269    @param kwargs:      subset selection by axis name.
270    @return:     subsettedFeature (ProfileSeriesFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)'''
271   
272   
273    #TODO   !!!!!!!!! Need to perform some sort of testing on the kwargs to check it is a profileseries request.
274    if outputdir is not None:
275        self.outputdir=outputdir
276    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
277        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
278    csml.csmllibs.csmlextra.checkDirExists(self.outputdir)
279    self.ncname=ncname   
280   
281    #perform the subset (note this included nearest neighbour searching, so may return a different set of kwargs
282    pathToSubsetNetCDF, crs,frame, timeName, times, strTimes, caltype, axisorder,units, fulldata, fillvalue, kwargs=self._subsetGrid(**kwargs) 
283
284    latName=crs.axes[crs.latAxis]
285    lonName=crs.axes[crs.lonAxis]
286    #Okay, got the data now. Need to write CSML feature and NetCDF files.
287    #Writing out the CSML feature
288   
289    # define domain/coverage  to use in 'value' attribute   
290    newdomain=csml.parser.ProfileSeriesDomain()
291   
292   
293    domainSubset, totalArraySize=csml.API.genSubset.subsetDomain(timeName,strTimes,self.domain, **kwargs)
294    cTT=csml.API.genSubset.getCoordTransformTable(domainSubset, crs, frame) 
295   
296    #Need to remove location (lat/lon) axes from the cTT
297    newGridOrds=[]
298    for gridOrd in cTT.gridOrdinates:
299        if gridOrd.coordAxisLabel.CONTENT not in  [latName,lonName]:
300            newGridOrds.append(gridOrd)
301    cTT.gridOrdinates=newGridOrds
302    axes=[]
303    for ord in newGridOrds:
304        axes.append(ord.coordAxisLabel.CONTENT)
305    cat = csml.csmllibs.csmlcrs.CRSCatalogue()
306    crsInfo=cat.determineCRS(knownCRSAxes=axes)   
307    crs=cat.getCRS(crsInfo[0].srsName)
308   
309    newdomain.coordTransformTable=cTT   
310    newdomain.id=csml.csmllibs.csmlextra.getRandomID()
311    #TODO need to remove 'location' from this srs
312    newdomain.srsName=crsInfo[0].srsName
313    axisLabels=csml.csmllibs.csmlextra.stringify(crsInfo[1].keys())
314    newdomain.axisLabels=axisLabels
315    newdomain.srsDimension=2
316    newdomain.dimension=2
317    env=csml.parser.GridEnvelope()
318    env.low=csml.parser.csString('0 0 0') #TODO
319    env.high=csml.parser.csString('0 0 0')
320    newdomain.limits=env
321    newdomain.aLabels=self.value.gridSeriesDomain.aLabels #todo 
322    rangeSet=csml.parser.RangeSet()
323    descid=csml.csmllibs.csmlextra.getRandomID()
324    descriptor=csml.parser.NetCDFExtract(id=descid,fileName=csml.parser.csString(ncname),variableName=self.name,arraySize=csml.parser.csString(totalArraySize))
325    rs=csml.parser.RangeSet()
326    va=csml.parser.ValueArray()
327    vc=csml.parser.MeasureOrNullList()
328    vc.href='#%s'%descid
329    vc.arcrole="http://ndg.nerc.ac.uk/xlinkUsage/insert#QuantityList"
330    vc.role="http://ndg.nerc.ac.uk/fileFormat/csmlStorageDescriptor"
331    vc.show='embed'
332    vc.insertedExtract=descriptor
333    try:
334        vc.uom=self.value.rangeSet.valueArray.valueComponent.uom
335    except:
336        vc.uom = 'unknown'  #TODO,
337   
338    va.valueComponent=vc
339    va.id=csml.csmllibs.csmlextra.getRandomID()
340    rs.valueArray=va
341    #profileseries coverage
342    cvg=csml.parser.ProfileSeriesCoverage()
343    cvg.id=csml.csmllibs.csmlextra.getRandomID()
344    cvg.rangeSet=rs
345    cvg.profileSeriesDomain=newdomain   
346    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
347   
348    #parameter, as before subsetting.
349    param = self.parameter
350   
351    #create 'location' attribute:
352    if (latName, lonName) is not (None, None):
353        locStr ='%s %s'%(kwargs[latName], kwargs[lonName])
354        loc=csml.parser.csString(locStr)
355    else:
356        loc = csml.parser.csString('unknown location')
357       
358   
359    #create a stand alone gridseries feature containing this coverage
360    subsettedFeature=csmlWrap.createProfileSeriesFeature(value=cvg,parameter=param, location=loc, featureID=csml.csmllibs.csmlextra.getRandomID(),name=self.name,description=self.description)
361 
362   
363    ### write netcdf using NCWriter class (wraps cdms) ###
364    nc=csml.csmllibs.netCDFWriter.NCwriter(pathToSubsetNetCDF)
365    ords=cTT.gridOrdinates
366    axislist=[]
367       
368    for a in axisorder:
369        axislist.append(self.gridnames[a])
370    stdname=param.getStandardName()
371    nc.genWriteVar(self.name.CONTENT,ords, times, caltype, axislist, units, stdname, fulldata, fillvalue, latitude=kwargs[latName], longitude=kwargs[lonName])
372    nc.closeFinishedFile()
373   
374    return subsettedFeature, pathToSubsetNetCDF, descriptor
375
376def subsetToProfile(self, outputdir = None, ncname='profile.nc',**kwargs):
377    '''Subset a GridSeriesFeature to a ProfileFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs.
378       
379        Is a two step process - subsets GridSeries to ProfileSeries, then ProfileSeries to Profile
380    @param outputdir:    name of directory to write data files to
381    @param ncname:     name of netcdf file to write out
382    @param kwargs:      subset selection by axis name.
383    @return:     subsettedFeature (ProfileFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)'''
384
385    odir=outputdir
386    tmpfile= csml.csmllibs.csmlextra.getRandomID()+'.nc'
387    profileSeries, pSfile, descriptor=self.subsetToProfileSeries(odir, tmpfile, **kwargs)   
388    del kwargs['latitude']    #TODO - need to remove excess kwargs based on the domain of the temporary profileSeries feature
389    del kwargs['longitude']
390    subsettedFeature, pathToSubsetNetCDF, descriptor =profileSeries.subsetToProfile(odir,ncname, **kwargs)   
391    tmpfile=odir + '/' + tmpfile
392    os.remove(tmpfile)
393    return subsettedFeature, pathToSubsetNetCDF, descriptor
394
395def subsetToPointSeries(self, outputdir =None, ncname='pointseries.nc',**kwargs):
396    '''Subset a GridSeriesFeature to a PointSeriesFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs.
397       
398        Is a two step process - subsets GridSeries to ProfileSeries, then ProfileSeries to PointSeries
399    @param outputdir:    name of directory to write data files to
400    @param ncname:     name of netcdf file to write out
401    @param kwargs:      subset selection by axis name.
402    @return:     subsettedFeature (PointSeriesFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)'''
403    tmpfile= csml.csmllibs.csmlextra.getRandomID()+'.nc'
404    profileSeries, pSfile, descriptor =self.subsetToProfileSeries(outputdir, tmpfile, **kwargs)   
405    del kwargs['latitude']    #TODO - need to remove excess kwargs based on the domain of the temporary profileSeries feature
406    del kwargs['longitude']
407    subsettedFeature, pathToSubsetNetCDF, descriptor=profileSeries.subsetToPointSeries(outputdir, ncname, **kwargs)   
408    tmpfile=outputdir + '/' + tmpfile
409    os.remove(tmpfile)
410    return subsettedFeature, pathToSubsetNetCDF, descriptor
Note: See TracBrowser for help on using the repository browser.