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

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

exposing getUnits methods etc

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