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

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

fixed problem with axis ordering in csml api

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        pathToSubsetCSML='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        for comp in self.rangeSet.aggregatedArray.component:
81            data=comp.getData(fileposition=listPosition, **kwargs)
82            self.files.append(comp.fileName.split()[listPosition])
83            #sys.exit()
84           
85            if fulldata ==[]:
86                fulldata = data.tolist()
87                print shape(fulldata)
88            #    sys.exit()
89            else:
90                for item in data.tolist():
91                    fulldata.append(item)
92        axisorder = data.getAxisIds()  #will need later!
93    #get the calendar type
94    try:
95        caltype, calunits = ops_AbstractFileExtract.__calendar(self.rangeSet.aggregatedArray.component[0].fileName.split()[0], 't')    #TODO should accept any time dim!!
96        csmltime.setcdtimeCalendar(caltype)
97    except:
98        csmltime.setcdtimeCalendar(csmltime.cdtime.DefaultCalendar)
99    ### define domain and rangeSet to use for feature in csml document####
100    domain=Parser.GridSeriesDomain()
101    domain.domainReference=Parser.TimePositionList(timePositions=strTimes) 
102    grid=Parser.Grid()
103    dc = self.getDomainComplement()
104    ordinates= []
105    i=0
106    valueStore=[]  # use the values again later to generate netcdf
107    arraySize=0
108    totalArraySize=1
109    for key in dc.keys():
110        arraySize=0
111        i=i+1
112        god=Parser.GridOrdinateDescription()
113        god.gridAxesSpanned='dim%s'%i
114        god.sequenceRule='+x+y+z'
115        god.definesAxis=key
116        straxisValues=''
117        if key in kwargs:
118            for val in dc[key]:
119                if val >= kwargs[key][0]:
120                    if val <= kwargs[key] [1]:
121                        arraySize=arraySize+1
122                        straxisValues=straxisValues+ str(val) + ', '
123        else: # this dimension has not been subsetted
124            for val in dc[key]:
125                arraySize=arraySize+1
126                straxisValues=straxisValues+ str(val) + ', '
127        totalArraySize=totalArraySize*arraySize
128        god.axisValues=straxisValues[:-2]
129        ordinates.append(god)
130    totalArraySize=totalArraySize*len(self.times)
131    grid.ordinates=ordinates
132    domain.domainComplement=grid
133    rangeSet=Parser.RangeSet()
134    rangeSet.arrayDescriptor=Parser.NetCDFExtract(id=self.id,fileName=pathToSubsetCSML,variableName=self.id,arraySize=[arraySize])
135
136    #### write csml document ##### -move this to the csmldocument module?
137    subsetCSML=CSMLDocument()
138    subsetCSML=subsetCSML.makeGridSeries(domain,rangeSet)
139    output=open(pathToSubsetCSML,'w')
140    output.write(subsetCSML)
141    output.close()
142
143    ### create and write netcdf - uses scientific python####
144    ncfile=NetCDFFile(pathToSubsetNetCDF,'w')
145    # create the dimensions       
146    ncfile.createDimension ( 'time', len(self.times))
147    time_var = ncfile.createVariable ( 'time', Float, ('time',) )
148    time_var.longname = 'time'
149    floatTimes=[]
150    print len(fulldata[0])
151    for time in self.times:
152        time=ops_AbstractFeature.__getCDtime(time).torel(calunits)
153        floatTimes.append(time.value)
154    time_var[:] =floatTimes[:]
155    time_var.units=calunits
156    time_var.calendar=caltype
157
158    for ordinate in ordinates:
159        ncfile.createDimension(ordinate.definesAxis, len(ordinate.axisValues.split())) 
160        item_var = ncfile.createVariable (ordinate.definesAxis, Float, (ordinate.definesAxis,) )
161        #convert to list
162        vals=[]
163        for val in ordinate.axisValues.split(','):
164            vals.append(float(val))
165        print ordinate.definesAxis
166        if ordinate.definesAxis=='latitude':
167            #this is just a fix till I work out what's going wrong
168            vals.reverse()
169        print vals
170        item_var[:]=vals[:]
171    #this needs reconsidering - do the shapes always match up??
172    for ordinate in ordinates:
173        if ordinate.definesAxis in axisorder:
174            print axisorder
175            #(ordinate.definesAxis)
176        #print ordinate.definesAxis
177    #print axisorder
178   
179   
180    if len(ordinates)==3:
181        #feature_var = ncfile.createVariable (self.id, Float, ('time',ordinates[1].definesAxis,ordinates[0].definesAxis,ordinates[2].definesAxis))
182        feature_var = ncfile.createVariable (self.id, Float, ('time',axisorder[1],axisorder[2],axisorder[3]))
183    elif len(grid.ordinates)==2:
184        #feature_var = ncfile.createVariable (self.id, Float, ('time',ordinates[0].definesAxis,ordinates[1].definesAxis))
185        feature_var = ncfile.createVariable (self.id, Float, ('time',axisorder[1],axisorder[2]))
186    print shape(feature_var)
187    print shape(fulldata)
188    #print fulldata
189    feature_var[:]=fulldata[:]
190    ncfile.close()
191
192    return pathToSubsetCSML, pathToSubsetNetCDF, totalArraySize
Note: See TracBrowser for help on using the repository browser.