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

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

more changes related to module reorganisation

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