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

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

fixed bug that was caused by sorting axes list unintentionally #2

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