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

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

fixed bug that was caused by sorting axes list unintentionally

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