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

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

subsetting by index initially working for rawfileextracts (more testing needed)

  • 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
12import pdb
13import cdms 
14import Numeric
15
16try:
17    import nappy 
18except ImportError:
19    pass
20    #print 'could not import NASAAmes interface'
21import string
22import sys
23import csml.csmllibs.csmltime
24# This is required to prevent Numeric arrays being truncated when printed.
25import MA
26MA.set_print_limit(0)
27
28# Generic mathematical expression solver, required by the raw and image
29# interfaces
30import NumericTransform
31
32try:
33    # Part of the PIL. Required for ImageFileInterface:
34    import Image
35except ImportError:
36    print 'could not import Image'
37    pass
38
39try:
40    # Required for RawFileInterface:
41    import struct
42except ImportError:
43    print 'could not import struct'
44    pass
45
46
47class DataInterface(object):
48        #Use DataInterface and setInterfaceType to instantiate the correct
49        #subclass for data
50        def __init__(self):
51                self.iface ='None'
52               
53        def setInterfaceType(self,interfaceType):
54                # function returns approprate data interface
55                #set interfaceType: should correspond to available datainterface subclasses
56                self.iface =interfaceType
57                if self.iface == 'nappy':
58                        return NappyInterface()
59                elif self.iface == 'cdunif':
60                        return cdunifInterface()
61                elif self.iface == 'raw':
62                        return RawFileInterface()
63                elif self.iface == 'image':
64                        return ImageFileInterface()
65               
66               
67        def getUnknownInterfaceType(self, filename):
68                #if the interface type is not known at the time of instantiation, then use
69                #this function to examine the file and return the correct interface (if it exists).
70                fileExtension = str(filename)[-3:]
71                if fileExtension == '.nc':
72                        return cdunifInterface()
73                if fileExtension == '.qxf':
74                        return cdunifInterface()
75                elif fileExtension == '.pp':
76                        return cdunifInterface()
77                elif fileExtension == 'ctl':
78                        return cdunifInterface()
79                elif fileExtension == 'xml':
80                        return cdmlInterface()
81                elif fileExtensions == 'png' or fileExtension == 'gif':
82                        return ImageFileInterface()
83                else:
84                        try:
85                                nappy.readFFI(filename) in [1001,1010,1020,2010,2110,2160,2310,3010,4010]
86                                return NappyInterface()
87                        except:
88                                print "Could not establish file type"
89                                print "See csmldataiface.py"
90                       
91       
92                                                       
93class AbstractDI(object):               
94        #Abstract data interface class
95        #does nothing, but contains templates for methods required for a data interface class
96        #individual interfaces (e.g NappyInterface) should override these methods
97
98        def __init__(self):
99                self.extractType=''
100                self.extractPrefix = ''
101                                       
102        def openFile(self, filename):
103           #opens file, must be overwritten by subclass
104            raise NotImplementedError 
105               
106        def closeFile(self):
107           #closes file, probably needs to be overwritten by subclass
108            try:
109                self.file.close()
110            except:
111                raise NotImplementedError
112         
113        def setAxis(self,axis):
114           #'set' the name of the current axis , must be overwritten by subclass
115           #this may just involve a stub (see NASAAmes interface) or may involve
116           #calling a real set method of the underlying api (see cdunif Interface)
117           raise NotImplementedError 
118                       
119        def getDataForAxis(self):
120            #return all data for axis, must be overwritten by subclass
121           raise NotImplementedError   
122       
123        def setVariable(self,varname):
124             #As for setAxis, 'set' the name of the current axis , must be overwritten by subclass
125             raise NotImplementedError
126       
127        def getDataForVar(self):
128            #return all data for variable, must be overwritten by subclass
129           raise NotImplementedError   
130       
131        def getSubsetOfDataForVar(self,**kwargs):
132            #return subset of data for variable, must be overwritten by subclass
133            raise NotImplementedError
134       
135       
136class NappyInterface(AbstractDI):       
137        # Data Interface for Nappy (NASA Ames Processing in Python)
138   
139        def __init__(self):
140                self.extractType='NASAAmesExtract'
141                self.extractPrefix = '_naextract_'
142                                         
143        def openFile(self, filename):
144                #print 'opening NA file: ' + str(filename)
145                self.file=nappy.openNAFile(filename)
146                #print 'reading data....'
147                #self.file.readData()
148                #print 'nappyopen ' + filename
149
150        def getListOfAxes(self):
151                axes=self.file.XNAME
152                #print 'before units stripped' + str(axes)
153                axes=self.stripunits(axes)
154                #print 'after units stripped' + str(axes)
155                return axes
156
157        def setAxis(self,axis):
158                axes = self.getListOfAxes()
159                self.axisstub=axes.index(axis)
160
161        def getAxisAttribute(self, att):
162                #need to do something here...? maybe
163                attValue=None
164                return attValue
165
166        def getTimeUnits(self):
167                axes = self.getListOfAxes()
168                for axis in axes:
169                        if string.find(string.upper(axis),'SECONDS SINCE') != -1:
170                                #found possible time axis.
171                                if axis[-3:]=='UTC':
172                                    units =string.lower(axis[:-4]) #hack!
173                                    units=units.replace('/','-') #try and clean it up
174                                else:
175                                    units=string.lower(axis)
176                                break
177                        elif string.find(string.upper(axis),'HOURS SINCE') != -1:
178                                #found possible time axis.
179                                units =(str(axis))
180                                break
181                        elif string.find(string.upper(axis),'DAYS SINCE') != -1:
182                                #found possible time axis.
183                                units =(str(axis))
184                                break
185                       
186                #revisit with udunits python library?
187                return units
188
189
190        def getDataForAxis(self):           
191            if self.file.X == None:
192                    #print 'reading data....'
193                    self.file.readData()
194            for x in self.file.X:
195                print type(x)
196            print x
197           
198            d=Numeric.array(self.file.X)
199            print d
200            if type(self.file.X[1])==list:
201            #if len(self.file.X) > 0:
202                    print '1'
203                    data = self.file.X[self.axisstub]
204            else:
205                    print '2'
206                    data =self.file.X
207            #print data
208            data=Numeric.array(data)
209            return data
210
211        def getSizeOfAxis(self,axis):
212
213                #check this function is okay
214                #is time always the first dimension in NA??
215                axes = self.getListOfAxes()
216                axisPosition=axes.index(axis)
217                #print "axis position" + str( axisPosition)
218                #print "NX" + str(self.file.NX)
219                try :
220                        axisSize=self.file.NX[axisPosition-1]
221                except:
222                        axisSize ='Unknown axis size'
223                return axisSize
224
225        def getListofVariables(self):
226                variableList=self.stripunits(self.file.VNAME)
227                return variableList
228
229        def setVariable(self,varname):
230                vlist=self.getListofVariables()
231                self.varstub=vlist.index(varname)
232
233        def getVariableAxes(self):
234                #hmm, now with Nasa Ames the Axis will be the same for all variables.
235                #so just use the getListOfAxes function again
236                #I think... check this!
237                varAxesList=self.getListOfAxes()
238                return varAxesList
239
240        def getVariableAttribute(self,attName):
241                if attName =='units':
242                        #strip the units (attribute) from the variable
243                        unitslist=self.getUnits(self.file.VNAME)
244                        attribValue = unitslist[self.varstub]
245                        try:
246                                attribValue = unitslist[self.varstub]
247                        except:
248                                attribValue = 'unknown'
249                else:
250                        attribValue = 'unknown'
251                return attribValue
252
253        def getDataForVar(self):
254            #NOTE TO SELF:
255            #Review this function (and in fact all of nasa ames data interface...)
256                if self.file.V == None:
257                        #print 'reading data....'
258                        self.file.readData()
259
260                try:
261                    if type(self.file.V[1])==list:
262                        data = self.file.V[self.varstub]
263                #else:
264                #       data =self.file.X
265                #       print data
266                    return data
267                except:
268                    data = self.file.X
269                   # print data
270                    return data
271
272        def getArraySizeOfVar(self):
273        #iterates through all dimensions in variable to get array size i.e a 3x3x3 array would have a size of 27
274
275                dimlist=self.file.NX
276                varsize =1
277                for item in dimlist:
278                        varsize = varsize * item
279                        #print "VARSISZE" + str(varsize)
280                return varsize
281
282        def getShapeOfVar(self):
283            #this should return a list.
284            varShape = []
285            for item in self.file.NX:
286                varShape.append(item)
287            return varShape
288
289        def getLowLimits(self):
290                lowlims = ""
291                for i in range (0, len(self.file.NX)):
292                        #for now, assume low limit is always of form 1 1 1 ..
293                        lowlims =lowlims + str(1)  +' '
294                return lowlims
295
296        def getHighLimits(self):
297                highlims = ""
298                for i in range (0, len(self.file.NX)):
299                        dimValue = self.file.NX[i]
300                        highlims =highlims  + str(dimValue) +' '
301                return highlims
302
303
304        def stripunits(self,listtostrip):
305                #strips units of measure from list
306                #eg ['Universal time (hours)', 'Altitude (km)', 'Latitude (degrees)', 'Longitude (degrees)']
307                #becomes ['Universal time', 'Altitude', 'Latitude', 'Longitude']
308                cleanlist = []
309                for item in listtostrip:
310                        openbracket=string.find(item,'(')
311                        if openbracket != -1:
312                                #if brackets exist, strip units.
313                                item=item[:openbracket-1]
314                        cleanlist.append(item)
315                return cleanlist
316
317        def getUnits(self,listwithunits):
318                #gets units from list
319                #eg ['Universal time (hours)', 'Altitude (km)', 'Latitude (degrees)', 'Longitude (degrees)']
320                #becomes ['hours', 'km', 'degrees', 'degrees']
321                unitlist=[]
322                for item in listwithunits:
323                        openbracket=string.find(item,'(')
324                        item = item[openbracket+1:-1]
325                        unitlist.append(item)
326                return unitlist
327
328        def getTimes(self):
329                print 'HELLO'
330                #This function attempts to determine the time axis and read the time data
331                #it may well not manage it.
332                axes = self.getListOfAxes()
333                for axis in axes:
334                        if string.find(string.upper(axis),'TIME') != -1:
335                                #found possible time axis.
336                                self.setAxis(axis)
337                                times=self.getDataForAxis()
338                                break
339                        elif string.find(string.upper(axis),'SECONDS SINCE') != -1:
340                                #found possible time axis.
341                                self.setAxis(axis)
342                                times=self.getDataForAxis()
343                                break
344                        elif string.find(string.upper(axis),'HOURS SINCE') != -1:
345                                #found possible time axis.
346                                self.setAxis(axis)
347                                times=self.getDataForAxis()
348                                break
349                        elif string.find(string.upper(axis),'DAYS SINCE') != -1:
350                                #found possible time axis.
351                                self.setAxis(axis)
352                                times=self.getDataForAxis()
353                                break
354                times=Numeric.array(times)
355                return times
356
357
358
359class cdunifInterface(AbstractDI):
360    #Data Interface for cdunif (netcdf & pp formats & grib (not tested with grib)
361
362    def __init__(self):
363        #these are just temporary values until we can determine whether the
364        #file is netcdf pp or grib
365        self.extractType='cdunifExtract'
366        self.extractPrefix = '_cdunifextract_'
367
368    def openFile(self, filename):
369        self.file=cdms.open(filename)
370
371        #now we have the file name can properly determine extractType/Prefix
372        fileExtension = str(filename)[-3:]
373        if fileExtension == '.nc':
374            self.extractType = 'NetCDFExtract'
375            self.extractPrefix = '_ncextract_'
376        elif fileExtension == '.qxf':
377            self.extractType = 'NetCDFExtract'
378            self.extractPrefix = '_ncextract_'
379        elif fileExtension == '.pp':
380            self.extractType  = 'PPExtract'
381            self.extractPrefix = '_ppextract_'
382        elif fileExtension == 'ctl':
383            self.extractType = 'GRIBExtract'
384            self.extractPrefix = '_gribextract_'
385        elif fileExtension == 'xml': 
386            self.extractType = 'NetCDFExtract'  #okay this isn't true, but ok for testing
387            self.extractPrefix = '_ncextract__' 
388    def getListOfAxes(self):
389        axes=self.file.dimensions.keys()
390        return axes
391
392    def getSizeOfAxis(self,axis):
393        axisSize=self.file.dimensions[axis]
394        return axisSize
395
396    def getListofVariables(self):
397        variableList=self.file.variables.keys()
398
399        # Hack to test if removing climatology_bounds fixes pywms bug
400        if 'climatology_bounds' in variableList:
401            variableList.remove('climatology_bounds')
402
403        return variableList
404
405    def setAxis(self,axis):
406        self.axisobj=self.file.getAxis(axis)
407
408    def getAxisAttribute(self, att):
409        attValue=self.axisobj.attributes[att]
410        return attValue
411   
412    def getTimeUnits(self):
413        #this does the same as getAxisAttribute, but is a separate function as different formats handle time differently.
414        return self.getAxisAttribute('units')
415
416    def getDataForAxis(self):
417        data = self.axisobj.getValue()
418        return data
419
420    def setVariable(self,varname):
421        self.varobj=self.file.variables[varname]
422
423    def getVariableAxes(self):
424        varAxesList=self.varobj.getAxisIds()
425        return varAxesList
426
427    def getVariableAttribute(self,attName):
428        #NEED TO REWRITE THIS!
429        if attName == 'long_name':
430            try:
431                attribValue = self.varobj.long_name
432            except:
433                attribValue='no_long_name_found'
434        elif attName == 'units':
435            try:
436                attribValue = self.varobj.units
437            except:
438                attribValue='no_units_found'
439        elif attName == '_FillValue':
440            try:
441                attribValue=self.varobj._FillValue
442                attribValue=attribValue.toscalar()
443            except:
444                try:
445                    attribValue=self.varobj.missing_value
446                    attribValue=attribValue.toscalar()
447                except:
448                    attribValue = None 
449        return attribValue
450
451    def getDataForVar(self):
452        data = self.varobj.getValue()
453        print type(data)
454        return data
455
456    def getSubsetOfDataForVar(self, **kwargs):     
457        #delete any slicing indices and use names
458        try:
459            del kwargs['upper']
460        except:
461            pass
462        try:
463            del kwargs['lower']
464        except:
465            pass
466       
467        #takes keyword args defining subset eg
468        #subset=getSubsetOfDataForVar(latitude=(0.,10.0), longitude=(90, 100.0))
469        subset=None
470        lonkey='longitude'
471        if lonkey=='longitude' in kwargs.keys():   #this test needs to be much more robust...!
472            if  type(kwargs[lonkey]) is tuple:
473                if kwargs[lonkey][0] > kwargs[lonkey][1]:
474                    #subsetting greenwich meridian around 0
475                    lonMin = kwargs[lonkey][0]
476                    lonMax =kwargs[lonkey][1]
477                    kwargs[lonkey]=(0.0, lonMax)
478                    subset1=self.file(self.varobj.id,**kwargs)
479                    kwargs[lonkey]=(lonMin,359.9999)
480                    print 'kwargs: %s'%kwargs
481                    try:
482                        subset2=self.file(self.varobj.id,**kwargs)
483                        longitudeAxis=subset1.getAxisIndex(lonkey)
484                        #concatenate arrays along longitude             
485                        subset = cdms.MV.concatenate([subset1,subset2],axis=longitudeAxis)
486                    except:
487                        subset=subset1
488        if type(subset) is not cdms.tvariable.TransientVariable:
489            subset=self.file(self.varobj.id,**kwargs)
490        return subset
491
492    def getArraySizeOfVar(self):
493    #iterates through all dimensions in variable to get array size i.e a 3x3x3 array would have a size of 27
494        var = self.varobj
495        size = var.shape
496        varsize = 1
497        for item in size:
498            varsize = item *varsize
499        return varsize
500
501    def getShapeOfVar(self):
502        varShape = []
503        for item in self.varobj.shape:
504            varShape.append(item)
505        return varShape
506
507    def getLowLimits(self):
508        dimNames = self.varobj.getAxisIds()
509        lowlims = ""
510        for i in range (1, len(dimNames)):
511            #for now, assume low limit is always of form 1 1 1 ..
512            lowlims =lowlims + str(1)  +' '
513        return lowlims
514
515    def getHighLimits(self):
516        dimNames = self.varobj.getAxisIds()
517        highlims = ""
518        for i in range (1, len(dimNames)):
519            dimValue = self.file.dimensions[dimNames[i]]
520            highlims =highlims  + str(dimValue) +' '
521        return highlims
522       
523   
524class cdmlInterface(cdunifInterface):
525    #this is more  or less the cdunif interface but a few methods have been overwritten
526    def __init__(self):
527        #this all needs to be revisited in csml v2.
528        self.extractType='cdmlExtract'
529        self.extractPrefix = '_cdmlextract_'
530       
531    def getListOfAxes(self):
532        axes=self.file.axes.keys() 
533        return axes
534
535    def getSizeOfAxis(self,axis):
536        axisSize=self.file.axes[axis].length
537        return axisSize
538
539
540class RawFileInterface(AbstractDI):
541
542   def __init__(self):
543      self.extractType   = 'rawExtract'
544      self.extractPrefix = '_rawextract_'
545
546
547   def openFile(self, filename):
548      self.file = open(filename, "rb")
549
550
551   def closeFile(self):
552      self.file.close()
553
554
555   # Read the data from the raw file into a multidimensional Numeric array
556   def readFile(self, **kwargs):
557      # Valid values for numericType and their associated struct format
558      # character:
559      numericLookup = {
560         'uint8':'B', 'uint16':'H', 'uint32':'I', 'uint64':'Q',
561          'int8':'b',  'int16':'h',  'int32':'i',  'int64':'q',
562         'float':'f', 'double':'d'
563      }
564
565      # Valid values for endianness and their associated struct format
566      # character:
567      endiannessLookup = {
568         'big':'>',  'little':'<',  'host':'='
569      }
570
571      # Build the desired format string:
572      try:
573         format_str  = endiannessLookup[kwargs['endianness']]
574         format_str += str(reduce(lambda a,b:a*b, kwargs['dimensions']))
575         format_str += numericLookup[kwargs['numericType']]
576      except KeyError, TypeError:
577         raise TypeError("Missing or invalid parameters. Required: endianness=str, dimensions=list and numericType=str")
578
579      # Read and unpack data into a Numeric array.
580      self.data = Numeric.array(struct.unpack(format_str, self.file.read()))
581      self.data.shape = tuple(map(int, kwargs['dimensions']))
582     
583      # If numericTransform or fillValue were provided, store them as
584      # attributes.
585      if 'numericTransform' in kwargs:
586         self.numericTransform = NumericTransform.infixExpression(kwargs['numericTransform'])
587      if 'fillValue' in kwargs:
588         self.fillValue = kwargs['fillValue']
589
590
591   # Return the fill value, if set, and transform if necessary:
592   def getFillValue(self):
593      # Both fillValue and numericTransform attributes may or may
594      # not exist...
595      try:
596         return self.numericTransform.solve( n = float(self.fillValue) )
597      except AttributeError:
598         try:
599            return self.fillValue
600         except AttributeError:
601            return None
602
603
604   # This is a just a special case of getSubsetOfDataForVar. To avoid
605   # duplication of code, just subset the entire array. (getSubset.. is
606   # optimised for this case)
607   def getDataForVar(self):
608      return self.getSubsetOfDataForVar(lower = (0,)*len(self.data.shape),
609                                        upper = self.data.shape)
610
611
612   # Subset the n-dimensional file based on array indices. Accepts parameters
613   # in the form of e.g.
614   #
615   # getSubsetOfDataForVar( lower=(0,0), upper=(1,2) )
616   #
617   # Note: The rank of the required subset, must exactly match the
618   # rank of the original data: len(lower) == len(upper) == rank of file
619   #
620   def getSubsetOfDataForVar(self, **kwargs):
621      # Assume subset parameters are passed as: lower=(0,0) upper=(512,512)
622      if 'lower' not in kwargs or 'upper' not in kwargs:
623         # Have not specified any subset parameters that we recognise, so raise
624         # an exception:
625         raise NotImplementedError("Only supports subsetting with lower+upper array indices")
626      elif not len(kwargs['lower']) == len(kwargs['upper']) == len(self.data.shape):
627         # Rank of subset is not the same as rank of full data array. so raise
628         # an exception:
629         raise NotImplementedError("Only supports subsets of same rank as full dataset")
630      elif Numeric.sometrue(Numeric.greater(kwargs['upper'], self.data.shape)):
631         # Requested upper bound of subset is beyond the size of the the full
632         # data array, so raise an exception
633         raise IndexException("Subset out of range")
634      elif Numeric.sometrue(Numeric.less( kwargs['upper'], Numeric.zeros(len(self.data.shape)))):
635         # Requested lower bound of subset is beyond the size of the the full
636         # data array, so raise an exception
637         raise IndexException("Subset out of range")
638      elif Numeric.sometrue(Numeric.less_equal(kwargs['upper'], kwargs['lower'])):
639         # lower bound <= upper_bound for at least one dimension, so raise an
640         # exception
641         raise IndexException("Upper bound less than lower bound")
642      elif tuple(kwargs['lower']) == (0,)*len(self.data.shape) and tuple(kwargs['upper']) == self.data.shape:
643         # Special case of requested subset == entire data file.
644         subset = self.data
645      else:
646         # We are okay to subset.
647
648         # I cant see any nice (and speedy) way of subsetting a Numeric
649         # array of arbitrary rank without the use of eval. By doing it
650         # this way, we can make use of the (possible) acceleration provided
651         # by Numeric/NumPy.
652         slices = []
653         for i in range(len(self.data.shape)):
654            lower = int(kwargs['lower'][i])
655            upper = int(kwargs['upper'][i]) 
656            slices.append(str(lower)+':'+str(upper))
657         subset = eval("self.data["+','.join(slices)+"]")
658
659      # Attempt to perform the numericTransform on the data array, if we get
660      # AttributeError, it most likely means that numericTransform was not
661      # specified, so return the data as-is
662      try:
663         return self.numericTransform.solve( n = subset )
664      except AttributeError:
665         return subset.copy()
666
667
668# Interface for reading data from image files. Requires PIL Image module.
669class ImageFileInterface(RawFileInterface):
670   def __init__(self):
671      self.extractType   = 'imageExtract'
672      self.extractPrefix = '_imageextract_'
673
674   def openFile(self, filename):
675      self.file = Image.open(filename)
676
677   def closeFile(self):
678      self.file = None #...Image does not seem to have a close() method.
679
680   def readFile(self, **kwargs):
681      # Convert the image to a Numeric array
682      self.data = Numeric.array(self.file.getdata())
683
684      # The order of dimensions in Numeric arrays do not correspond with
685      # the order of dimensions in traditional image files, so we need to
686      # reverse the dimensions before shaping the array.
687      dimensions = map(int,self.file.size)
688      dimensions.reverse()
689      self.data.shape = tuple(dimensions) # self.data.shape = (Y,X)
690
691      if 'numericTransform' in kwargs:
692         # numericTransform was specified, so compile the expression:
693         self.numericTransform = NumericTransform.infixExpression(kwargs['numericTransform'])
694      if 'fillValue' in kwargs:
695         self.fillValue = kwargs['fillValue']
Note: See TracBrowser for help on using the repository browser.