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

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

changed time to t

Line 
1''' ops_GridSeriesFeature  contains operations for GridSeriesFeatures'''
2from API import *
3from CSMLDocument import *
4from Scientific.IO.NetCDF import *  #use this instead of cdms for now for it's simple write interface..
5from Numeric import *
6
7
8def testmethod(self):
9    #print 'testmethod for gridseries feature'
10    return 'testmethod - gridseries'
11
12def getAllowedSubsettings(self):
13    return ['subsetToGridSeries']
14
15def getDomainReference(self):
16    #Inspects a time position list for the domain reference.
17    #TODO, does not handle a file extract in place of a list.
18    if isinstance(self.domain.domainReference,Parser.TimePositionList):
19        time = {}
20        time['t'] = self.domain.domainReference.timePositions
21        domainref  = time
22    return domainref
23
24def getDomainComplement(self):
25    #This will return a list containing one or more ordinates:
26    #currently in form [Name, values]
27    #assumes ordinate.axisValues is a file extract
28    #TODO axisValues may be inline -see also domainReference for similar problem
29    domaincomp ={}
30    dc = self.domain.domainComplement
31    #dc should be a grid!
32    if isinstance(dc, Parser.Grid):
33        for ordinate in dc.ordinates:
34            domaincomp[ordinate.definesAxis]=ordinate.axisValues.getData()
35    return domaincomp
36
37def getDomain(self):
38    #returns both the domain reference axes and domain compliment axes in a single domain dictionary
39    #axes are in no particular order
40    domain = {}
41    dr=getDomainReference(self)
42    dc=getDomainComplement(self)
43    for key in dc.keys():
44        domain[key]=dc[key]
45    for key in dr.keys():
46        domain[key]=dr[key]
47    return domain
48
49
50
51
52def subsetToGridSeries(self, timeSubset,  csmlpath=None, ncpath=None,**kwargs):
53    if csmlpath is not None:
54        pathToSubsetCSML = csmlpath
55    else:
56        pathToSubsetCSML='temp.xml'
57    if ncpath is not None:
58        pathToSubsetNetCDF=ncpath
59    else:
60        pathToSubsetNetCDF='temp.nc'
61
62    domainref = getDomainReference(self)
63    self.times=timeSubset
64    self.files=[]
65    strTimes=''
66    fulldata=[]
67    if len(self.times) == 2:
68        tone=ops_AbstractFeature.__getCDtime(self.times[0])
69        ttwo=ops_AbstractFeature.__getCDtime(self.times[1])
70        dr=getDomainReference(self)
71        self.times=[]
72        for time in dr['t'].split():
73            timeok=ops_AbstractFeature.__compareTimes(tone,time,ttwo)
74            if timeok ==1:
75                self.times.append(time)
76    for time in self.times:
77        listPosition=domainref['t'].split().index(time)
78        strTimes= strTimes + ' ' + time
79        #print self.rangeSet.aggregatedArray.component
80        timedim=self.rangeSet.aggregatedArray.component[0].variableName
81       
82        for comp in self.rangeSet.aggregatedArray.component:
83            data=comp.getData(fileposition=listPosition, **kwargs)
84            self.files.append(comp.fileName.split()[listPosition])
85            #sys.exit()
86
87            if fulldata ==[]:
88                fulldata = data.tolist()
89         #       print shape(fulldata)
90            #    sys.exit()
91            else:
92                for item in data.tolist():
93                    fulldata.append(item)
94        axisorder = data.getAxisIds()  #will need later!
95    #get the calendar type
96    try:
97        timedim ='t'   #TODO should accept any time dim!!
98        caltype, calunits = ops_AbstractFileExtract.__calendar(self.rangeSet.aggregatedArray.component[0].fileName.split()[0], timedim)   
99        csmltime.setcdtimeCalendar(caltype)
100    except:
101        csmltime.setcdtimeCalendar(csmltime.cdtime.DefaultCalendar)
102    ### define domain and rangeSet to use for feature in csml document####
103    domain=Parser.GridSeriesDomain()
104    domain.domainReference=Parser.TimePositionList(timePositions=strTimes)
105    grid=Parser.Grid()
106    dc = self.getDomainComplement()
107    ordinates= []
108    i=0
109    valueStore=[]  # use the values again later to generate netcdf
110    arraySize=0
111    totalArraySize=1
112    for key in dc.keys():
113        arraySize=0
114        i=i+1
115        god=Parser.GridOrdinateDescription()
116        god.gridAxesSpanned='dim%s'%i
117        god.sequenceRule='+x+y+z'
118        god.definesAxis=key
119        straxisValues=''
120        if key in kwargs:
121            for val in dc[key]:
122                if val >= kwargs[key][0]:
123                    if val <= kwargs[key] [1]:
124                        arraySize=arraySize+1
125                        straxisValues=straxisValues+ str(val) + ', '
126        else: # this dimension has not been subsetted
127            for val in dc[key]:
128                arraySize=arraySize+1
129                straxisValues=straxisValues+ str(val) + ', '
130        totalArraySize=totalArraySize*arraySize
131        god.axisValues=straxisValues[:-2]
132        ordinates.append(god)
133    totalArraySize=totalArraySize*len(self.times)
134    grid.ordinates=ordinates
135    domain.domainComplement=grid
136    rangeSet=Parser.RangeSet()
137    rangeSet.arrayDescriptor=Parser.NetCDFExtract(id=self.id,fileName=pathToSubsetNetCDF,variableName=self.id,arraySize=[arraySize])
138
139    #### write csml document #####
140    csmldoc=CSMLDocument("mydoc", "mymetadata")
141    csmldoc.addGridSeriesFeature(domain,rangeSet,datasetID="A",featureID="B",description="C")
142    csmldoc=csmldoc.consolidate()
143    output=open(pathToSubsetCSML,'w')
144    output.write(csmldoc)
145    output.close()
146
147    ### create and write netcdf - uses scientific python####
148    ncfile=NetCDFFile(pathToSubsetNetCDF,'w')
149    # create the dimensions
150    ncfile.createDimension ( 't', len(self.times))
151    time_var = ncfile.createVariable ( 't', Float, ('t',) )
152    time_var.longname = 't'
153    floatTimes=[]
154    #print len(fulldata[0])
155    for time in self.times:
156        time=ops_AbstractFeature.__getCDtime(time).torel(calunits)
157        floatTimes.append(time.value)
158    time_var[:] =floatTimes[:]
159    time_var.units=calunits
160    time_var.calendar=caltype
161
162    for ordinate in ordinates:
163        ncfile.createDimension(ordinate.definesAxis, len(ordinate.axisValues.split()))
164        item_var = ncfile.createVariable (ordinate.definesAxis, Float, (ordinate.definesAxis,) )
165        #convert to list
166        vals=[]
167        for val in ordinate.axisValues.split(','):
168            vals.append(float(val))
169        #print ordinate.definesAxis
170        if ordinate.definesAxis=='latitude':
171            #this is just a fix till I work out what's going wrong
172            vals.reverse()
173        #print vals
174        item_var[:]=vals[:]
175    #this needs reconsidering - do the shapes always match up??
176
177
178    if len(ordinates)==3:
179        feature_var = ncfile.createVariable (self.id, Float,('t',axisorder[1],axisorder[2],axisorder[3]))
180    elif len(grid.ordinates)==2:
181        feature_var =ncfile.createVariable (self.id, Float,('t',axisorder[1],axisorder[2]))
182    feature_var[:]=fulldata[:]
183    ncfile.close()
184    return pathToSubsetCSML,pathToSubsetNetCDF, totalArraySize
Note: See TracBrowser for help on using the repository browser.