source: TI02-CSML/branches/csml-cdms2/API/ops_PointSeriesFeature.py @ 3627

# Adapted for numpy/ma/cdms2 by convertcdms.py
''' ops_PointSeriesFeature  contains operations for PointSeriesFeatures. These methods are attached to PointSeriesFeature instances at runtime.'''
import csml,sys
import cdms2 as cdms, MV2 as MV, numpy.oldnumeric.ma as MA

def testmethod(self):
    ''' testmethod for gridseries feature'''
    print 'testmethod for pointSeries feature'
    return 'testmethod pointSeries'

def getAllowedSubsettings(self):
    '''get the allowed operations
    @return:    list of operations'''
    return ['subsetToPointSeries']
    
def getDomain(self):
    ''' returns domain as a dictionary of ordinates {name: [values], ...} '''
    self.domain={}
    valList=[]
    for val in self.value.pointSeriesDomain.timePositionList.CONTENT.split():
        valList.append(val)
    self.domain['time']=valList
    [lat,lon]=self.location.CONTENT.split()
    self.domain['longitude']=lon
    self.domain['latitude']=lat

    return self.domain    
   
def getLongitudeAxis(self):
    ''' get the index of the longitude axis'''
    return 'longitude' #dummy value not used, subsetting is by time only

def getLatitudeAxis(self):
    ''' get the index of the latitude axis'''
    return 'latitude' #dummy value not used,  subsetting is by time only

def getTimeAxis(self):
    ''' get the index of the time axis'''
    return 'times'
    
def subsetToPointSeries(self, outputdir=None, ncname='pointseries.nc' ,times=None):
    '''Subset a PointSeriesFeature to a PointSeriesFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs.
    
    @param outputdir:    name of directory to write data files to
    @param ncname:     name of netcdf file to write out
    @param times:       selection by times
    @return:     subsettedFeature (PointSeriesFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)'''
    if outputdir is not None:
        self.outputdir=outputdir
    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
      
    if self.outputdir is not None:
        pathToSubsetNetCDF=self.outputdir + '/' +ncname
    else: 
        pathToSubsetNetCDF=ncname
    if self.outputdir is not None:
        csml.csmllibs.csmlextra.checkDirExists(self.outputdir)
  
    if type(times) is str:
        selection = tuple(csml.csmllibs.csmlextra.listify(times))
    else:
        selection = times  
    tlist=[]    
    for val in self.value.pointSeriesDomain.timePositionList.CONTENT.split():
        tlist.append(val)
   
    indices=[]
    if  len(selection) ==2:       
        selectedTimes=[]
        for domaintime in tlist:
            if selection[0] == domaintime:
                minIndex=tlist.index(domaintime)
                selectedTimes.append(domaintime)
            if domaintime > selection[0]:
                if domaintime < selection[1]:
                    selectedTimes.append(domaintime)          
            if selection[1] == domaintime:
                maxIndex=tlist.index(domaintime)
                selectedTimes.append(domaintime)          
                break   
                    
    elif len(selection)==1:   #point only
        for domaintime in tlist:
            if selection[0] == domaintime:
                
                minIndex=tlist.index(domaintime)
                maxIndex=minIndex
                selectedTimes=selection
                break
    
    
    
    #fileList=[]
        
    if hasattr(self.value.rangeSet, 'quantityList'):        
        #inline
        data=self.value.rangeSet.quantityList.CONTENT.split()[minIndex:maxIndex+1]
        fillvalue=None
        units=self.value.rangeSet.quantityList.uom
    else:
        #file extract
        if hasattr(self.value.rangeSet, 'valueArray'):
            fextract=self.value.rangeSet.valueArray.valueComponent.insertedExtract.components
            uom = self.value.rangeSet.valueArray.valueComponent.insertedExtract.uom
        else:
            fextract=self.value.rangeSet.arrayDescriptor.components
            uom = self.value.rangeSet.arrayDescriptor.uom
        #for f in fextract.fileList.fileNames.CONTENT.split():
            #fileList.append(f)
        
        data, fillvalue=fextract.getDataFromChunks(minIndex, maxIndex)
        units=[uom.CONTENT]
   
    #Now write out the CSML Feature:
    #create a PointSeriesCoverage
    cvg=csml.parser.PointSeriesCoverage()
    cvg.id=csml.csmllibs.csmlextra.getRandomID()
    
    #create a TimeSeriesDomain (appropriate domain for a PointSeriesFeature
    ptsd=csml.parser.TimeSeries()
    ptsd.id=csml.csmllibs.csmlextra.getRandomID()
    #create (and populate) a TimePositionList. Using keyword arguements for conciseness.
    ptsd.timePositionList=csml.parser.TimePositionList(CONTENT=csml.csmllibs.csmlextra.stringify(selectedTimes))
        
    sdid=csml.csmllibs.csmlextra.getRandomID() # for the storage descriptor
    #create a RangeSet
    rs=csml.parser.RangeSet()
    va=csml.parser.ValueArray()
    vc=csml.parser.MeasureOrNullList()
    vc.href='#%s'%sdid
    vc.arcrole="http://ndg.nerc.ac.uk/xlinkUsage/insert#QuantityList"
    vc.role="http://ndg.nerc.ac.uk/fileFormat/csmlStorageDescriptor"
    vc.show='embed'
    try:
        vc.uom=self.value.rangeSet.valueArray.valueComponent.uom
    except:
        vc.uom = 'unknown'  #TODO,
    va.valueComponent=vc
    va.id=csml.csmllibs.csmlextra.getRandomID()
    rs.valueArray=va
    
    
    #Add the domain and rangeSet as attributes of the coverage
    cvg.pointSeriesDomain=ptsd
    cvg.rangeSet=rs
    
    totalArraySize=10
    if hasattr(self, 'name'):
        varname=self.name.CONTENT
    else:
        varname=self.id
    descriptor=csml.parser.NetCDFExtract(id=sdid,fileName=csml.parser.csString(ncname),variableName=csml.parser.csString(varname),arraySize=csml.parser.csString(totalArraySize))

    #the parameter of the feature is of type Phenomenon, here href creates "xlink:href=..."
    param=self.parameter
    
      
    #create a stand alone pointseries feature containing this coverage
    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
    subsettedFeature=csmlWrap.createPointSeriesFeature(value=cvg,parameter=param,location=self.location,featureID=self.id,description=self.description)
      
    
    #Now write out the NetCDF File:   
    nc=cdms.open(pathToSubsetNetCDF,'w')
    var=MV.array(data)
    var.id=varname
    var.name=varname
    var.units=units[0] # hopefully there is just one unit in this list..
    setattr(var, 'missing_value' ,fillvalue)
    
    floatTimes=[]
    
    tOne=csml.csmllibs.csmltime.getCDtime(selection[0])
    tbase=csml.csmllibs.csmltime.getBaseUnits(tOne)
    for time in selectedTimes:
        time=csml.csmllibs.csmltime.getCDtime(time).torel(tbase)
        floatTimes.append(time.value)
    timeAx=cdms.createAxis(floatTimes)
    timeAx.designateTime()
    timeAx.id='time'
    timeAx.units=tbase  
    var.setAxis(0,timeAx)      
    nc.write(var)
    nc.close()   
    print 'NetCDF file written to %s'%pathToSubsetNetCDF
    return subsettedFeature, pathToSubsetNetCDF, descriptor
    
    
def subsetToPoint(self, outputdir=None, ncname='point.nc',time=None):
    '''Subset a PointSeriesFeature to a PointFeature. Performs the subset (note this includes nearest neighbour searching, so may return a different set of kwargs.
        
    @param outputdir:    name of directory to write data files to
    @param ncname:     name of netcdf file to write out
    @param time:       selection by single time
    @return:     subsettedFeature (PointFeature instance)  pathToSubsetNetCDF (filepath), descriptor (array descriptor instance)'''
    if outputdir is not None:
        self.outputdir=outputdir
    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR 
    pseriesFeature, pathToSubsetNetCDF, descriptor=self.subsetToPointSeries(self.outputdir, ncname,time)          
    #now need to take this point series feature containing one time and rewrite it as a point feature.
    pointFeature=csml.parser.PointFeature()
    pointFeature.location=pseriesFeature.location
    pointFeature.parameter=pseriesFeature.parameter
    pointFeature.description=pseriesFeature.description
    cvg=csml.parser.PointCoverage()
    domain=csml.parser.PointDomain()
    domain.id=csml.csmllibs.csmlextra.getRandomID()
    domain.pointMember=csml.parser.csString(pseriesFeature.location.CONTENT)
    #create (and populate) a TimePositionList. Using keyword arguements for conciseness.
    cvg.pointDomain=domain
    pointFeature.value=cvg
    pointFeature.time=csml.parser.csString(time)
    return pointFeature, pathToSubsetNetCDF, descriptor