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

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

added code to assist in preservation of CF (and other) attributes. Not completely working atm

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1#!/usr/bin/env python
2# Adapted for numpy/ma/cdms2 by convertcdms.py
3
4'''
5**************************************************************************************
6csmldataIface.py
7contains classes for interfacing with various files
8currently supports cdunif (NetCDF, PP, Grib(untested)) And Nappy (NASAAmes)
9use by instantiating the factory class: DataInterface
10last updated 28 January 2008
11Dominic Lowe, BADC
12**************************************************************************************
13'''
14
15#caching module:
16from open_file_cache import OpenFileCache
17
18import pdb
19import cdms2 as cdms 
20import numpy.oldnumeric as Numeric
21
22try:
23    import nappy 
24except ImportError:
25    print "NASAAmes interface not available. CSML will still work, but won't support NASA Ames files"
26import string
27import sys
28import csml.csmllibs.csmltime
29# This is required to prevent Numeric arrays being truncated when printed.
30import numpy.oldnumeric.ma as MA
31#MA.set_print_limit(0)
32import numpy; numpy.set_printoptions(threshold=0)
33
34# Generic mathematical expression solver, required by the raw and image
35# interfaces
36import numpy.oldnumeric as NumericTransform
37
38try:
39    # Part of the PIL. Required for ImageFileInterface:
40    import Image
41except ImportError:
42    print "Python Image Library not available. CSML will still work, but won't support image files"
43
44try:
45    # Required for RawFileInterface:
46    import struct
47except ImportError:
48    print 'could not import struct'
49    pass
50
51cache=OpenFileCache()
52cache.limit=40
53
54class DataInterface(object):
55        '''DataInterface is the entry point for all data types. Use DataInterface and setInterfaceType to instantiate the correct
56        subclass for data''' 
57        def __init__(self):
58                self.iface ='None'
59               
60        def setInterfaceType(self,interfaceType):
61                '''returns approprate data interface
62                @param interfaceType:     should correspond to available datainterface subclasses
63                @return:    a subclass of DataInterface e.g. a ImageFileInterface'''
64                self.iface =interfaceType
65                if self.iface == 'nappy':
66                        return NappyInterface()
67                elif self.iface == 'cdunif':
68                        return cdunifInterface()
69                elif self.iface == 'raw':
70                        return RawFileInterface()
71                elif self.iface == 'image':
72                        return ImageFileInterface()
73               
74               
75        def getUnknownInterfaceType(self, filename):
76                '''if the interface type is not known at the time of instantiation, then use
77                this method to examine the file and return the correct interface (if it exists).
78                @param filename:      name of data file'''
79                fileExtension = str(filename)[-3:]
80                if fileExtension == '.nc':
81                        return cdunifInterface()
82                if fileExtension == 'qxf':
83                        return cdunifInterface()
84                elif fileExtension == '.pp':
85                        return cdunifInterface()
86                elif fileExtension == 'ctl':
87                        return cdunifInterface()
88                elif fileExtension == 'xml':
89                        return cdmlInterface()
90                elif fileExtensions == 'png' or fileExtension == 'gif':
91                        return ImageFileInterface()
92                else:
93                        try:
94                                nappy.readFFI(filename) in [1001,1010,1020,2010,2110,2160,2310,3010,4010]
95                                return NappyInterface()
96                        except:
97                                print "Could not establish file type"
98                                print "See csmldataiface.py"
99                       
100       
101                                                       
102class AbstractDI(object):               
103        '''Abstract data interface class
104        does nothing, but contains templates for methods required for a data interface class
105        individual interfaces (e.g NappyInterface) should override these methods. This is the minimum
106        set required for data access and subsetting.
107        Note that if scanning of data (for CSML creation) is required more methods may be needed - see cdunif interface'''
108
109        def __init__(self):
110                ''' When subclassed: unique 'type'  and 'prefix' must be defined'''
111                self.extractType=''
112                self.extractPrefix = ''
113                                       
114        def openFile(self, filename):
115           ''' When subclassed: opens file, sets self.file'''
116           raise NotImplementedError 
117               
118        def closeFile(self):
119           '''closes file, probably needs to be overwritten by subclass'''
120           try:
121               self.file.close()
122           except:
123               raise NotImplementedError
124         
125        def setAxis(self,axis):
126               '''set' the name of the current axis , must be overwritten by subclass
127           this may just involve a stub (see NASAAmes interface) or may involve
128           calling a real set method of the underlying api (see cdunif Interface)'''
129               raise NotImplementedError 
130                       
131        def getDataForAxis(self):
132               ''' return all data for axis, must be overwritten by subclass '''
133               raise NotImplementedError       
134       
135        def setVariable(self,varname):
136             '''As for setAxis, 'set' the name of the current axis , must be overwritten by subclass'''
137             raise NotImplementedError
138       
139        def getDataForVar(self):
140               '''return all data for variable, must be overwritten by subclass'''
141               raise NotImplementedError       
142       
143        def getSubsetOfDataForVar(self,**kwargs):
144            '''return subset of data for variable, must be overwritten by subclass'''
145            raise NotImplementedError
146       
147       
148class NappyInterface(AbstractDI):       
149        ''' Data Interface for Nappy (NASA Ames Processing in Python)'''
150   
151        def __init__(self):
152                self.extractType='NASAAmesExtract'
153                self.extractPrefix = '_naextract_'
154                                         
155        def openFile(self, filename):
156                ''' open file'''
157                self.file=nappy.openNAFile(filename)
158
159        def getListOfAxes(self):
160                ''' gets list of axes in file'''
161                axes=self.file.XNAME
162                #print 'before units stripped' + str(axes)
163                axes=self.stripunits(axes)
164                #print 'after units stripped' + str(axes)
165                return axes
166
167        def setAxis(self,axis):
168                '''sets selected axis in Data Interface'''
169                axes = self.getListOfAxes()
170                self.axisstub=axes.index(axis)
171
172        def getAxisAttribute(self, att):
173                ''' not implemented yet -should return axis attribute'''
174                #need to do something here...? maybe
175                attValue=None
176                return attValue
177
178        def getTimeUnits(self):
179                ''' returns units of time axis'''
180                axes = self.getListOfAxes()
181                for axis in axes:
182                        if string.find(string.upper(axis),'SECONDS SINCE') != -1:
183                                #found possible time axis.
184                                if axis[-3:]=='UTC':
185                                    units =string.lower(axis[:-4]) #hack!
186                                    units=units.replace('/','-') #try and clean it up
187                                else:
188                                    units=string.lower(axis)
189                                break
190                        elif string.find(string.upper(axis),'HOURS SINCE') != -1:
191                                #found possible time axis.
192                                units =(str(axis))
193                                break
194                        elif string.find(string.upper(axis),'DAYS SINCE') != -1:
195                                #found possible time axis.
196                                units =(str(axis))
197                                break
198                       
199                #revisit with udunits python library?
200                return units
201
202
203        def getDataForAxis(self):           
204            ''' gets data for an axis'''
205            if self.file.X == None:
206                    #print 'reading data....'
207                    self.file.readData()
208            for x in self.file.X:
209                pass
210                #print type(x)
211           
212           
213            d=Numeric.array(self.file.X)
214           
215            if type(self.file.X[1])==list:
216            #if len(self.file.X) > 0:
217                    data = self.file.X[self.axisstub]
218            else:
219                    data =self.file.X
220            #print data
221            data=Numeric.array(data)
222            return data
223
224        def getSizeOfAxis(self,axis):
225                ''' gets size of an axis'''
226                #check this function is okay
227                #is time always the first dimension in NA??
228                axes = self.getListOfAxes()
229                axisPosition=axes.index(axis)
230                #print "axis position" + str( axisPosition)
231                #print "NX" + str(self.file.NX)
232                try :
233                        axisSize=self.file.NX[axisPosition-1]
234                except:
235                        axisSize ='Unknown axis size'
236                return axisSize
237
238        def getListofVariables(self):
239                '''gets list of variables in the file'''
240                variableList=self.stripunits(self.file.VNAME)
241                return variableList
242
243        def setVariable(self,varname):
244                ''' sets the specified variable to be the selected variable'''
245                vlist=self.getListofVariables()
246                self.varstub=vlist.index(varname)
247
248        def getVariableAxes(self):
249                ''' Get axis names for a  variable '''
250                #hmm, now with Nasa Ames the Axis will be the same for all variables.
251                #so just use the getListOfAxes function again
252                #I think... check this!
253                varAxesList=self.getListOfAxes()
254                return varAxesList
255
256        def getVariableAttribute(self,attName):
257                ''' get an attribute for a variable'''
258                if attName =='units':
259                        #strip the units (attribute) from the variable
260                        unitslist=self.getUnits(self.file.VNAME)
261                        attribValue = unitslist[self.varstub]
262                        try:
263                                attribValue = unitslist[self.varstub]
264                        except:
265                                attribValue = 'unknown'
266                else:
267                        attribValue = 'unknown'
268                return attribValue
269
270        def getDataForVar(self):
271            ''' get data for a  variable'''
272            #NOTE TO SELF:
273            #Review this function (and in fact all of nasa ames data interface...)
274            if self.file.V == None:
275                    #print 'reading data....'
276                    self.file.readData()
277
278            try:
279                if type(self.file.V[1])==list:
280                    data = self.file.V[self.varstub]
281            #else:
282            #   data =self.file.X
283            #   print data
284                return data
285            except:
286                data = self.file.X
287                # print data
288                return data
289
290        def getArraySizeOfVar(self):
291            '''iterates through all dimensions in variable to get array size i.e a 3x3x3 array would have a size of 27'''
292
293            dimlist=self.file.NX
294            varsize =1
295            for item in dimlist:
296                    varsize = varsize * item
297                    #print "VARSISZE" + str(varsize)
298            return varsize
299
300        def getShapeOfVar(self):
301            '''this returns a list with the shape of the variable'''
302            varShape = []
303            for item in self.file.NX:
304                varShape.append(item)
305            return varShape
306
307        def getLowLimits(self):
308            ''' gets low limits of the data'''
309            lowlims = ""
310            for i in range (0, len(self.file.NX)):
311                    #for now, assume low limit is always of form 1 1 1 ..
312                    lowlims =lowlims + str(1)  +' '
313            return lowlims
314
315        def getHighLimits(self):
316            ''' gets high limits of the data'''
317            highlims = ""
318            for i in range (0, len(self.file.NX)):
319                    dimValue = self.file.NX[i]
320                    highlims =highlims  + str(dimValue) +' '
321            return highlims
322
323
324        def stripunits(self,listtostrip):
325            ''' strips units of measure from list
326            eg ['Universal time (hours)', 'Altitude (km)', 'Latitude (degrees)', 'Longitude (degrees)']
327            becomes ['Universal time', 'Altitude', 'Latitude', 'Longitude']'''
328            cleanlist = []
329            for item in listtostrip:
330                    openbracket=string.find(item,'(')
331                    if openbracket != -1:
332                            #if brackets exist, strip units.
333                            item=item[:openbracket-1]
334                    cleanlist.append(item)
335            return cleanlist
336
337        def getUnits(self,listwithunits):
338            '''gets units from list
339            eg ['Universal time (hours)', 'Altitude (km)', 'Latitude (degrees)', 'Longitude (degrees)']
340            becomes ['hours', 'km', 'degrees', 'degrees']'''
341            unitlist=[]
342            for item in listwithunits:
343                    openbracket=string.find(item,'(')
344                    item = item[openbracket+1:-1]
345                    unitlist.append(item)
346            return unitlist
347
348        def getTimes(self):
349            '''This function attempts to determine the time axis and read the time data
350            it may well not manage it.'''
351            axes = self.getListOfAxes()
352            for axis in axes:
353                if string.find(string.upper(axis),'TIME') != -1:
354                        #found possible time axis.
355                        self.setAxis(axis)
356                        times=self.getDataForAxis()
357                        break
358                elif string.find(string.upper(axis),'SECONDS SINCE') != -1:
359                        #found possible time axis.
360                        self.setAxis(axis)
361                        times=self.getDataForAxis()
362                        break
363                elif string.find(string.upper(axis),'HOURS SINCE') != -1:
364                        #found possible time axis.
365                        self.setAxis(axis)
366                        times=self.getDataForAxis()
367                        break
368                elif string.find(string.upper(axis),'DAYS SINCE') != -1:
369                        #found possible time axis.
370                        self.setAxis(axis)
371                        times=self.getDataForAxis()
372                        break
373            times=Numeric.array(times)
374            return times
375
376
377
378class cdunifInterface(AbstractDI):
379    '''Data Interface for cdunif (netcdf & pp formats & grib (not tested with grib)'''
380
381    def __init__(self):
382        '''these are just temporary values until we can determine whether the
383        file is netcdf pp or grib'''
384        self.extractType='cdunifExtract'
385        self.extractPrefix = '_cdunifextract_'
386   
387    def openFile(self, filename):
388        ''' open file'''
389        self._filename=filename
390        self.file=cache.openFile(cdms.open,filename)
391        #now we have the file name can properly determine extractType/Prefix
392        fileExtension = str(filename)[-3:]
393        if fileExtension == '.nc':
394            self.extractType = 'NetCDFExtract'
395            self.extractPrefix = '_ncextract_'
396        elif fileExtension == '.qxf':
397            self.extractType = 'NetCDFExtract'
398            self.extractPrefix = '_ncextract_'
399        elif fileExtension == '.pp':
400            self.extractType  = 'PPExtract'
401            self.extractPrefix = '_ppextract_'
402        elif fileExtension == 'ctl':
403            self.extractType = 'GRIBExtract'
404            self.extractPrefix = '_gribextract_'
405        elif fileExtension == 'xml': 
406            self.extractType = 'NetCDFExtract'  #okay this isn't true, but ok for testing
407            self.extractPrefix = '_ncextract__' 
408   
409    def closeFile(self):
410        pass #let the cache manage this
411   
412    def getListOfAxes(self):
413        ''' returns list of dimensions/axes'''
414        axes=self.file.dimensions.keys()
415        return axes
416
417    def getSizeOfAxis(self,axis):
418        ''' returns size of axis'''
419        axisSize=self.file.dimensions[axis]
420        return axisSize
421
422    def getListofVariables(self):
423        ''' returns list of variables'''
424        variableList=self.file.variables.keys()
425
426        # Hack to test if removing climatology_bounds fixes pywms bug
427        if 'climatology_bounds' in variableList:
428            variableList.remove('climatology_bounds')
429
430        return variableList
431
432    def setAxis(self,axis):
433        ''' set axis'''
434        if not hasattr(self, 'file'):
435            raise Exception, 'Could not open/find underlying file: %s  -  If you are the system maintainer check file paths and permissions'%self._filename
436        self.axisobj=self.file.getAxis(axis)
437
438           
439
440    def getAxisAttribute(self, att):
441        ''' get attribute of set axis'''
442        attValue=self.axisobj.attributes[att]
443        return attValue
444   
445    def getTimeUnits(self):
446        ''' get time units - only works if set axis is time...'''
447        #this does the same as getAxisAttribute, but is a separate function as different formats handle time differently.
448        return self.getAxisAttribute('units')
449
450    def getDataForAxis(self):
451        ''' get data for entire axis'''
452        data = self.axisobj.getValue()
453        return data
454
455    def setVariable(self,varname):
456        ''' sets variable '''
457        self.varobj=self.file.variables[varname]
458
459    def getVariableAxes(self):
460        ''' gets axes for set variable'''
461        varAxesList=self.varobj.getAxisIds()
462        return varAxesList
463   
464    def getVariableAttributes(self):
465        '''gets names of all attributes belonging to variable'
466        returns list'''
467        return self.varobj.listattributes()
468
469    def getVariableAttribute(self,attName):
470        ''' gets requested attribute for set variable'''
471        if attName == '_FillValue':
472            try:
473                attribValue=self.varobj._FillValue
474                attribValue=attribValue.item()
475            except:
476                try:
477                    attribValue=self.varobj.missing_value
478                    attribValue=attribValue.item()
479                except:
480                    attribValue = None 
481        else:
482            attribValue = self.varobj.getattribute(attName)
483        return attribValue
484
485    def getDataForVar(self):
486        ''' get data for set variable'''
487        data = self.varobj.getValue()
488        return data
489   
490    def _fixLongitudeRequest(self, **kwargs):
491        ''' fixes problem of longitude requests taking different forms'''
492        lonkey='longitude'
493        if lonkey in kwargs.keys():   #this test needs to be much more robust...!
494            if  type(kwargs[lonkey]) is tuple:
495                newkws=[]
496                for val in kwargs[lonkey]:
497                    newval=((val + 180) % 360) - 180
498                    newkws.append(newval)
499               
500                 #lon = ((t_lon + 180) % 360) - 180
501                #lonmin=kwargs[lonkey][0] - (kwargs[lonkey][0]/180)*360.0
502                #lonmax=kwargs[lonkey][1] - (kwargs[lonkey][1]/180)*360.0
503               
504                #if newkws[0] > newkws[1]:
505                    #flip=[newkws[1], newkws[0]]
506                    #newkws=flip
507                kwargs[lonkey]=tuple(newkws)
508        return kwargs
509               
510    def _componentTimes(self, times):
511        ''' handles time formats'''
512        if type(times) is tuple:
513            comptimes=[]
514            for t in times:
515                compt=csml.csmllibs.csmltime.getCDtime(t)
516                comptimes.append(compt)
517            result=tuple(comptimes)
518        else: #just a single time
519            result=csml.csmllibs.csmltime.getCDtime(times)
520        return result
521       
522   
523    def getSubsetOfDataForVar(self, **kwargs):     
524        ''' Main subsetting method. Kwargs contains dictionary of requested subset.'''
525        #put any slicing indices aside for later and use names       
526        try:
527            upper=kwargs['upper']
528            del kwargs['upper']
529        except:
530            pass
531        try:
532            lower=kwargs['lower']
533            del kwargs['lower']
534        except:
535            pass
536        if 'time' in kwargs:
537            kwargs['time']=self._componentTimes(kwargs['time'])
538        #kwargs=self._fixLongitudeRequest(**kwargs)
539        #return self.file(self.varobj.id,**kwargs)
540        try:
541            #takes keyword args defining subset eg
542            #subset=getSubsetOfDataForVar(latitude=(0.,10.0), longitude=(90, 100.0))
543            subset=None
544            lonkey=None
545            for possiblelonkey in kwargs.keys():
546                if possiblelonkey in ['longitude', 'lon', 'longitude0', 'longitude1']:
547                    lonkey=possiblelonkey
548            if lonkey:   #this test needs to be much more robust... the problem is - how to establish which is the longitude axis without axes to the csml doc
549                if  type(kwargs[lonkey]) is tuple:
550                    if kwargs[lonkey][0] > kwargs[lonkey][1]:
551                        #subsetting greenwich meridian around 0
552                        lonMin = kwargs[lonkey][0]
553                        lonMax =kwargs[lonkey][1]
554                        kwargs[lonkey]=(0.1, lonMax)
555                        subset1=self.file(self.varobj.id,**kwargs)
556                        kwargs[lonkey]=(lonMin,360)
557                        try:
558                            subset2=self.file(self.varobj.id,**kwargs)
559                            longitudeAxis=subset1.getAxisIndex(lonkey)
560                            #concatenate arrays along longitude             
561                            subset = cdms.MV2.concatenate([subset2,subset1],axis=longitudeAxis)
562                        except:
563                            subset=subset1
564            if type(subset) is not cdms.tvariable.TransientVariable:
565                subset=self.file(self.varobj.id,**kwargs)
566        except:             
567                return None
568                #try and slice with indices instead
569        return subset
570   
571   
572    def getArraySizeOfVar(self):
573        '''iterates through all dimensions in variable to get array size i.e a 3x3x3 array would have a size of 27'''
574        var = self.varobj
575        size = var.shape
576        varsize = 1
577        for item in size:
578            varsize = item *varsize
579        return varsize
580
581    def getShapeOfVar(self):
582        ''' get shape of set variable '''
583        varShape = []
584        for item in self.varobj.shape:
585            varShape.append(item)
586        return varShape
587
588    def getLowLimits(self):
589        ''' get gml low limits of set variable'''
590        dimNames = self.varobj.getAxisIds()
591        lowlims = ""
592        for i in range (1, len(dimNames)):
593            #for now, assume low limit is always of form 1 1 1 ..
594            lowlims =lowlims + str(1)  +' '
595        return lowlims
596
597    def getHighLimits(self):
598        ''' get gml high limits of set variable'''
599        dimNames = self.varobj.getAxisIds()
600        highlims = ""
601        for i in range (1, len(dimNames)):
602            dimValue = self.file.dimensions[dimNames[i]]
603            highlims =highlims  + str(dimValue) +' '
604        return highlims
605       
606   
607class cdmlInterface(cdunifInterface):
608    ''' this is more  or less the cdunif interface but a few methods have been overwritten to read CDML'''
609    def __init__(self):
610        #this all needs to be revisited in csml v2.
611        self.extractType='cdmlExtract'
612        self.extractPrefix = '_cdmlextract_'
613       
614    def getListOfAxes(self):
615        ''' get list of axis in file'''
616        axes=self.file.axes.keys() 
617        return axes
618
619    def getSizeOfAxis(self,axis):
620        ''' get size of set axis'''
621        axisSize=self.file.axes[axis].length
622        return axisSize
623
624
625class RawFileInterface(AbstractDI):
626    ''' interface for raw file types'''
627    def __init__(self):
628        self.extractType   = 'rawExtract'
629        self.extractPrefix = '_rawextract_'
630 
631 
632    def openFile(self, filename):
633        self.file = open(filename, "rb")
634
635    def closeFile(self):
636            self.file.close()
637
638
639    # Read the data from the raw file into a multidimensional Numeric array
640    def readFile(self, **kwargs):
641        # Determine the numeric type:
642        if 'numericType' in kwargs:
643            try:
644                numericTypeCode = {
645                    'uint8':Numeric.UInt8,
646                    'uint16':Numeric.UInt16,
647                    'uint32':Numeric.UInt32,
648                    'int8':Numeric.Int8,
649                    'int16':Numeric.Int16,
650                    'int32':Numeric.Int32,
651                    'float':Numeric.Float32,
652                    'double':Numeric.Float64
653                }[kwargs['numericType']]
654            except KeyError:
655                raise TypeError("Invalid numericType: " + str(kwargs['numericType']))
656        else:
657            raise KeyError("numericType is mandatory.")
658       
659        # Read the file into a numpy array:
660        self.data = Numeric.fromstring(self.file.read(), numericTypeCode)
661        # If necessary, swap the byte order:
662        if 'endianess' in kwargs:
663            if ((kwargs['endianness'] == 'little' and not Numeric.LittleEndian) or (kwargs['endianness'] == 'big' and Numeric.LittleEndian)):
664                self.data = self.data.byteswap()   
665        # Declare the shape of the array:
666        dimensions = map(int,kwargs['dimensions'])
667        dimensions.reverse()
668        self.data.shape = tuple(dimensions)
669        # If numericTransform or fillValue were provided, store them as
670        # attributes.
671        if 'numericTransform' in kwargs:
672            self.numericTransform = NumericTransform.infixExpression(kwargs['numericTransform'])
673        if 'fillValue' in kwargs:
674            self.fillValue = kwargs['fillValue']
675   
676   # Return the fill value, if set, and transform if necessary:
677    def getFillValue(self):
678      # Both fillValue and numericTransform attributes may or may
679      # not exist...
680        try:
681            return self.numericTransform.solve( n = float(self.fillValue) )
682        except AttributeError:
683            try:
684                return self.fillValue
685            except AttributeError:
686                return None
687
688
689   # This is a just a special case of getSubsetOfDataForVar. To avoid
690   # duplication of code, just subset the entire array. (getSubset.. is
691   # optimised for this case)
692    def getDataForVar(self):
693        return self.getSubsetOfDataForVar(lower = (0,)*len(self.data.shape),
694                                        upper = self.data.shape)
695
696
697   # Subset the n-dimensional file based on array indices. Accepts parameters
698   # in the form of e.g.
699   #
700   # getSubsetOfDataForVar( lower=(0,0), upper=(1,2) )
701   #
702   # Note: The rank of the required subset, must exactly match the
703   # rank of the original data: len(lower) == len(upper) == rank of file
704   #
705    def getSubsetOfDataForVar(self, **kwargs):
706        # Assume subset parameters are passed as: lower=(0,0) upper=(512,512)
707        if 'lower' not in kwargs or 'upper' not in kwargs:
708            # Have not specified any subset parameters that we recognise, so raise
709            # an exception:
710            raise NotImplementedError("Only supports subsetting with lower+upper array indices")
711        elif not len(kwargs['lower']) == len(kwargs['upper']) == len(self.data.shape):
712            # Rank of subset is not the same as rank of full data array. so raise
713            # an exception:
714            raise NotImplementedError("Only supports subsets of same rank as full dataset")
715        elif Numeric.sometrue(Numeric.greater(kwargs['upper'], self.data.shape)):
716            # Requested upper bound of subset is beyond the size of the the full
717            # data array, so raise an exception
718            raise IndexError("Subset out of range")
719        elif Numeric.sometrue(Numeric.less( kwargs['upper'], Numeric.zeros(len(self.data.shape)))):
720            # Requested lower bound of subset is beyond the size of the the full
721            # data array, so raise an exception
722            raise IndexError("Subset out of range")
723        elif Numeric.sometrue(Numeric.less_equal(kwargs['upper'], kwargs['lower'])):
724            # lower bound <= upper_bound for at least one dimension, so raise an
725            # exception
726            raise IndexError("Upper bound less than lower bound")
727        elif tuple(kwargs['lower']) == (0,)*len(self.data.shape) and tuple(kwargs['upper']) == self.data.shape:
728            # Special case of requested subset == entire data file.
729            subset = self.data
730        else:
731            # We are okay to subset.
732   
733            # I cant see any nice (and speedy) way of subsetting a Numeric
734            # array of arbitrary rank without the use of eval. By doing it
735            # this way, we can make use of the (possible) acceleration provided
736            # by Numeric/NumPy.
737            slices = []
738            for i in range(len(self.data.shape)):
739                lower = int(kwargs['lower'][i])
740                upper = int(kwargs['upper'][i]) 
741                slices.append(str(lower)+':'+str(upper))
742            subset = eval("self.data["+','.join(slices)+"]")
743   
744        # Attempt to perform the numericTransform on the data array, if we get
745        # AttributeError, it most likely means that numericTransform was not
746        # specified, so return the data as-is
747        try:
748            return self.numericTransform.solve( n = subset )
749        except AttributeError:
750            return subset.copy()
751
752class ImageFileInterface(RawFileInterface):
753    '''Interface for reading data from image files. Requires PIL Image module.'''
754    def __init__(self):
755        self.extractType   = 'imageExtract'
756        self.extractPrefix = '_imageextract_'
757   
758    def image2array(self,im):
759        #Adapted from code by Fredrik Lundh, http://www.pythonware.com
760        #http://effbot.org/zone/pil-numpy.htm
761            if im.mode not in ("L", "F"):
762                raise ValueError, "can only convert single-layer images"
763            if im.mode == "L":
764                a = Numeric.fromstring(im.tostring(), Numeric.UnsignedInt8)
765            else:
766                a = Numeric.fromstring(im.tostring(), Numeric.Float32)
767            a.shape = im.size[1], im.size[0]
768            return a
769
770    def openFile(self, filename):
771        self.file = Image.open(filename)
772
773    def closeFile(self):
774        self.file = None #...Image does not seem to have a close() method.
775
776    def readFile(self, **kwargs):
777        # Convert the image to a Numeric array
778     
779        self.data=self.image2array(self.file)
780        #slower method:
781        #self.data = Numeric.array(self.file.getdata())
782
783        if 'numericTransform' in kwargs:
784            # numericTransform was specified, so compile the expression:
785            self.numericTransform = NumericTransform.infixExpression(kwargs['numericTransform'])
786        if 'fillValue' in kwargs:
787            self.fillValue = kwargs['fillValue']
788
789class HDF5Interface(AbstractDI):
790    '''Data Interface for HDF5'''
791    def __init__(self):
792        self.extractType='hdfExtract'
793        self.extractPrefix='_hdfextract_'
794   
795    def openFile(self, filename):
796        '''some code to open the file'''
797        pass
798
799    def setAxis(self,axis):
800        '''some code to set an axis to be queried, may not need to do much, depending on your format'''
801        pass
802       
803    def getDataForAxis(self):
804        '''some code to return the values for an axis'''
805        return data
806
807    def setVariable(self,varname):
808        '''some code to set a variable to be queried, may not need to do much, depending on your format'''
809        pass
810
811    def getDataForVar(self):
812        '''some code to return all values for a variable'''
813        return data
814
815    def getSubsetOfDataForVar(self, **kwargs):
816        '''takes keyword args defining subset eg
817        subset=getSubsetOfDataForVar(latitude=(0.,10.0), longitude=(90, 100.0), ...)
818        and returns a subset of data for tha variable '''
819        return data
820
821    def closeFile(self):
822        '''some code to close the file'''
823        pass
Note: See TracBrowser for help on using the repository browser.