source: TI02-CSML/branches/csml-pml/csmllibs/csmldataiface.py @ 2461

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/TI02-CSML/branches/csml-pml/csmllibs/csmldataiface.py@2461
Revision 2461, 22.8 KB checked in by mggr, 14 years ago (diff)

(mggr checking in for mhen and hoping not to mess up): ongoing development

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