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

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

(mggr committing for mhen): 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   def readFile(self, **kwargs):
538      # index into an array that determines the type based on various
539      # combinations of settings in the CSML.  This is very clever, but
540      # rather obscure.. blame mhen :p
541      format_index = 0
542
543      if kwargs['type'] == 'float':
544         format_index |= 12 # Float: Set bit 2 (float) and bit 3 (signed)
545      elif kwargs['signedness'] != 'unsigned':
546         format_index |= 8  # Signed Int: set bit 3 (signed)
547         
548      # Set bits 0 and 1 depending on whether the data are 8, 16, 32 or 64
549      # bits wide.
550      if kwargs['depth'] in [8,16,32,64]:
551         format_index |= [8,16,32,64].index(kwargs['depth'])
552         
553      # Set the endianness of the format string. (Set to Host-Endian if not
554      # defined
555      if kwargs['endianness'] == 'big':
556         format_str = '>'
557      elif kwargs['endianness'] == 'little':
558         format_str = '<'
559      else:
560         format_str = '='
561         
562      # Set the size of the data to read (product of dimensions)
563      format_str += str(reduce(lambda a,b:a*b, kwargs['dimensions']))
564     
565      # Lookup the type character
566      format_str += "BHIQ----bhiq--fd"[format_index]
567
568      self.data = Numeric.array(struct.unpack(format_str, self.file.read()))
569      self.data.shape = tuple(map(int, kwargs['dimensions']))
570     
571      if 'numericTransform' in kwargs:
572         # if a numericTransform was supplied, then compile it:
573         self.numericTransform = NumericTransform.infixExpression(kwargs['numericTransform'])
574
575      if 'fillValue' in kwargs:
576         self.fillValue = kwargs['fillValue']
577
578
579
580   # Allows data to be indexed as if it were a multidimensional array:
581   def index(self, *indices):
582      index1d = 0
583      factor  = 1
584      for n in range(len(indices)):
585         index1d += factor * indices[n]
586         factor  *= self.dimensions[n]
587      return self.data[index1d]
588
589
590   # Return the fill value, if set, and transform if necessary:
591   def getFillValue(self):
592      if(self.fillValue):
593         if(self.numericTransform):
594            return self.numericTransform.solve( n = float(self.fillValue) )
595         else:
596            return self.fillValue
597      else:
598         return None
599
600
601   def getDataForVar(self):
602      if self.numericTransform:
603         # Solve the numeric transform and return the resulting array of floats
604         return self.numericTransform.solve( n = self.data )
605      else:
606         return self.data
607
608
609   def getSubsetOfDataForVar(self, **kwargs):
610      # Assume subset parameters are passed as: lower=(0,0) upper=(512,512)
611      if 'lower' not in kwargs or 'upper' not in kwargs:
612         raise NotImplementedError("Only supports subsetting with lower+upper array indices")
613      elif len(kwargs['lower']) != len(self.data.shape) or len(kwargs['upper']) != len(self.data.shape):
614         raise NotImplementedError("Only supports subsets of same dimensionality as full dataset")
615      else:
616         if not self.data:
617            self.readFile()
618         subset = []
619
620         # Recursive function for the purpose of iterating over an array of
621         # arbitrary dimensionality. (Iterate over the desired subset, and
622         # append each data item encountered to the subset array)
623         ptr = [None] * len(self.data.shape)
624         def iterate(n = 0):
625            for ptr[n] in range(kwargs['lower'][n], kwargs['upper'][n] + 1):
626               if n < len(self.data.shape) - 1:
627                  iterate(n + 1)
628               else:
629                  subset.insert(0, self.index(*ptr))
630         iterate()
631
632      if self.numericTransform:
633         # Solve the numeric transform and return the resulting array of floats
634         return  self.numericTransform.solve( n = subset )
635      else:
636         return subset
637
638
639# Interface from reading data from image files. Requires PIL Image module
640class ImageFileInterface(RawFileInterface):
641   def __init__(self):
642      self.extractType   = 'imageExtract'
643      self.extractPrefix = '_imageextract_'
644
645   def openFile(self, filename):
646      self.file = Image.open(filename)
647
648   def closeFile(self):
649      self.file = None
650
651   def readFile(self, **kwargs):
652      # Convert the image to a Numeric array
653      self.data = Numeric.array(im.getdata())
654      self.data.shape = tuple(map(int, im.size))
655
656      if 'numericTransform' in kwargs:
657         # numericTransform was specified, so compile the expression:
658         self.numericTransform = NumericTransform.infixExpression(kwargs['numericTransform'])
659      if 'fillValue' in kwargs:
660         self.fillValue = kwargs['fillValue']
Note: See TracBrowser for help on using the repository browser.