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

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

added code to preserve CF projection variables and subset them accordingly

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