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

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

fixed subsetting bug when handling times, and also removed dependence on longitude/latitude names

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