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

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

more module moves
Line 
1''' ops_GridSeriesFeature  contains operations for GridSeriesFeatures'''
2
3from Numeric import *
4import csml.parser
5import csml.csmllibs.csmltime
6import csml.csmllibs.csmlDocument
7import csml.API.ops_AbstractFeature
8
9
10def testmethod(self):
11    #print 'testmethod for gridseries feature'
12    return 'testmethod - gridseries'
13
14def getAllowedSubsettings(self):
15    return ['subsetToGridSeries']
16
17def getDomainReference(self):
18    #Inspects a time position list for the domain reference.
19    #TODO, does not handle a file extract in place of a list.
20    if isinstance(self.domain.domainReference,csml.parser.TimePositionList):
21        time = {}
22        time['t'] = self.domain.domainReference.timePositions
23        domainref  = time
24    return domainref
25
26def getDomainComplement(self):
27    #This will return a list containing one or more ordinates:
28    #currently in form [Name, values]
29    #assumes ordinate.axisValues is a file extract
30    #TODO axisValues may be inline -see also domainReference for similar problem
31    domaincomp ={}
32    dc = self.domain.domainComplement
33    #dc should be a grid!
34    if isinstance(dc, csml.parser.Grid):
35        for ordinate in dc.ordinates:
36            domaincomp[ordinate.definesAxis]=ordinate.axisValues.getData()
37    return domaincomp
38
39def getDomain(self):
40    #returns both the domain reference axes and domain compliment axes in a single domain dictionary
41    #axes are in no particular order
42    domain = {}
43    dr=getDomainReference(self)
44    dc=getDomainComplement(self)
45    for key in dc.keys():
46        domain[key]=dc[key]
47    for key in dr.keys():
48        domain[key]=dr[key]
49    return domain
50
51
52
53
54def subsetToGridSeries(self, timeSubset,  csmlpath=None, ncpath=None,**kwargs):
55    if csmlpath is not None:
56        pathToSubsetCSML = csmlpath
57    else:
58        pathToSubsetCSML='temp.xml'
59    if ncpath is not None:
60        pathToSubsetNetCDF=ncpath
61    else:
62        pathToSubsetNetCDF='temp.nc'
63   
64    #deal with longitude requests
65    #if the request is in -ve,+ve eg (-30,30) but the data is in (0,360) need to handle this by changing the args.
66    dc = self.getDomainComplement()
67    for key in dc.keys():
68        if key == 'longitude': #how do we test if it is longitude properly?
69            print dc[key]
70            for val in dc[key]:
71                if val < 0:
72                    pass
73                else:
74                    if kwargs[key][0] < 0:
75                        kwargs[key]=(kwargs[key][0]+360,kwargs[key][1])
76                    if kwargs[key][1] < 0:
77                        kwargs[key]=(kwargs[key][0],kwargs[key][1]+360)
78         
79    domainref = getDomainReference(self)
80    self.times=timeSubset
81    self.files=[]
82    strTimes=''
83    fulldata=[]
84    if len(self.times) == 2:
85        try:
86            tone=ops_AbstractFeature.__getCDtime(self.times[0])
87        except:
88            tone=self.times[0]
89        try:
90            ttwo=ops_AbstractFeature.__getCDtime(self.times[1])
91        except:
92            ttwo=self.times[1]
93        dr=getDomainReference(self)
94        self.times=[]
95        for time in dr['t'].split():
96            timeok=csml.API.ops_AbstractFeature.__compareTimes(tone,time,ttwo)
97            if timeok ==1:
98                self.times.append(time)
99   
100    #to keep track of files that have already been fetched. eg. if multiple times are in a single file only need to get data from that file once...
101    numFiles=len(self.rangeSet.aggregatedArray.component[0].fileName.split())
102    timeToFileRatio=len(domainref['t'].split())/numFiles
103    filesFetched=[]
104    #get data:
105    for time in self.times:
106        listPosition=domainref['t'].split().index(time)
107        strTimes= strTimes + ' ' + time
108        timedim=self.rangeSet.aggregatedArray.component[0].variableName
109        for comp in self.rangeSet.aggregatedArray.component:
110            filePos=(listPosition)/timeToFileRatio
111            if filePos in filesFetched:
112                continue #already got data from this file, try next time
113            data=comp.getData(fileposition=filePos, times=self.times, **kwargs)
114            self.files.append(comp.fileName.split()[filePos])
115            if fulldata ==[]:
116               fulldata = data.tolist()
117            else:
118                for item in data.tolist():
119                    fulldata.append(item)
120            filesFetched.append(filePos)
121        axisorder = data.getAxisIds()  #will need later!
122        print 'filesFetched %s'%filesFetched
123    try:
124        caltype=self.domain.domainReference.frame.split(':',1)[0]
125        calunits=self.domain.domainReference.frame.split(':',1)[1]
126        csml.csmllibs.csmltime.setcdtimeCalendar(caltype)
127    except:
128        csml.csmllibs.csmltime.setcdtimeCalendar(csmltime.cdtime.DefaultCalendar)
129    ### define domain and rangeSet to use for feature in csml document####
130    domain=csml.parser.GridSeriesDomain()
131    domain.domainReference=csml.parser.TimePositionList(timePositions=strTimes)
132    grid=csml.parser.Grid()
133    #dc = self.getDomainComplement()
134    ordinates= []
135    i=0
136    valueStore=[]  # use the values again later to generate netcdf
137    arraySize=0
138    totalArraySize=1
139    for key in dc.keys():
140        arraySize=0
141        i=i+1
142        god=csml.parser.GridOrdinateDescription()
143        god.gridAxesSpanned='dim%s'%i
144        god.sequenceRule='+x+y+z'
145        god.definesAxis=key
146        straxisValues=''
147       
148        #now deal with each argument:
149        if key in kwargs:
150            if kwargs[key][0] < kwargs[key][1]:   
151                for val in dc[key]:
152                    if val >= kwargs[key][0]:
153                        if val <= kwargs[key] [1]:
154                            arraySize=arraySize+1
155                            straxisValues=straxisValues+ str(val) + ', '
156            else:#this if deals with selections such as longitude (330,30) where the lower limit is 'greater' than the upper limit in a mathematical sense.
157                    for val in dc[key]:
158                        if val >= kwargs[key][0]:
159                            arraySize=arraySize+1
160                            straxisValues=straxisValues+ str(val) + ', '
161                    for val in dc[key]:
162                        if val <= kwargs[key] [1]:
163                            arraySize=arraySize+1
164                            straxisValues=straxisValues+ str(val) + ', '
165        else: # this dimension has not been subsetted at all
166            for val in dc[key]:
167                arraySize=arraySize+1
168                straxisValues=straxisValues+ str(val) + ', '
169        totalArraySize=totalArraySize*arraySize
170        god.axisValues=straxisValues[:-2]
171        ordinates.append(god)
172    totalArraySize=totalArraySize*len(self.times)
173    grid.ordinates=ordinates
174    domain.domainComplement=grid
175    rangeSet=csml.parser.RangeSet()
176    rangeSet.arrayDescriptor=csml.parser.NetCDFExtract(id=self.id,fileName=pathToSubsetNetCDF,variableName=self.id,arraySize=[arraySize])
177
178    #### write csml document #####
179    csmldoc=csml.csmllibs.csmlDocument.CSMLDocument("mydoc", "mymetadata")
180    csmldoc.addGridSeriesFeature(domain,rangeSet,datasetID="A",featureID="B",description="C")
181    csmldoc=csmldoc.consolidate()
182    output=open(pathToSubsetCSML,'w')
183    output.write(csmldoc)
184    output.close()
185
186    ### write netcdf using NCWriter class (wraps cdms) ###
187    nc=NetCDFWriter.NCwriter(pathToSubsetNetCDF)
188    floatTimes=[]
189    for time in self.times:
190        time=ops_AbstractFeature.__getCDtime(time).torel(calunits)
191        floatTimes.append(time.value)
192    nc.addAxis('t',floatTimes,isTime=1,units=calunits,calendar=caltype)
193    for ordinate in ordinates:
194        lon,lat=None,None
195        if ordinate.definesAxis=='longitude':
196            lon=1
197        if ordinate.definesAxis=='latitude':
198            lat=1
199        #convert to list
200        vals=[]
201        for val in ordinate.axisValues.split(','):
202            vals.append(float(val))
203        nc.addAxis(ordinate.definesAxis,vals,isLon=lon,isLat=lat,units='')#to do, units attribute for CF compliance
204    if len(ordinates)==3:
205        axes=['t',axisorder[1],axisorder[2],axisorder[3]]
206    elif len(grid.ordinates)==2:
207        axes=['t',axisorder[1],axisorder[2]]
208    nc.addVariable(fulldata,self.id, axes,units='') #to do, units attribute for CF compliance
209    nc.closeFinishedFile()
210    return pathToSubsetCSML,pathToSubsetNetCDF, totalArraySize
Note: See TracBrowser for help on using the repository browser.