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

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

much of the generic subsetting code moved out of ops_gridseries into genSubset.py

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