source: TI02-CSML/trunk/parser/API/ops_GridSeriesFeature.py @ 1196

''' ops_GridSeriesFeature  contains operations for GridSeriesFeatures'''
from API import *
from CSMLDocument import *
from Scientific.IO.NetCDF import *  #use this instead of cdms for now for it's simple write interface..
from Numeric import *


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

def getAllowedSubsettings(self):
    return ['subsetToGridSeries']

def getDomainReference(self):
    #Inspects a time position list for the domain reference.
    #TODO, does not handle a file extract in place of a list.
    if isinstance(self.domain.domainReference,Parser.TimePositionList):
        time = {}
        time['t'] = self.domain.domainReference.timePositions
        domainref  = time
    return domainref

def getDomainComplement(self):
    #This will return a list containing one or more ordinates:
    #currently in form [Name, values]
    #assumes ordinate.axisValues is a file extract
    #TODO axisValues may be inline -see also domainReference for similar problem
    domaincomp ={}
    dc = self.domain.domainComplement
    #dc should be a grid!
    if isinstance(dc, Parser.Grid):
        for ordinate in dc.ordinates:
            domaincomp[ordinate.definesAxis]=ordinate.axisValues.getData()
    return domaincomp

def getDomain(self):
    #returns both the domain reference axes and domain compliment axes in a single domain dictionary
    #axes are in no particular order
    domain = {}
    dr=getDomainReference(self)
    dc=getDomainComplement(self)
    for key in dc.keys():
        domain[key]=dc[key]
    for key in dr.keys():
        domain[key]=dr[key]
    return domain




def subsetToGridSeries(self, timeSubset,  csmlpath=None, ncpath=None,**kwargs):
    if csmlpath is not None:
        pathToSubsetCSML = csmlpath
    else:
        pathToSubsetCSML='temp.xml'
    if ncpath is not None:
        pathToSubsetNetCDF=ncpath
    else:
        pathToSubsetNetCDF='temp.nc'

    domainref = getDomainReference(self)
    self.times=timeSubset
    self.files=[]
    strTimes=''
    fulldata=[]
    if len(self.times) == 2:
        tone=ops_AbstractFeature.__getCDtime(self.times[0])
        ttwo=ops_AbstractFeature.__getCDtime(self.times[1])
        dr=getDomainReference(self)
        self.times=[]
        for time in dr['t'].split():
            timeok=ops_AbstractFeature.__compareTimes(tone,time,ttwo)
            if timeok ==1:
                self.times.append(time)
    for time in self.times:
        listPosition=domainref['t'].split().index(time)
        strTimes= strTimes + ' ' + time
        #print self.rangeSet.aggregatedArray.component
        for comp in self.rangeSet.aggregatedArray.component:
            data=comp.getData(fileposition=listPosition, **kwargs)
            self.files.append(comp.fileName.split()[listPosition])
            #sys.exit()

            if fulldata ==[]:
                fulldata = data.tolist()
         #       print shape(fulldata)
            #    sys.exit()
            else:
                for item in data.tolist():
                    fulldata.append(item)
        axisorder = data.getAxisIds()  #will need later!
    #get the calendar type
    try:
        caltype, calunits = ops_AbstractFileExtract.__calendar(self.rangeSet.aggregatedArray.component[0].fileName.split()[0], 't')    #TODO should accept any time dim!!
        csmltime.setcdtimeCalendar(caltype)
    except:
        csmltime.setcdtimeCalendar(csmltime.cdtime.DefaultCalendar)
    ### define domain and rangeSet to use for feature in csml document####
    domain=Parser.GridSeriesDomain()
    domain.domainReference=Parser.TimePositionList(timePositions=strTimes)
    grid=Parser.Grid()
    dc = self.getDomainComplement()
    ordinates= []
    i=0
    valueStore=[]  # use the values again later to generate netcdf
    arraySize=0
    totalArraySize=1
    for key in dc.keys():
        arraySize=0
        i=i+1
        god=Parser.GridOrdinateDescription()
        god.gridAxesSpanned='dim%s'%i
        god.sequenceRule='+x+y+z'
        god.definesAxis=key
        straxisValues=''
        if key in kwargs:
            for val in dc[key]:
                if val >= kwargs[key][0]:
                    if val <= kwargs[key] [1]:
                        arraySize=arraySize+1
                        straxisValues=straxisValues+ str(val) + ', '
        else: # this dimension has not been subsetted
            for val in dc[key]:
                arraySize=arraySize+1
                straxisValues=straxisValues+ str(val) + ', '
        totalArraySize=totalArraySize*arraySize
        god.axisValues=straxisValues[:-2]
        ordinates.append(god)
    totalArraySize=totalArraySize*len(self.times)
    grid.ordinates=ordinates
    domain.domainComplement=grid
    rangeSet=Parser.RangeSet()
    rangeSet.arrayDescriptor=Parser.NetCDFExtract(id=self.id,fileName=pathToSubsetNetCDF,variableName=self.id,arraySize=[arraySize])

    #### write csml document ##### -move this to the csmldocument module?
    subsetCSML=CSMLDocument()
    subsetCSML=subsetCSML.makeGridSeries(domain,rangeSet)
    output=open(pathToSubsetCSML,'w')
    output.write(subsetCSML)
    output.close()

    ### create and write netcdf - uses scientific python####
    ncfile=NetCDFFile(pathToSubsetNetCDF,'w')
    # create the dimensions
    ncfile.createDimension ( 'time', len(self.times))
    time_var = ncfile.createVariable ( 'time', Float, ('time',) )
    time_var.longname = 'time'
    floatTimes=[]
    #print len(fulldata[0])
    for time in self.times:
        time=ops_AbstractFeature.__getCDtime(time).torel(calunits)
        floatTimes.append(time.value)
    time_var[:] =floatTimes[:]
    time_var.units=calunits
    time_var.calendar=caltype

    for ordinate in ordinates:
        ncfile.createDimension(ordinate.definesAxis, len(ordinate.axisValues.split()))
        item_var = ncfile.createVariable (ordinate.definesAxis, Float, (ordinate.definesAxis,) )
        #convert to list
        vals=[]
        for val in ordinate.axisValues.split(','):
            vals.append(float(val))
        #print ordinate.definesAxis
        if ordinate.definesAxis=='latitude':
            #this is just a fix till I work out what's going wrong
            vals.reverse()
        #print vals
        item_var[:]=vals[:]
    #this needs reconsidering - do the shapes always match up??



    if len(ordinates)==3:
        feature_var = ncfile.createVariable (self.id,
Float,('time',axisorder[1],axisorder[2],axisorder[3]))
    elif len(grid.ordinates)==2:
        feature_var =ncfile.createVariable (self.id,
Float,('time',axisorder[1],axisorder[2]))
    feature_var[:]=fulldata[:]
    ncfile.close()
    return pathToSubsetCSML,pathToSubsetNetCDF, totalArraySize