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

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

more restructuring of subsetting code

Line 
1''' ops_GridSeriesFeature  contains operations for GridSeriesFeatures'''
2
3import csml.parser
4import csml.csmllibs.csmltime
5import csml.csmllibs.csmldocument
6import csml.API.ops_AbstractFeature
7import csml.API.genSubset
8import csml.csmllibs.netCDFWriter
9import csmlutils
10
11
12import sys  #remove later
13
14def testmethod(self):
15    #print 'testmethod for gridseries feature'
16    return 'testmethod - gridseries'
17
18def getAllowedSubsettings(self):
19    return ['subsetToGridSeries']  #other operations
20
21def getDomain(self):
22    #returns domain as a dictionary of ordinates {name: [values], ...}
23    self.domain={}
24    for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
25        name=gridOrd.coordAxisLabel.CONTENT
26        if hasattr(gridOrd.coordAxisValues, '__insertedExtract'):
27            self.domain[name]=gridOrd.coordAxisValues.__insertedExtract.getData()
28        else:
29            vals=gridOrd.coordAxisValues.coordinateList.CONTENT
30            valList=[]
31            for val in vals.split(): 
32                valList.append(val)
33            self.domain[name]=valList
34    return self.domain
35
36   
37def subsetToGridSeries(self, csmlpath=None, ncpath=None,**kwargs):
38    #set self.domain:
39    self.getDomain()
40   
41    #get the CRS from a  the  catalogue
42    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
43    crs=cat.getCRS(self.value.gridSeriesDomain.srsName) 
44   
45    #non-feature specific setup code, mainly handles the time dimension/calendar
46    pathToSubsetNetCDF, kwargs, timeName, calunits, caltype, times=csml.API.genSubset._genericSubset(self, csmlpath, ncpath, self.domain, kwargs)
47     
48    #Find the original time name (timeAxis) for the file.
49    for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
50        if gridOrd.coordAxisLabel.CONTENT==timeName:
51            timeAxis = gridOrd.gridAxesSpanned.CONTENT  #only works for regular grids atm
52        else: 
53            print 'crs not supported'
54             
55
56   #get the ratio of times to files
57    timeToFileRatio = csml.API.genSubset._getTimeToFileRatio(self, self.domain, timeName)
58       
59           
60    #list 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...
61    filesFetched=[]
62       
63    ##Get names of variables in file and relate them to the subset selection
64    selection={}
65    for gridOrd in self.value.gridSeriesDomain.coordTransformTable.gridOrdinates:
66        selection[gridOrd.gridAxesSpanned.CONTENT]=kwargs[gridOrd.coordAxisLabel.CONTENT]
67   
68   
69    self.files=[]
70    strTimes=''
71    fulldata=[]
72   
73   
74    for time in times:
75        listPosition=self.domain[timeName].index(time)
76        strTimes= strTimes + ' ' + time
77        for comp in csmlutils.listify(self.value.rangeSet.valueArray.valueComponent.quantityList.__insertedExtract.components): 
78            filePos=int(float(listPosition)/timeToFileRatio)
79            if filePos in filesFetched:
80                continue #already got data from this file, try next time
81            print selection
82            data=comp.getData(fileposition=filePos, **selection)
83            print data
84           
85               
86            self.files.append(comp.fileList.fileNames.CONTENT[filePos]) #TODO, get the right file name
87            if fulldata ==[]:
88                fulldata = data.tolist()
89                print fulldata
90            else:
91                for item in data.tolist():
92                    fulldata.append(item)
93            filesFetched.append(filePos)
94        axisorder = data.getAxisIds()  #will need later!
95
96    #elif hasattr(self.value.rangeSet, 'datablock'): #not tested
97        #pass
98   
99    #Okay, got the data now. Need to write out CSML and NetCDF files...
100 
101
102    #Writing out the CSML file
103    # define domain/coverage  to use in 'value' attribute   
104   
105    domain=csml.parser.GridSeriesDomain()
106    grid=csml.parser.GridCoordinatesTable()
107    domainSubset, totalArraySize=csml.API.genSubset._subsetDomain(timeName,strTimes,self.domain, **kwargs)
108    cTT=csml.parser.GridCoordinatesTable()
109    ords =[]
110    for key in domainSubset.keys():
111        go=csml.parser.GridOrdinateDescription()
112        go.coordAxisLabel=csml.parser.csString(key)
113        go.gridAxesSpanned=csml.parser.csString(key)
114        go.coordAxisValues = csml.parser.SpatialOrTemporalPositionList()
115        if key==crs.axes[crs.timeAxis]:
116            go.coordAxisValues.timePositionList=csml.parser.csString(domainSubset[key]) #self.domain.key placeholder
117        else:
118            go.coordAxisValues.coordinateList=csml.parser.csString(domainSubset[key]) #self.domain.key placeholder
119        ords.append(go) #go needs a few more properties setting
120    cTT.gridOrdinates=ords
121    grid=self.value.gridSeriesDomain.coordTransformTable
122    domain.coordTransformTable=cTT
123    rangeSet=csml.parser.RangeSet()
124    rangeSet.arrayDescriptor=csml.parser.NetCDFExtract(id=self.id,fileName=csml.parser.csString(pathToSubsetNetCDF),variableName=csml.parser.csString(self.id),arraySize=csml.parser.csString(totalArraySize))
125    cvg=csml.parser.GridSeriesCoverage()
126    cvg.rangeSet=rangeSet
127    cvg.gridSeriesDomain=domain   
128    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
129    subsettedFeature=csmlWrap.createGridSeriesFeature(value=cvg,datasetID="A",featureID="B",description="C")
130   
131   
132    ### write netcdf using NCWriter class (wraps cdms) ###
133    nc=csml.csmllibs.netCDFWriter.NCwriter(pathToSubsetNetCDF)
134    floatTimes=[]
135    for time in times:
136        time=csml.API.ops_AbstractFeature.__getCDtime(time).torel(calunits)
137        floatTimes.append(time.value)
138    nc.addAxis('t',floatTimes,isTime=1,units=calunits,calendar=caltype)
139    for ord in ords:
140        print ord.coordAxisLabel.CONTENT
141        vals=[]
142        lon,lat=None,None
143        if ord.coordAxisLabel.CONTENT=='Time':
144            continue
145        else:
146            for val in ord.coordAxisValues.coordinateList.CONTENT.split(','):
147                if val != ' ':
148                    vals.append(float(val)) 
149                    print vals
150        if ord.coordAxisLabel.CONTENT=='Lon':
151            lon=1
152            name='longitude'
153        elif ord.coordAxisLabel.CONTENT=='Lat':
154            lat=1
155            name='latitude'
156        else:
157            name=ord.coordAxisLabel.CONTENT
158 
159        nc.addAxis(name,vals,isLon=lon,isLat=lat,units='')#to do, units attribute for CF compliance
160    if len(ords)==3:
161        axes=['t',axisorder[1],axisorder[2]]
162    elif len(ords)==2:
163        axes=['t',axisorder[1]]
164    #print fulldata
165    print axes
166    nc.addVariable(fulldata,self.id, axes,units='') #to do, units attribute for CF compliance
167    print 'added'
168    nc.closeFinishedFile()
169    return subsettedFeature, pathToSubsetNetCDF
170   
Note: See TracBrowser for help on using the repository browser.