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

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

fixed bug where longitude requests were not always working as expected

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