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

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

fixing ID/naming issues and paths

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