source: TI02-CSML/trunk/csml/csmllibs/csmldataiface.py @ 1594

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

several bug fixes to subsetting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1#!/usr/bin/env python
2
3#**************************************************************************************
4#csmldataIface.py
5#contains classes for interfacing with various files
6#currently supports cdunif (NetCDF, PP, Grib(untested)) And Nappy (NASAAmes)
7#use by instantiating the factory class: DataInterface
8#v0.00 30th November 2005
9#Dominic Lowe, BADC
10#**************************************************************************************
11
12
13import cdms 
14import nappy 
15import string
16import sys
17
18class DataInterface(object):
19        #Use DataInterface and setInterfaceType to instantiate the correct
20        #subclass for data
21        def __init__(self):
22                self.iface ='None'
23               
24        def setInterfaceType(self,interfaceType):
25                # function returns approprate data interface
26                #set interfaceType: should correspond to available datainterface subclasses
27                self.iface =interfaceType
28                if self.iface == 'nappy':
29                        return NappyInterface()
30                elif self.iface == 'cdunif':
31                        return cdunifInterface()
32               
33               
34        def getUnknownInterfaceType(self, filename):
35                #if the interface type is not known at the time of instantiation, then use
36                #this function to examine the file and return the correct interface (if it exists).
37                fileExtension = str(filename)[-3:]
38                if fileExtension == '.nc':
39                        return cdunifInterface()
40                if fileExtension == '.qxf':
41                        return cdunifInterface()
42                elif fileExtension == '.pp':
43                        return cdunifInterface()
44                elif fileExtension == 'ctl':
45                        return cdunifInterface()
46                elif fileExtension == 'xml':
47                        return cdunifInterface()
48                else:
49                        try:
50                                nappy.readFFI(filename) in [1001,1010,1020,2010,2110,2160,2310,3010,4010]
51                                return NappyInterface()
52                        except:
53                                print "Could not establish file type"
54                                print "See csmldataiface.py"
55                       
56       
57                                                       
58class AbstractDI(object):               
59        #Abstract data interface class
60        #does nothing, but contains templates for methods required for a data interface class
61        #individual interfaces (e.g NappyInterface) should override these methods
62
63        def __init__(self):
64                self.extractType=''
65                self.extractPrefix = ''
66                                       
67        def openFile(self, filename):
68                #opens file
69                self.file=''
70               
71        def closeFile(self):
72                #closes file
73                self.file.close()
74                       
75               
76               
77class NappyInterface(AbstractDI):       
78        # Data Interface for Nappy (NASA Ames Processing in Python)
79   
80        def __init__(self):
81                self.extractType='NASAAmesExtract'
82                self.extractPrefix = '_naextract_'
83                                         
84        def openFile(self, filename):
85                #print 'opening NA file: ' + str(filename)
86                self.file=nappy.openNAFile(filename)
87                #print 'reading data....'
88                #self.file.readData()
89                #print 'nappyopen ' + filename
90
91        def getListOfAxes(self):
92                axes=self.file.XNAME
93                #print 'before units stripped' + str(axes)
94                axes=self.stripunits(axes)
95                #print 'after units stripped' + str(axes)
96                return axes
97
98        def setAxis(self,axis):
99                axes = self.getListOfAxes()
100                self.axisstub=axes.index(axis)
101
102        def getAxisAttribute(self, att):
103                        #need to do something here...? maybe
104                pass
105                return attValue
106
107        def getTimeUnits(self):
108                axes = self.getListOfAxes()
109                for axis in axes:
110                        if string.find(string.upper(axis),'SECONDS SINCE') != -1:
111                                #found possible time axis.
112                                if axis[-3:]=='UTC':
113                                    units =string.lower(axis[:-4]) #hack!
114                                    units=units.replace('/','-') #try and clean it up
115                                else:
116                                    units=string.lower(axis)
117                                break
118                        elif string.find(string.upper(axis),'HOURS SINCE') != -1:
119                                #found possible time axis.
120                                units =(str(axis))
121                                break
122                        elif string.find(string.upper(axis),'DAYS SINCE') != -1:
123                                #found possible time axis.
124                                units =(str(axis))
125                                break
126                       
127                #revisit with udunits python library?
128                return units
129
130
131        def getDataForAxis(self):
132
133                if self.file.X == None:
134                        #print 'reading data....'
135                        self.file.readData()
136
137                if type(self.file.X[1])==list:
138                #if len(self.file.X) > 0:
139                        data = self.file.X[self.axisstub]
140                else:
141                        data =self.file.X
142                return data
143
144        def getSizeOfAxis(self,axis):
145
146                #check this function is okay
147                #is time always the first dimension in NA??
148                axes = self.getListOfAxes()
149                axisPosition=axes.index(axis)
150                #print "axis position" + str( axisPosition)
151                #print "NX" + str(self.file.NX)
152                try :
153                        axisSize=self.file.NX[axisPosition-1]
154                except:
155                        axisSize ='Unknown axis size'
156                return axisSize
157
158        def getListofVariables(self):
159                variableList=self.stripunits(self.file.VNAME)
160                return variableList
161
162        def setVariable(self,varname):
163                vlist=self.getListofVariables()
164                self.varstub=vlist.index(varname)
165
166        def getVariableAxes(self):
167                #hmm, now with Nasa Ames the Axis will be the same for all variables.
168                #so just use the getListOfAxes function again
169                #I think... check this!
170                varAxesList=self.getListOfAxes()
171                return varAxesList
172
173        def getVariableAttribute(self,attName):
174                if attName =='units':
175                        #strip the units (attribute) from the variable
176                        unitslist=self.getUnits(self.file.VNAME)
177                        attribValue = unitslist[self.varstub]
178                        try:
179                                attribValue = unitslist[self.varstub]
180                        except:
181                                attribValue = 'unknown'
182                else:
183                        attribValue = 'unknown'
184                return attribValue
185
186        def getDataForVar(self):
187            #NOTE TO SELF:
188            #Review this function (and in fact all of nasa ames data interface...)
189                if self.file.V == None:
190                        #print 'reading data....'
191                        self.file.readData()
192
193                try:
194                    if type(self.file.V[1])==list:
195                        data = self.file.V[self.varstub]
196                #else:
197                #       data =self.file.X
198                #       print data
199                    return data
200                except:
201                    data = self.file.X
202                   # print data
203                    return data
204
205        def getArraySizeOfVar(self):
206        #iterates through all dimensions in variable to get array size i.e a 3x3x3 array would have a size of 27
207
208                dimlist=self.file.NX
209                varsize =1
210                for item in dimlist:
211                        varsize = varsize * item
212                        #print "VARSISZE" + str(varsize)
213                return varsize
214
215        def getShapeOfVar(self):
216            #this should return a list.
217            varShape = []
218            for item in self.file.NX:
219                varShape.append(item)
220            return varShape
221
222        def getLowLimits(self):
223                lowlims = ""
224                for i in range (0, len(self.file.NX)):
225                        #for now, assume low limit is always of form 1 1 1 ..
226                        lowlims =lowlims + str(1)  +' '
227                return lowlims
228
229        def getHighLimits(self):
230                highlims = ""
231                for i in range (0, len(self.file.NX)):
232                        dimValue = self.file.NX[i]
233                        highlims =highlims  + str(dimValue) +' '
234                return highlims
235
236
237        def stripunits(self,listtostrip):
238                #strips units of measure from list
239                #eg ['Universal time (hours)', 'Altitude (km)', 'Latitude (degrees)', 'Longitude (degrees)']
240                #becomes ['Universal time', 'Altitude', 'Latitude', 'Longitude']
241                cleanlist = []
242                for item in listtostrip:
243                        openbracket=string.find(item,'(')
244                        if openbracket != -1:
245                                #if brackets exist, strip units.
246                                item=item[:openbracket-1]
247                        cleanlist.append(item)
248                return cleanlist
249
250        def getUnits(self,listwithunits):
251                #gets units from list
252                #eg ['Universal time (hours)', 'Altitude (km)', 'Latitude (degrees)', 'Longitude (degrees)']
253                #becomes ['hours', 'km', 'degrees', 'degrees']
254                unitlist=[]
255                for item in listwithunits:
256                        openbracket=string.find(item,'(')
257                        item = item[openbracket+1:-1]
258                        unitlist.append(item)
259                return unitlist
260
261        def getTimes(self):
262                #This function attempts to determine the time axis and read the time data
263                #it may well not manage it.
264                axes = self.getListOfAxes()
265                for axis in axes:
266                        if string.find(string.upper(axis),'TIME') != -1:
267                                #found possible time axis.
268                                self.setAxis(axis)
269                                times=self.getDataForAxis()
270                                break
271                        elif string.find(string.upper(axis),'SECONDS SINCE') != -1:
272                                #found possible time axis.
273                                self.setAxis(axis)
274                                times=self.getDataForAxis()
275                                break
276                        elif string.find(string.upper(axis),'HOURS SINCE') != -1:
277                                #found possible time axis.
278                                self.setAxis(axis)
279                                times=self.getDataForAxis()
280                                break
281                        elif string.find(string.upper(axis),'DAYS SINCE') != -1:
282                                #found possible time axis.
283                                self.setAxis(axis)
284                                times=self.getDataForAxis()
285                                break
286                return times
287
288
289
290class cdunifInterface(AbstractDI):
291    #Data Interface for cdunif (netcdf & pp formats & grib (not tested with grib)
292
293    def __init__(self):
294        #these are just temporary values until we can determine whether the
295        #file is netcdf pp or grib
296        self.extractType='cdunifExtract'
297        self.extractPrefix = '_cdunifextract_'
298
299    def openFile(self, filename):
300        self.file=cdms.open(filename)
301        #print 'cdunifopen ' + filename
302
303        #now we have the file name can properly determine extractType/Prefix
304        fileExtension = str(filename)[-3:]
305        if fileExtension == '.nc':
306            self.extractType = 'NetCDFExtract'
307            self.extractPrefix = '_ncextract_'
308        elif fileExtension == '.qxf':
309            self.extractType = 'NetCDFExtract'
310            self.extractPrefix = '_ncextract_'
311        elif fileExtension == '.pp':
312            self.extractType  = 'PPExtract'
313            self.extractPrefix = '_ppextract_'
314        elif fileExtension == 'ctl':
315            self.extractType = 'GRIBExtract'
316            self.extractPrefix = '_gribextract_'
317        elif fileExtension == 'xml': 
318            self.extractType = 'NetCDFExtract'  #okay this isn't true, but ok for testing
319            self.extractPrefix = '_ncextract__' 
320    def getListOfAxes(self):
321        axes=self.file.dimensions.keys()
322        return axes
323
324    def getSizeOfAxis(self,axis):
325        axisSize=self.file.dimensions[axis]
326        return axisSize
327
328    def getListofVariables(self):
329        variableList=self.file.variables.keys()
330        return variableList
331
332    def setAxis(self,axis):
333        self.axisobj=self.file.getAxis(axis)
334
335    def getAxisAttribute(self, att):
336        attValue=self.axisobj.attributes[att]
337        return attValue
338
339    def getTimeUnits(self):
340        #this does the same as getAxisAttribute, but is a separate function as different formats handle time differently.
341        return self.getAxisAttribute('units')
342
343    def getDataForAxis(self):
344        data = self.axisobj.getValue()
345        return data
346
347    def setVariable(self,varname):
348        self.varobj=self.file.variables[varname]
349
350    def getVariableAxes(self):
351        varAxesList=self.varobj.getAxisIds()
352        return varAxesList
353
354    def getVariableAttribute(self,attName):
355
356        #atts=self.varobj.getattributes(attName)
357        #print atts
358        if attName == 'long_name':
359            try:
360                attribValue = self.varobj.long_name
361            except:
362                attribValue=''
363        elif attName == 'units':
364            try:
365                attribValue = self.varobj.units
366            except:
367                attribValue=''
368
369        return attribValue
370
371    def getDataForVar(self):
372        data = self.varobj.getValue()
373        return data
374
375    def getSubsetOfDataForVar(self, **kwargs):
376        #takes keyword args defining subset eg
377        #subset=getSubsetOfDataForVar(latitude=(0.,10.0), longitude=(90, 100.0))
378        for key in kwargs:
379            if key == 'longitude':   #this test needs to be more robust...
380                if kwargs[key][0] > kwargs[key][1]:
381                    #subsetting greenwich meridian around 0
382                    lonMin = kwargs[key][0]
383                    lonMax =kwargs[key][1]
384                    kwargs[key]=(0.0, lonMax)
385                    sel=cdms.selectors.Selector(**kwargs)
386                    subset1=self.file(self.varobj.name,sel)
387                    kwargs[key]=(lonMin,359.9999)
388                    sel=cdms.selectors.Selector(**kwargs)
389                    subset2=self.file(self.varobj.name,sel)
390                    #concatenate arrays along longitude
391                    longitudeAxis=subset1.getAxisIndex('longitude') # this needs to be more robust test.
392                    subset = cdms.MV.concatenate([subset1,subset2],axis=longitudeAxis)
393                else:
394                    sel=cdms.selectors.Selector(**kwargs)
395                    subset=self.file(self.varobj.name,sel)
396        data = subset
397        return data
398
399    def getArraySizeOfVar(self):
400    #iterates through all dimensions in variable to get array size i.e a 3x3x3 array would have a size of 27
401        var = self.varobj
402        size = var.shape
403        varsize = 1
404        for item in size:
405            varsize = item *varsize
406        return varsize
407
408    def getShapeOfVar(self):
409        varShape = []
410        for item in self.varobj.shape:
411            varShape.append(item)
412        return varShape
413
414    def getLowLimits(self):
415        dimNames = self.varobj.getAxisIds()
416        lowlims = ""
417        for i in range (1, len(dimNames)):
418            #for now, assume low limit is always of form 1 1 1 ..
419            lowlims =lowlims + str(1)  +' '
420        return lowlims
421
422    def getHighLimits(self):
423        dimNames = self.varobj.getAxisIds()
424        highlims = ""
425        for i in range (1, len(dimNames)):
426            dimValue = self.file.dimensions[dimNames[i]]
427            highlims =highlims  + str(dimValue) +' '
428        return highlims
429       
430   
431
432       
Note: See TracBrowser for help on using the repository browser.