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

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

Global attribtues now preserved except when they are data dependent CF attributes (which need recalculating). So far this is only implemented for gridseries to gridseries

  • 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 getAxisAttributes(self):
446        ''' get all attribute names as list '''
447        return self.axisobj.attributes.keys()
448       
449    def getTimeUnits(self):
450        ''' get time units - only works if set axis is time...'''
451        #this does the same as getAxisAttribute, but is a separate function as different formats handle time differently.
452        return self.getAxisAttribute('units')
453
454    def getDataForAxis(self):
455        ''' get data for entire axis'''
456        data = self.axisobj.getValue()
457        return data
458
459    def setVariable(self,varname):
460        ''' sets variable '''
461        self.varobj=self.file.variables[varname]
462
463    def getVariableAxes(self):
464        ''' gets axes for set variable'''
465        varAxesList=self.varobj.getAxisIds()
466        return varAxesList
467   
468    def getVariableAttributes(self):
469        '''gets names of all attributes belonging to variable'
470        returns list'''
471        return self.varobj.listattributes()
472
473    def getVariableAttribute(self,attName):
474        ''' gets requested attribute for set variable'''
475        if attName == '_FillValue':
476            try:
477                attribValue=self.varobj._FillValue
478                attribValue=attribValue.item()
479            except:
480                try:
481                    attribValue=self.varobj.missing_value
482                    attribValue=attribValue.item()
483                except:
484                    attribValue = None 
485        else:
486            attribValue = self.varobj.getattribute(attName)
487        return attribValue
488
489    def getDataForVar(self):
490        ''' get data for set variable'''
491        data = self.varobj.getValue()
492        return data
493   
494    def _fixLongitudeRequest(self, **kwargs):
495        ''' fixes problem of longitude requests taking different forms'''
496        lonkey='longitude'
497        if lonkey in kwargs.keys():   #this test needs to be much more robust...!
498            if  type(kwargs[lonkey]) is tuple:
499                newkws=[]
500                for val in kwargs[lonkey]:
501                    newval=((val + 180) % 360) - 180
502                    newkws.append(newval)
503               
504                 #lon = ((t_lon + 180) % 360) - 180
505                #lonmin=kwargs[lonkey][0] - (kwargs[lonkey][0]/180)*360.0
506                #lonmax=kwargs[lonkey][1] - (kwargs[lonkey][1]/180)*360.0
507               
508                #if newkws[0] > newkws[1]:
509                    #flip=[newkws[1], newkws[0]]
510                    #newkws=flip
511                kwargs[lonkey]=tuple(newkws)
512        return kwargs
513               
514    def _componentTimes(self, times):
515        ''' handles time formats'''
516        if type(times) is tuple:
517            comptimes=[]
518            for t in times:
519                compt=csml.csmllibs.csmltime.getCDtime(t)
520                comptimes.append(compt)
521            result=tuple(comptimes)
522        else: #just a single time
523            result=csml.csmllibs.csmltime.getCDtime(times)
524        return result
525       
526   
527    def getSubsetOfDataForVar(self, **kwargs):     
528        ''' Main subsetting method. Kwargs contains dictionary of requested subset.'''
529        #put any slicing indices aside for later and use names       
530        try:
531            upper=kwargs['upper']
532            del kwargs['upper']
533        except:
534            pass
535        try:
536            lower=kwargs['lower']
537            del kwargs['lower']
538        except:
539            pass
540        if 'time' in kwargs:
541            kwargs['time']=self._componentTimes(kwargs['time'])
542        #kwargs=self._fixLongitudeRequest(**kwargs)
543        #return self.file(self.varobj.id,**kwargs)
544        try:
545            #takes keyword args defining subset eg
546            #subset=getSubsetOfDataForVar(latitude=(0.,10.0), longitude=(90, 100.0))
547            subset=None
548            lonkey=None
549            for possiblelonkey in kwargs.keys():
550                if possiblelonkey in ['longitude', 'lon', 'longitude0', 'longitude1']:
551                    lonkey=possiblelonkey
552            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
553                if  type(kwargs[lonkey]) is tuple:
554                    if kwargs[lonkey][0] > kwargs[lonkey][1]:
555                        #subsetting greenwich meridian around 0
556                        lonMin = kwargs[lonkey][0]
557                        lonMax =kwargs[lonkey][1]
558                        kwargs[lonkey]=(0.1, lonMax)
559                        subset1=self.file(self.varobj.id,**kwargs)
560                        kwargs[lonkey]=(lonMin,360)
561                        try:
562                            subset2=self.file(self.varobj.id,**kwargs)
563                            longitudeAxis=subset1.getAxisIndex(lonkey)
564                            #concatenate arrays along longitude             
565                            subset = cdms.MV2.concatenate([subset2,subset1],axis=longitudeAxis)
566                        except:
567                            subset=subset1
568            if type(subset) is not cdms.tvariable.TransientVariable:
569                subset=self.file(self.varobj.id,**kwargs)
570        except:             
571                return None
572                #try and slice with indices instead
573        return subset
574   
575   
576    def getArraySizeOfVar(self):
577        '''iterates through all dimensions in variable to get array size i.e a 3x3x3 array would have a size of 27'''
578        var = self.varobj
579        size = var.shape
580        varsize = 1
581        for item in size:
582            varsize = item *varsize
583        return varsize
584
585    def getShapeOfVar(self):
586        ''' get shape of set variable '''
587        varShape = []
588        for item in self.varobj.shape:
589            varShape.append(item)
590        return varShape
591
592    def getLowLimits(self):
593        ''' get gml low limits of set variable'''
594        dimNames = self.varobj.getAxisIds()
595        lowlims = ""
596        for i in range (1, len(dimNames)):
597            #for now, assume low limit is always of form 1 1 1 ..
598            lowlims =lowlims + str(1)  +' '
599        return lowlims
600
601    def getHighLimits(self):
602        ''' get gml high limits of set variable'''
603        dimNames = self.varobj.getAxisIds()
604        highlims = ""
605        for i in range (1, len(dimNames)):
606            dimValue = self.file.dimensions[dimNames[i]]
607            highlims =highlims  + str(dimValue) +' '
608        return highlims
609       
610    def getGlobalAttributes(self):
611        '''return dictionary of 'global' attributes'''
612        atts= self.file.attributes.keys()
613        globalatts={}
614        for i in range(0,len(atts)):
615            globalatts[atts[i]]=self.file.attributes[atts[i]]
616        return globalatts
617       
618   
619class cdmlInterface(cdunifInterface):
620    ''' this is more  or less the cdunif interface but a few methods have been overwritten to read CDML'''
621    def __init__(self):
622        #this all needs to be revisited in csml v2.
623        self.extractType='cdmlExtract'
624        self.extractPrefix = '_cdmlextract_'
625       
626    def getListOfAxes(self):
627        ''' get list of axis in file'''
628        axes=self.file.axes.keys() 
629        return axes
630
631    def getSizeOfAxis(self,axis):
632        ''' get size of set axis'''
633        axisSize=self.file.axes[axis].length
634        return axisSize
635
636
637class RawFileInterface(AbstractDI):
638    ''' interface for raw file types'''
639    def __init__(self):
640        self.extractType   = 'rawExtract'
641        self.extractPrefix = '_rawextract_'
642 
643 
644    def openFile(self, filename):
645        self.file = open(filename, "rb")
646
647    def closeFile(self):
648            self.file.close()
649
650
651    # Read the data from the raw file into a multidimensional Numeric array
652    def readFile(self, **kwargs):
653        # Determine the numeric type:
654        if 'numericType' in kwargs:
655            try:
656                numericTypeCode = {
657                    'uint8':Numeric.UInt8,
658                    'uint16':Numeric.UInt16,
659                    'uint32':Numeric.UInt32,
660                    'int8':Numeric.Int8,
661                    'int16':Numeric.Int16,
662                    'int32':Numeric.Int32,
663                    'float':Numeric.Float32,
664                    'double':Numeric.Float64
665                }[kwargs['numericType']]
666            except KeyError:
667                raise TypeError("Invalid numericType: " + str(kwargs['numericType']))
668        else:
669            raise KeyError("numericType is mandatory.")
670       
671        # Read the file into a numpy array:
672        self.data = Numeric.fromstring(self.file.read(), numericTypeCode)
673        # If necessary, swap the byte order:
674        if 'endianess' in kwargs:
675            if ((kwargs['endianness'] == 'little' and not Numeric.LittleEndian) or (kwargs['endianness'] == 'big' and Numeric.LittleEndian)):
676                self.data = self.data.byteswap()   
677        # Declare the shape of the array:
678        dimensions = map(int,kwargs['dimensions'])
679        dimensions.reverse()
680        self.data.shape = tuple(dimensions)
681        # If numericTransform or fillValue were provided, store them as
682        # attributes.
683        if 'numericTransform' in kwargs:
684            self.numericTransform = NumericTransform.infixExpression(kwargs['numericTransform'])
685        if 'fillValue' in kwargs:
686            self.fillValue = kwargs['fillValue']
687   
688   # Return the fill value, if set, and transform if necessary:
689    def getFillValue(self):
690      # Both fillValue and numericTransform attributes may or may
691      # not exist...
692        try:
693            return self.numericTransform.solve( n = float(self.fillValue) )
694        except AttributeError:
695            try:
696                return self.fillValue
697            except AttributeError:
698                return None
699
700
701   # This is a just a special case of getSubsetOfDataForVar. To avoid
702   # duplication of code, just subset the entire array. (getSubset.. is
703   # optimised for this case)
704    def getDataForVar(self):
705        return self.getSubsetOfDataForVar(lower = (0,)*len(self.data.shape),
706                                        upper = self.data.shape)
707
708
709   # Subset the n-dimensional file based on array indices. Accepts parameters
710   # in the form of e.g.
711   #
712   # getSubsetOfDataForVar( lower=(0,0), upper=(1,2) )
713   #
714   # Note: The rank of the required subset, must exactly match the
715   # rank of the original data: len(lower) == len(upper) == rank of file
716   #
717    def getSubsetOfDataForVar(self, **kwargs):
718        # Assume subset parameters are passed as: lower=(0,0) upper=(512,512)
719        if 'lower' not in kwargs or 'upper' not in kwargs:
720            # Have not specified any subset parameters that we recognise, so raise
721            # an exception:
722            raise NotImplementedError("Only supports subsetting with lower+upper array indices")
723        elif not len(kwargs['lower']) == len(kwargs['upper']) == len(self.data.shape):
724            # Rank of subset is not the same as rank of full data array. so raise
725            # an exception:
726            raise NotImplementedError("Only supports subsets of same rank as full dataset")
727        elif Numeric.sometrue(Numeric.greater(kwargs['upper'], self.data.shape)):
728            # Requested upper bound of subset is beyond the size of the the full
729            # data array, so raise an exception
730            raise IndexError("Subset out of range")
731        elif Numeric.sometrue(Numeric.less( kwargs['upper'], Numeric.zeros(len(self.data.shape)))):
732            # Requested lower bound of subset is beyond the size of the the full
733            # data array, so raise an exception
734            raise IndexError("Subset out of range")
735        elif Numeric.sometrue(Numeric.less_equal(kwargs['upper'], kwargs['lower'])):
736            # lower bound <= upper_bound for at least one dimension, so raise an
737            # exception
738            raise IndexError("Upper bound less than lower bound")
739        elif tuple(kwargs['lower']) == (0,)*len(self.data.shape) and tuple(kwargs['upper']) == self.data.shape:
740            # Special case of requested subset == entire data file.
741            subset = self.data
742        else:
743            # We are okay to subset.
744   
745            # I cant see any nice (and speedy) way of subsetting a Numeric
746            # array of arbitrary rank without the use of eval. By doing it
747            # this way, we can make use of the (possible) acceleration provided
748            # by Numeric/NumPy.
749            slices = []
750            for i in range(len(self.data.shape)):
751                lower = int(kwargs['lower'][i])
752                upper = int(kwargs['upper'][i]) 
753                slices.append(str(lower)+':'+str(upper))
754            subset = eval("self.data["+','.join(slices)+"]")
755   
756        # Attempt to perform the numericTransform on the data array, if we get
757        # AttributeError, it most likely means that numericTransform was not
758        # specified, so return the data as-is
759        try:
760            return self.numericTransform.solve( n = subset )
761        except AttributeError:
762            return subset.copy()
763
764class ImageFileInterface(RawFileInterface):
765    '''Interface for reading data from image files. Requires PIL Image module.'''
766    def __init__(self):
767        self.extractType   = 'imageExtract'
768        self.extractPrefix = '_imageextract_'
769   
770    def image2array(self,im):
771        #Adapted from code by Fredrik Lundh, http://www.pythonware.com
772        #http://effbot.org/zone/pil-numpy.htm
773            if im.mode not in ("L", "F"):
774                raise ValueError, "can only convert single-layer images"
775            if im.mode == "L":
776                a = Numeric.fromstring(im.tostring(), Numeric.UnsignedInt8)
777            else:
778                a = Numeric.fromstring(im.tostring(), Numeric.Float32)
779            a.shape = im.size[1], im.size[0]
780            return a
781
782    def openFile(self, filename):
783        self.file = Image.open(filename)
784
785    def closeFile(self):
786        self.file = None #...Image does not seem to have a close() method.
787
788    def readFile(self, **kwargs):
789        # Convert the image to a Numeric array
790     
791        self.data=self.image2array(self.file)
792        #slower method:
793        #self.data = Numeric.array(self.file.getdata())
794
795        if 'numericTransform' in kwargs:
796            # numericTransform was specified, so compile the expression:
797            self.numericTransform = NumericTransform.infixExpression(kwargs['numericTransform'])
798        if 'fillValue' in kwargs:
799            self.fillValue = kwargs['fillValue']
800
801class HDF5Interface(AbstractDI):
802    '''Data Interface for HDF5'''
803    def __init__(self):
804        self.extractType='hdfExtract'
805        self.extractPrefix='_hdfextract_'
806   
807    def openFile(self, filename):
808        '''some code to open the file'''
809        pass
810
811    def setAxis(self,axis):
812        '''some code to set an axis to be queried, may not need to do much, depending on your format'''
813        pass
814       
815    def getDataForAxis(self):
816        '''some code to return the values for an axis'''
817        return data
818
819    def setVariable(self,varname):
820        '''some code to set a variable to be queried, may not need to do much, depending on your format'''
821        pass
822
823    def getDataForVar(self):
824        '''some code to return all values for a variable'''
825        return data
826
827    def getSubsetOfDataForVar(self, **kwargs):
828        '''takes keyword args defining subset eg
829        subset=getSubsetOfDataForVar(latitude=(0.,10.0), longitude=(90, 100.0), ...)
830        and returns a subset of data for tha variable '''
831        return data
832
833    def closeFile(self):
834        '''some code to close the file'''
835        pass
Note: See TracBrowser for help on using the repository browser.