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

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

Global attribtues now preserved except when they are data dependent CF attributes (which need recalculating). So far this is only implemented for gridseries to gridseries

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    #TODO: this should be a separate method
255    sliceinfo=csml.csmllibs.csmlpreservation.getDomainSlice(domainSubset, self.domain)
256    for pv in preservedData.projectionVariables:
257        pvslice=[]
258        for pvaxis in pv.axes:
259            nc.addAxis(pvaxis,pv.axes[pvaxis][sliceinfo[pvaxis][0]:sliceinfo[pvaxis][1]], pv.attributes) #of form: nc.addAxis(axisname, data)
260            pvslice.append(tuple([sliceinfo[pvaxis][0],sliceinfo[pvaxis][1]]))
261
262        #add the sliced (or not) variables.
263        if len(pvslice) > 0:
264            nc.addVariable(pv.data(*pvslice), pv.name, pv.axes.keys(),**pv.attributes)
265        else:
266            nc.addVariable(pv.data, pv.name, pv.axes.keys(),**pv.attributes)
267   
268    if preservedData.globalAttributes !={}:
269        nc.addGlobalAttributes(preservedData.globalAttributes)
270    nc.closeFinishedFile()
271    print 'NetCDF file written to %s'%pathToSubsetNetCDF
272    return subsettedFeature, pathToSubsetNetCDF, descriptor
273   
274
275def subsetToProfileSeries(self, outputdir=None, ncname='profileseries.nc',**kwargs):
276    '''Subset a GridSeriesFeature to a ProfileSeriesFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs
277    @param outputdir:    name of directory to write data files to
278    @param ncname:     name of netcdf file to write out
279    @param kwargs:      subset selection by axis name.
280    @return:     subsettedFeature (ProfileSeriesFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)'''
281   
282   
283    #TODO   !!!!!!!!! Need to perform some sort of testing on the kwargs to check it is a profileseries request.
284    if outputdir is not None:
285        self.outputdir=outputdir
286    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
287        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
288    csml.csmllibs.csmlextra.checkDirExists(self.outputdir)
289    self.ncname=ncname   
290   
291    #perform the subset (note this included nearest neighbour searching, so may return a different set of kwargs
292    pathToSubsetNetCDF, crs,frame, timeName, times, strTimes, caltype, axisorder,units, fulldata, fillvalue, kwargs=self._subsetGrid(**kwargs) 
293
294    latName=crs.axes[crs.latAxis]
295    lonName=crs.axes[crs.lonAxis]
296    #Okay, got the data now. Need to write CSML feature and NetCDF files.
297    #Writing out the CSML feature
298   
299    # define domain/coverage  to use in 'value' attribute   
300    newdomain=csml.parser.ProfileSeriesDomain()
301   
302   
303    domainSubset, totalArraySize=csml.API.genSubset.subsetDomain(timeName,strTimes,self.domain, **kwargs)
304    cTT=csml.API.genSubset.getCoordTransformTable(domainSubset, crs, frame) 
305   
306    #Need to remove location (lat/lon) axes from the cTT
307    newGridOrds=[]
308    for gridOrd in cTT.gridOrdinates:
309        if gridOrd.coordAxisLabel.CONTENT not in  [latName,lonName]:
310            newGridOrds.append(gridOrd)
311    cTT.gridOrdinates=newGridOrds
312    axes=[]
313    for ord in newGridOrds:
314        axes.append(ord.coordAxisLabel.CONTENT)
315    cat = csml.csmllibs.csmlcrs.CRSCatalogue()
316    crsInfo=cat.determineCRS(knownCRSAxes=axes)   
317    crs=cat.getCRS(crsInfo[0].srsName)
318   
319    newdomain.coordTransformTable=cTT   
320    newdomain.id=csml.csmllibs.csmlextra.getRandomID()
321    #TODO need to remove 'location' from this srs
322    newdomain.srsName=crsInfo[0].srsName
323    axisLabels=csml.csmllibs.csmlextra.stringify(crsInfo[1].keys())
324    newdomain.axisLabels=axisLabels
325    newdomain.srsDimension=2
326    newdomain.dimension=2
327    env=csml.parser.GridEnvelope()
328    env.low=csml.parser.csString('0 0 0') #TODO
329    env.high=csml.parser.csString('0 0 0')
330    newdomain.limits=env
331    newdomain.aLabels=self.value.gridSeriesDomain.aLabels #todo 
332    rangeSet=csml.parser.RangeSet()
333    descid=csml.csmllibs.csmlextra.getRandomID()
334    descriptor=csml.parser.NetCDFExtract(id=descid,fileName=csml.parser.csString(ncname),variableName=self.name,arraySize=csml.parser.csString(totalArraySize))
335    rs=csml.parser.RangeSet()
336    va=csml.parser.ValueArray()
337    vc=csml.parser.MeasureOrNullList()
338    vc.href='#%s'%descid
339    vc.arcrole="http://ndg.nerc.ac.uk/xlinkUsage/insert#QuantityList"
340    vc.role="http://ndg.nerc.ac.uk/fileFormat/csmlStorageDescriptor"
341    vc.show='embed'
342    vc.insertedExtract=descriptor
343    try:
344        vc.uom=self.value.rangeSet.valueArray.valueComponent.uom
345    except:
346        vc.uom = 'unknown'  #TODO,
347   
348    va.valueComponent=vc
349    va.id=csml.csmllibs.csmlextra.getRandomID()
350    rs.valueArray=va
351    #profileseries coverage
352    cvg=csml.parser.ProfileSeriesCoverage()
353    cvg.id=csml.csmllibs.csmlextra.getRandomID()
354    cvg.rangeSet=rs
355    cvg.profileSeriesDomain=newdomain   
356    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
357   
358    #parameter, as before subsetting.
359    param = self.parameter
360   
361    #create 'location' attribute:
362    if (latName, lonName) is not (None, None):
363        locStr ='%s %s'%(kwargs[latName], kwargs[lonName])
364        loc=csml.parser.csString(locStr)
365    else:
366        loc = csml.parser.csString('unknown location')
367       
368   
369    #create a stand alone gridseries feature containing this coverage
370    subsettedFeature=csmlWrap.createProfileSeriesFeature(value=cvg,parameter=param, location=loc, featureID=csml.csmllibs.csmlextra.getRandomID(),name=self.name,description=self.description)
371 
372   
373    ### write netcdf using NCWriter class (wraps cdms) ###
374    nc=csml.csmllibs.netCDFWriter.NCwriter(pathToSubsetNetCDF)
375    ords=cTT.gridOrdinates
376    axislist=[]
377       
378    for a in axisorder:
379        axislist.append(self.gridnames[a])
380    stdname=param.getStandardName()
381    nc.genWriteVar(self.name.CONTENT,ords, times, caltype, axislist, units, stdname, fulldata, fillvalue, latitude=kwargs[latName], longitude=kwargs[lonName])
382    nc.closeFinishedFile()
383   
384    return subsettedFeature, pathToSubsetNetCDF, descriptor
385
386def subsetToProfile(self, outputdir = None, ncname='profile.nc',**kwargs):
387    '''Subset a GridSeriesFeature to a ProfileFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs.
388       
389        Is a two step process - subsets GridSeries to ProfileSeries, then ProfileSeries to Profile
390    @param outputdir:    name of directory to write data files to
391    @param ncname:     name of netcdf file to write out
392    @param kwargs:      subset selection by axis name.
393    @return:     subsettedFeature (ProfileFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)'''
394
395    odir=outputdir
396    tmpfile= csml.csmllibs.csmlextra.getRandomID()+'.nc'
397    profileSeries, pSfile, descriptor=self.subsetToProfileSeries(odir, tmpfile, **kwargs)   
398    del kwargs['latitude']    #TODO - need to remove excess kwargs based on the domain of the temporary profileSeries feature
399    del kwargs['longitude']
400    subsettedFeature, pathToSubsetNetCDF, descriptor =profileSeries.subsetToProfile(odir,ncname, **kwargs)   
401    tmpfile=odir + '/' + tmpfile
402    os.remove(tmpfile)
403    return subsettedFeature, pathToSubsetNetCDF, descriptor
404
405def subsetToPointSeries(self, outputdir =None, ncname='pointseries.nc',**kwargs):
406    '''Subset a GridSeriesFeature to a PointSeriesFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs.
407       
408        Is a two step process - subsets GridSeries to ProfileSeries, then ProfileSeries to PointSeries
409    @param outputdir:    name of directory to write data files to
410    @param ncname:     name of netcdf file to write out
411    @param kwargs:      subset selection by axis name.
412    @return:     subsettedFeature (PointSeriesFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)'''
413    tmpfile= csml.csmllibs.csmlextra.getRandomID()+'.nc'
414    profileSeries, pSfile, descriptor =self.subsetToProfileSeries(outputdir, tmpfile, **kwargs)   
415    del kwargs['latitude']    #TODO - need to remove excess kwargs based on the domain of the temporary profileSeries feature
416    del kwargs['longitude']
417    subsettedFeature, pathToSubsetNetCDF, descriptor=profileSeries.subsetToPointSeries(outputdir, ncname, **kwargs)   
418    tmpfile=outputdir + '/' + tmpfile
419    os.remove(tmpfile)
420    return subsettedFeature, pathToSubsetNetCDF, descriptor
Note: See TracBrowser for help on using the repository browser.