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

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

initial version of PML methods (will need testing) - mggr checking in for mhen

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