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

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

subsetting by index initially working for rawfileextracts (more testing needed)

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