source: TI02-CSML/trunk/csml/API/ops_ProfileSeriesFeature.py @ 3096

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

subsetting of profileseries, pointseries retested and fixed

Line 
1''' ops_ProfileSeriesFeature  contains operations for ProfileSeriesFeatures'''
2import csml
3
4def testmethod(self):
5    print 'testmethod for profileSeries feature'
6    return 'testmethod profileseries'
7
8
9def getDomain(self):
10    #returns domain as a dictionary of ordinates {name: [values], ...}
11    #For ProfileSeries this is a combination of the location attribute and the time ordinate.
12    self.domain={}
13    self.gridnames={}
14    for gridOrd in self.value.profileSeriesDomain.coordTransformTable.gridOrdinates:
15        name=gridOrd.coordAxisLabel.CONTENT
16        self.gridnames[gridOrd.gridAxesSpanned.CONTENT]=name
17        #handle time axis
18        if hasattr(gridOrd.coordAxisValues, 'timePositionList'):  #..it should do..
19            if hasattr(gridOrd.coordAxisValues, 'insertedExtract'):
20                self.domain[name], fill, axisorder, units=gridOrd.timePositionList.insertedExtract.getData()
21            else:
22                vals=gridOrd.coordAxisValues.timePositionList.CONTENT             
23                valList=[]
24                for val in vals.split(): 
25                    valList.append(val)
26                self.domain[name]=valList
27
28        #handle other axes , eg pressure levels, but ignore location
29        #TODO should use lonAxis at latAxis from CRS rather than latitude, longitude.
30        if hasattr(gridOrd.coordAxisValues, 'coordinateList'): 
31            if gridOrd.coordAxisLabel not in ['latitude', 'longitude']:
32                if hasattr(gridOrd.coordAxisValues, 'insertedExtract'):
33                    self.domain[name], fill, axisorder, units=gridOrd.coordinateList.insertedExtract.getData()
34                else:
35                    vals=gridOrd.coordAxisValues.coordinateList.CONTENT             
36                    valList=[]
37                    for val in vals.split(): 
38                        valList.append(float(val))
39                    self.domain[name]=valList
40    return self.domain
41
42def _subsetProfileSeries(self, **kwargs):
43    '''this takes a temporal selection from a ProfileSeries object  - i.e. returns one or more 'profiles', can be used to make another ProfileSeriesFeature or a Profile feature in the case of a single selection'''
44           
45    #set self.domain:
46    self.getDomain()     
47   
48    #get the CRS from a  the  catalogue
49    cat=csml.csmllibs.csmlcrs.CRSCatalogue()
50    crs=cat.getCRS(self.value.profileSeriesDomain.srsName)
51   
52    #non-feature specific setup code, mainly handles the time dimension/calendar
53    pathToSubsetNetCDF, kwargs, timeAxis, timeName, caltype, times=csml.API.genSubset.genericSubset(self, self.outputdir,self.ncname, self.domain, kwargs)
54     
55     
56    ##Get names of variables in file and relate them to the subset selection
57    selection={}
58   
59   
60    for gridOrd in self.value.profileSeriesDomain.coordTransformTable.gridOrdinates:
61        try:
62            selection[gridOrd.gridAxesSpanned.CONTENT]=kwargs[gridOrd.coordAxisLabel.CONTENT]
63        except KeyError:
64            allValues=tuple(self.domain[gridOrd.coordAxisLabel.CONTENT])
65    strTimes, axisorder, units, fulldata, fillvalue =csml.API.genSubset.getTheData(self, selection, times, timeName, timeAxis)
66    return pathToSubsetNetCDF, crs, timeName, times, strTimes, caltype, axisorder,units, fulldata, fillvalue, kwargs
67           
68           
69def subsetToProfile(self, outputdir=None,ncname ='profile.nc',**kwargs):
70    #perform the subset
71    if outputdir is not None:
72        self.outputdir=outputdir
73    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
74        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
75    self.ncname=ncname
76    pathToSubsetNetCDF, crs, timeName, times, strTimes,  caltype, axisorder,units, fulldata, fillvalue, kwargs=self._subsetProfileSeries(**kwargs) 
77   
78    try:
79        lat=crs.axes[crs.latAxis]
80    except:
81        lat='latitude'
82    try:
83        lon=crs.axes[crs.lonAxis]
84    except:
85        lon='longitude'
86    time=crs.axes[crs.timeAxis]
87   
88    for axis in axisorder:
89        if axis not in [time, lat, lon]:
90            #this should be the 'profile' axis
91            if len(self.domain[axis]) ==1:
92                prof=self.domain[axis][0]
93            else:
94                prof=self.domain[axis]
95           
96            profName=axis
97   
98   
99    #get the totalArraySize
100    domainSubset, totalArraySize=csml.API.genSubset.subsetDomain(timeName,strTimes,self.domain, **kwargs)
101
102    #Okay, got the data now. Need to write CSML feature and NetCDF files.
103    #Writing out the CSML feature
104    # define domain/coverage  to use in 'value' attribute   
105   
106    newdomain=csml.parser.csString(prof)
107   
108    #rangeSet.arrayDescriptor=csml.parser.NetCDFExtract(id=self.id,fileName=csml.parser.csString(pathToSubsetNetCDF),variableName=csml.parser.csString(self.id),arraySize=csml.parser.csString(totalArraySize))
109    descid=csml.csmllibs.csmlextra.getRandomID()
110    descriptor=csml.parser.NetCDFExtract(id=descid,fileName=csml.parser.csString(ncname),variableName=self.name,arraySize=csml.parser.csString(totalArraySize))
111    rs=csml.parser.RangeSet()
112    va=csml.parser.ValueArray()
113    vc=csml.parser.MeasureOrNullList()
114    vc.href='#%s'%descid
115    vc.arcrole="http://ndg.nerc.ac.uk/xlinkUsage/insert#QuantityList"
116    vc.role="http://ndg.nerc.ac.uk/fileFormat/csmlStorageDescriptor"
117    vc.show='embed'
118    vc.insertedExtract=descriptor
119    try:
120        vc.uom=self.value.rangeSet.valueArray.valueComponent.uom
121    except:
122        vc.uom = 'unknown'  #TODO,
123   
124    va.valueComponent=vc
125    va.id=csml.csmllibs.csmlextra.getRandomID()
126    rs.valueArray=va
127   
128
129   
130   
131   
132    #profileseries coverage
133    cvg=csml.parser.ProfileSeriesCoverage()
134    cvg.rangeSet=rs
135    cvg.profileSeriesDomain=newdomain   
136   
137    #parameter, as before subsetting.
138    param = self.parameter
139    stdname=param.getStandardName()
140       
141    #create 'location' attribute:
142    loc=self.location  #locataion attribute of profile is same as that of parent profileseries
143   
144    #create 'time' attribute
145    tName=crs.axes[crs.timeAxis]
146    try:
147        t=kwargs[tName]
148    except:
149        t=csml.parser.csString('unknown')
150               
151    #create a stand alone profile feature containing this coverage
152    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
153    subsettedFeature=csmlWrap.createProfileFeature(value=cvg,parameter=param,location=loc, time=t, featureID=self.id,description=self.description)
154 
155    ### write netcdf using NCWriter class (wraps cdms) ###
156    nc=csml.csmllibs.netCDFWriter.NCwriter(pathToSubsetNetCDF)
157    axislist=[]
158    for a in axisorder:
159        if self.gridnames.has_key(a):
160            axislist.append(self.gridnames[a])
161        else:
162            axislist.append(a)
163    ords=None 
164    #need to pass in all dimensions except time as the kwargs..
165    otherdims={}
166    otherdims[profName]=prof
167    otherdims['latitude']=float(loc.CONTENT.split()[0])
168    otherdims['longitude']=float(loc.CONTENT.split()[1])
169    nc.genWriteVar(self.name.CONTENT,ords, times,  caltype, axislist, units, stdname, fulldata, fillvalue, **otherdims)
170    nc.closeFinishedFile()
171    print 'NetCDF file written to %s'%pathToSubsetNetCDF
172    return subsettedFeature, pathToSubsetNetCDF,descriptor
173
174def subsetToPointSeries(self, outputdir=None,  ncname='pointseries.nc',**kwargs):
175    #this operation needs more testing.
176    #perform the subset
177    if outputdir is not None:
178        self.outputdir=outputdir
179    elif csml.API.csmlContainer.globalContainer.OUTPUTDIR is not None:
180        self.outputdir=csml.API.csmlContainer.globalContainer.OUTPUTDIR
181    self.ncname=ncname
182   
183    pathToSubsetNetCDF, crs, timeName, times, strTimes, caltype, axisorder,units, fulldata, fillvalue, kwargs=self._subsetProfileSeries(**kwargs)   
184    time=crs.axes[crs.timeAxis]
185    lon='longitude'
186    lat='latitude'
187    for axis in axisorder:
188        if axis not in [time, lat, lon]:
189            #another axis
190            if len(self.domain[axis]) ==1:
191                extraAxis=self.domain[axis][0]
192            else:
193                extraAxis=self.domain[axis]
194            extraName=axis
195   
196    newdomain=csml.parser.TimeSeries()
197    newdomain.timePositionList=csml.parser.csString(csml.csmllibs.csmlextra.cleanString(str(times)))
198   
199    rangeSet=csml.parser.RangeSet()
200    rangeSet.arrayDescriptor=csml.parser.NetCDFExtract(id=self.id,fileName=csml.parser.csString(pathToSubsetNetCDF),variableName=csml.parser.csString(self.id),arraySize=csml.parser.csString(len(times)))
201   
202    #gridseries coverage
203    cvg=csml.parser.PointSeriesCoverage()
204    cvg.rangeSet=rangeSet
205    cvg.pointSeriesDomain=newdomain   
206   
207    #parameter, as before subsetting.
208    param = self.parameter
209    stdname=param.getStandardName()
210       
211    #create 'location' attribute:
212    loc=self.location  #locataion attribute of profile is same as that of parent profileseries
213   
214               
215    #create a stand alone point series feature containing this coverage
216    csmlWrap=csml.csmllibs.csmlfeaturewrap.CSMLWrapper()
217    subsettedFeature=csmlWrap.createPointSeriesFeature(value=cvg,parameter=param,location=loc, featureID=self.id,description=self.description)
218 
219    ### write netcdf using NCWriter class (wraps cdms) ###
220    nc=csml.csmllibs.netCDFWriter.NCwriter(pathToSubsetNetCDF)
221    axislist=[]
222    for a in axisorder:
223        if self.gridnames.has_key(a):
224            axislist.append(self.gridnames[a])
225        else:
226            axislist.append(a)
227    ords=None
228    #need to pass in all dimensions except time as the kwargs..
229    otherdims={}
230    otherdims[extraName]=extraAxis
231    otherdims['latitude']=float(loc.CONTENT.split()[0])
232    otherdims['longitude']=float(loc.CONTENT.split()[1])   
233    nc.genWriteVar(self.name.CONTENT,ords, times, caltype, axislist, units, stdname,fulldata, fillvalue, **otherdims)
234    nc.closeFinishedFile()
235    print 'NetCDF file written to %s'%pathToSubsetNetCDF
236    return subsettedFeature, pathToSubsetNetCDF, rangeSet.arrayDescriptor
237
Note: See TracBrowser for help on using the repository browser.