Changeset 4270 for cows/trunk/cows


Ignore:
Timestamp:
02/10/08 12:08:51 (11 years ago)
Author:
domlowe
Message:

moving wms_csmllayer code

Location:
cows/trunk/cows/service/imps
Files:
2 edited
1 moved

Legend:

Unmodified
Added
Removed
  • cows/trunk/cows/service/imps/csmlbackend/wms_csmllayer.py

    r4269 r4270  
     1""" 
     2implementation of ILayerMapper, ILayer, IDimension, ILayerSlab interfaces, as defined in wms_iface.py 
     3 
     4""" 
     5import os, string 
     6import csml 
     7try: 
     8    import cdms2 as cdms 
     9except: 
     10    import cdms 
     11import Image 
     12from copy import copy 
     13from pywms.render_imp import RGBARenderer 
     14from matplotlib import cm 
     15import genutil 
     16from pylons import config, request, session  #config must have tmpfilebuffer and csmlstore values 
     17import cdtime 
     18import logging 
     19log = logging.getLogger(__name__) 
     20import ConfigParser 
     21 
     22from cows.service.wxs_iface import ILayerMapper 
     23from cows.service.wms_iface import IwmsLayer, IwmsDimension, IwmsLayerSlab 
     24 
     25try: 
     26    from ndgUtils import ndgObject, ndgRetrieve 
     27except: 
     28    log.warning("ndgUtils library could not be loaded, files in the eXist database won't be available, although you should still be able to access files from the csmlstore directory referenced in the ini file.") 
     29 
     30 
     31 
     32 
     33class CSMLLayerMapper(ILayerMapper): 
     34    """ 
     35    Map keyword arguments to a collection of layers. 
     36    Supports the retrieval of sets of layers according to arbitary 
     37    keyword/value pairs. 
     38    Implements  ILayerMapper  
     39     
     40    """ 
     41    def __init__(self): 
     42        self.layermapcache={} 
     43     
     44    def _getInfo(self, feature): 
     45        ''' given a csml feature, return info about the layer/feature 
     46        @return:   title, abstract, dimensions, units, crss ''' 
     47 
     48        try: 
     49            title=feature.name.CONTENT 
     50        except: 
     51            title='' 
     52        try: 
     53            abstract=feature.description.CONTENT 
     54        except: 
     55            abstract=title 
     56         
     57        units=feature.getDomainUnits() 
     58        dimensions={} 
     59        tmpunits=copy(units) 
     60        tmpunits.reverse() 
     61        domain = feature.getDomain() 
     62        for dim in feature.getAxisLabels(): 
     63            nextdim=CSMLDimension(domain, dim, tmpunits.pop()) 
     64            if dim not in ['latitude', 'longitude']: 
     65                dimensions[dim]=nextdim 
     66        crs=feature.getNativeCRS() 
     67        crss=[self._crscat.getCRS(crs).twoD] 
     68        if 'EPSG:4326' in crss: 
     69            crss.append('CRS:84') 
     70            crss.append('WGS84') 
     71        return title, abstract, dimensions, units, crss 
     72         
     73         
     74 
     75    def _get_csml_doc(self, fileoruri): 
     76        """ 
     77        Gets a csml document from file or exist when passed a file name or uri          
     78                
     79        Note, access control is not implemented on file access, only on exist access. 
     80         
     81        """ 
     82        if string.find(fileoruri,'__NDG-A0__') == -1: 
     83            #it's a local file not an identifier, read from local csmlstore 
     84            file=fileoruri 
     85            csml_dir = config['csmlstore'] 
     86            path = os.path.join(csml_dir, file) 
     87            if os.path.exists(path+'.csml'): 
     88                f = path+'.csml' 
     89            elif os.path.exists(path+'.xml'): 
     90                f = path +'.xml' 
     91            else: 
     92                raise ValueError("Cannot find CSML file %s" % file) 
     93             
     94            d = csml.parser.Dataset() 
     95            d.parse(f) 
     96         
     97        else: 
     98            #it's an NDG identifier, get the document from exist. 
     99            uri=fileoruri  
     100            uriN=ndgObject(uri) 
     101            conflocation=config['ndgconfig'] 
     102            cf=ConfigParser.ConfigParser() 
     103            cf.read(conflocation) 
     104            status,x=ndgRetrieve(uriN, cf, discovery=0) 
     105            d=csml.parser.Dataset() 
     106            d.parse(x)                          
     107        return d 
     108     
     109    def map(self, **kwargs): 
     110        """ 
     111        Given csml.parser.Dataset object list the names of 
     112        all layers available. 
     113         
     114        @return: A mapping of layer names to ILayer implementations. 
     115        @raise ValueError: If no layers are available for these keywords.  
     116        """ 
     117        fileoruri=kwargs['fileoruri'] 
     118        if fileoruri in self.layermapcache.keys(): 
     119            #we've accessed this layer map before, get it from the cache dictionary 
     120            return self.layermapcache[fileoruri] 
     121          
     122        ds = self._get_csml_doc(fileoruri) 
     123        layermap={} 
     124        self._crscat=csml.csmllibs.csmlcrs.CRSCatalogue() 
     125        for feature in csml.csmllibs.csmlextra.listify(ds.featureCollection.featureMembers): 
     126            title, abstract, dimensions, units, crss=self._getInfo(feature) 
     127            layermap[feature.id]=CSMLLayer(title,abstract, dimensions, units, crss, feature) 
     128        if len(layermap) > 0: 
     129            self.layermapcache[fileoruri]=layermap 
     130            return layermap 
     131        else: 
     132            raise ValueError 
     133 
     134 
     135class CSMLLayer(IwmsLayer): 
     136    """ 
     137     representing a WMS layer.    Implements ILayer 
     138 
     139    @ivar title: The layer title.  As seen in the Capabilities document. 
     140    @ivar abstract:  Abstract as seen in the Capabilities document. 
     141    @ivar dimensions: A dictionary of IDimension objects. 
     142    @ivar units: A string describing the units. 
     143    @ivar crss: A sequence of SRS/CRSs supported by this layer. 
     144 
     145    @todo: Do we need minValue/maxValue? 
     146 
     147    """ 
     148     
     149    def __init__(self, title, abstract, dimensions, units, crss, feature): 
     150        self.featureInfoFormats=None #NotImplemented 
     151        self.title=title 
     152        self.abstract=abstract 
     153        self.dimensions=dimensions 
     154        self.units=units 
     155        self.crss=crss 
     156        self._feature=feature 
     157        self.legendSize=(30,100) 
     158        bb= self._feature.getCSMLBoundingBox().getBox() 
     159        #convert 0 - 360 to -180, 180 as per common WMS convention 
     160        if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361: 
     161            bb[0], bb[2]=-180, 180 
     162        self.wgs84BBox = bb 
     163        self.featureInfoFormats = ['text/html'] 
     164        try: 
     165            self.wgs84BBox = self.getBBox('EPSG:4326') 
     166        except: 
     167            raise ValueError("Layer must provide a bounding box in EPSG:4326 " 
     168                             "coordinates for compatibility with WMS-1.3.0") 
     169        self.featureinfofilecache={} #used for caching netcdf file in getFeatureInfo 
     170         
     171    def getBBox(self, crs): 
     172        """ 
     173        @return: A 4-typle of the bounding box in the given coordinate 
     174            reference system. 
     175        """ 
     176        #bb= self._feature.getCSMLBoundingBox().getBox() 
     177        #convert 0 - 360 to -180, 180 as per common WMS convention 
     178        #if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361: 
     179        #    bb[0], bb[2]=-180, 180 
     180        #self.wgs84BBox = bb 
     181        return self.wgs84BBox 
     182        #raise NotImplementedError 
     183         
     184    def getSlab(self, crs, dimValues=None, renderOpts={}): 
     185        """ 
     186        Creates a slab of the layer in a particular CRS and set of 
     187        dimensions. 
     188 
     189        @param crs: The coordinate reference system. 
     190        @param dimValues: A mapping of dimension names to dimension values 
     191            as specified in the IDimension.extent 
     192        @param renderOpts: A generic mapping object for passing rendering 
     193            options 
     194        @return: An object implementing ILayerSlab 
     195        #create netcdf for whole lat/lon for given dimValues, use to init slab 
     196        """ 
     197 
     198        log.debug('getSlab(%s, %s)' % (crs, dimValues)) 
     199         
     200        #unicode conversion 
     201        for dimval in dimValues: 
     202            if dimval != 'time': 
     203                dimValues[dimval]=float(dimValues[dimval]) 
     204            else: 
     205                #remove any trailing Zs from time string 
     206                if dimValues[dimval] [-1:] in ['Z', 'z']: 
     207                    dimValues[dimval]=dimValues[dimval][:-1] 
     208        if type(self._feature) == csml.parser.GridSeriesFeature: 
     209            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc' 
     210            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues) 
     211            #for now have to read netcdf back from  
     212            #disk (limitiation of CSML api) 
     213            netcdf=cdms.open(result[1]) 
     214            #and then delete the temporary file 
     215            os.system('rm %s'%result[1]) 
     216            bbox=self.getBBox(crs) 
     217            return CSMLLayerSlab(netcdf, self, crs, dimValues, renderOpts, bbox) 
     218        else: 
     219            raise NotImplementedError 
     220         
     221    def getCacheKey(self, crs, dimValues=None, renderOpts={}): 
     222        """ 
     223        Create a unique key for use in caching a slab. 
     224 
     225        The intention here is that most of the work should be done when 
     226        instantiating an ILayerSlab object.  These can be cached by the 
     227        server for future use.  The server will first call getCacheKey() 
     228        for the slab creation arguments and if the key is in it's cache 
     229        it will use a pre-generated ILayerSlab object. 
     230 
     231        """ 
     232        dimList = list(dimValues.items()) 
     233        dimList.sort() 
     234        return '%s:%s' % (crs, dimList) 
     235 
     236    def getFeatureInfo(self, format, crs, point, dimValues): 
     237        """ 
     238        Return a response string descibing the feature at a given 
     239        point in a given CRS. 
     240 
     241        Currently only "html" is supported as output format 
     242 
     243        @param format: One of self.featureInfoFormats.  Defines which 
     244            format the response will be in. 
     245        @param crs: One of self.crss 
     246        @param point: a tuple (x, y) in the supplied crs of the point 
     247            being selected. 
     248        @param dimValues: A mapping of dimension names to dimension values. 
     249        @return: A string containing the response. 
     250 
     251        """ 
     252         
     253        #cached netcdf is indexed by a tuple of feature id and dimvalues - i.e. typically a particular time and Z value for that feature. 
     254        #look in dictionary for cached copy, and if so use that as the file object. 
     255        dictindex=str((self._feature.id, dimValues)) 
     256        if dictindex in self.featureinfofilecache: 
     257            log.debug('calling cache') 
     258            f=self.featureinfofilecache[dictindex] 
     259        else: #else, use the csml api to subset the feature afresh 
     260            log.debug('not calling cache') 
     261            randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc' 
     262            result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues) 
     263            #for now have to read netcdf back from disk (limitation of CSML api) 
     264            f=cdms.open(result[1]) 
     265            #append to cache: 
     266            self.featureinfofilecache[dictindex]=f 
     267            #and then delete the temporary file 
     268            os.system('rm %s'%result[1]) 
     269         
     270        netcdf = f(self.title)  #netcdf here is a cdms transient variable 
     271         
     272        
     273        #Now grab the netCDF object for the point specified. 
     274        #The reason for the 'cob' option is so that if the grid the data  
     275        #is defined on does not have a grid point at the point specified,  
     276        #we should  still get the nearest location 
     277         
     278        t_point = netcdf(latitude=(point[1], point[1], 'cob'), longitude=(point[0], point[0], 'cob')) 
     279        #now get the value recorded at this location 
     280        value = t_point.getValue().tolist() 
     281        log.debug(value) 
     282        log.debug(t_point.fill_value()) 
     283        #and the fill_value too 
     284        fill_value = t_point.fill_value() 
     285        #value is actually embedded in a multi dimensional list,  
     286        #so we need to extract the actual value from the list 
     287        while type(value) is list: 
     288                value = value[0] 
     289 
     290        #now check if the value is actually the fill_value rather than  
     291        #a value recorded at the point specified 
     292        log.debug('%s %s' % (value, fill_value)) 
     293        if (2*fill_value) == value: 
     294                value = "No value found at position: "+str(point[1])+", "+str(point[0]) 
     295        else: 
     296                value = "Value found at position: "+str(point[1])+", "+str(point[0])+" is: "+str(value) 
     297        # finally return the value 
     298        return value 
     299 
     300 
     301class CSMLDimension(IwmsDimension): 
     302    """ 
     303    implements IDimension 
     304    @ivar units: The units string. 
     305    @ivar extent: Sequence of extent values. 
     306 
     307    """ 
     308     
     309    def __init__(self, domain, dimname, unit): 
     310        self.units = unit 
     311        self.extent = [] 
     312        #patch to handle current limitations of multiple time dimension scanning in csml.  
     313        if string.lower(self.units)[:10] in ['days_since', 'seconds_si', 'minutes_si', 'hours_sinc','months _sin', 'years_sinc']: 
     314            if type(domain[dimname][0]) is not str   : 
     315                tunits=self.units.replace('_', ' ') 
     316                for val in domain[dimname]: 
     317                    csmltime= csml.csmllibs.csmltime.UDtimeToCSMLtime(cdtime.reltime(float(val), tunits).tocomp()) 
     318                    self.extent.append(csmltime) 
     319            else: 
     320                for val in domain[dimname]: 
     321                    self.extent.append(str(val)) 
     322        else: 
     323            for val in domain[dimname]: 
     324                self.extent.append(str(val)) 
     325        #for time axis replace units with iso string 
     326        if dimname == 'time': 
     327            self.units='ISO8601' 
     328        
     329class CSMLLayerSlab(IwmsLayerSlab): 
     330    """ 
     331    Implements LayerSlab 
     332    Represents a particular horizontal slice of a WMS layer. 
     333 
     334    ILayerSlab objects are designed to be convenient to cache. 
     335    They should be pickleable to enable memcached support in the future. 
     336 
     337    @ivar layer: The source ILayer instance. 
     338    @ivar crs: The coordinate reference system. 
     339    @ivar dimValues: A mapping of dimension values of this view. 
     340    @ivar renderOpts: The renderOpts used to create this view. 
     341    @ivar bbox: The bounding box as a 4-tuple. 
     342    """ 
     343     
     344    def __init__(self, netcdf, layer, crs, dimValues, renderOpts, bbox): 
     345        self._netcdf=netcdf 
     346        self.layer = layer 
     347        self.crs = crs 
     348        self.dimValues = dimValues 
     349        self.renderOpts=renderOpts 
     350        self.bbox=bbox 
     351         
     352        #set colour map for ALL images from this slab 
     353        v=self._netcdf(layer.title) 
     354        tvar=v(squeeze=1) 
     355        #get the min and max values to use for the colourmapping. 
     356        #change the fill values to ensure they aren't picked up as false max/mins. 
     357        tvar.missing_value=999999999 
     358        value=tvar.getValue() 
     359        self.minval=min(min(l) for l in value) 
     360        tvar.missing_value=-0.000000001 
     361        value=tvar.getValue() 
     362        self.maxval=max(max(l) for l in value) 
     363 
     364         
     365    def getImage(self, bbox, width, height): 
     366        """ 
     367        Create an image of a sub-bbox of a given size. 
     368 
     369        @ivar bbox: A bbox 4-tuple. 
     370        @ivar width: width in pixels.`   
     371        @ivar height: height in pixels. 
     372        @return: A PIL Image object. 
     373 
     374        """ 
     375 
     376        log.debug('getImage(%s, %s, %s)' % (bbox, width, height)) 
     377         
     378 
     379        cmap=eval(config['colourmap']) # renderOpts is hook for colourmap, for now use config 
     380        grid=Grid(self.layer, self._netcdf, bbox, width, height) 
     381        #how to handle varmin,varmax? ..read array? 
     382        #minval, maxval=genutil.minmax(grid.value) 
     383        #minval=min(min(l) for l in grid.value) 
     384        #maxval=max(max(l) for l in grid.value) 
     385        minval=self.minval 
     386        maxval=self.maxval 
     387        renderer=RGBARenderer(minval, maxval)          
     388        return renderer.renderGrid(grid, bbox, width, height, cmap) 
     389     
     390class Grid(object): 
     391    """A class encapsulating a simple regularly spaced, rectilinear 
     392    grid.  This is the only type of grid pywms is expected to 
     393    understand and adaptors should be provided to connect to 
     394    underlying implementations such as cdms or csml. 
     395 
     396    @cvar crs: Coordinate reference system 
     397 
     398    @ivar x0: coordinate of the lower left corner. 
     399    @ivar y0: coordinate of the lower left corner. 
     400    @ivar dx: The x grid spacing. 
     401    @ivar dy: The y grid spacing. 
     402    @ivar nx: The number of x grid points. 
     403    @ivar ny: The number of y grid points. 
     404    @ivar value: A masked array of the grid values. 
     405    @ivar ix: The dimension index of longidude in value 
     406    @ivar iy: The dimension index of latitude in value 
     407    @ivar long_name: The name of the field. 
     408    @ivar units: The units of the field. 
     409    """ 
     410    def __init__(self, layer, netcdf, bbox, width, height): 
     411        #we know the axes are called latitude and longitude as the CSML code has written it: 
     412        v=netcdf(layer.title) 
     413        tvar=v(latitude=(bbox[1], bbox[3]), longitude=(bbox[0],bbox[2]),squeeze=1) 
     414        order=tvar.getOrder() 
     415        #array of data 
     416        self.value=tvar.getValue() 
     417        if order == 'xy': 
     418            self.ix=0 
     419            self.iy=1 
     420        else: 
     421            self.ix=1 
     422            self.iy=0 
     423        lat = tvar.getLatitude() 
     424        lon = tvar.getLongitude() 
     425        self.x0=lon[0] 
     426        self.y0=lat[0] 
     427        self.dx=abs(lon[0]-lon[1]) 
     428        self.dy=abs(lat[0]-lat[1]) 
     429         
     430     
     431 
     432         
     433 
     434        self.nx=len(lon) 
     435        self.ny=len(lat) 
     436        self.long_name=tvar.id  #TODO, get long name from feature 
     437        self.units=tvar.units 
  • cows/trunk/cows/service/imps/wms_csmllayer.py

    r4266 r4270  
    1 """ 
    2 implementation of ILayerMapper, ILayer, IDimension, ILayerSlab interfaces, as defined in wms_iface.py 
    3  
    4 """ 
    5 import os, string 
    6 import csml 
    7 try: 
    8     import cdms2 as cdms 
    9 except: 
    10     import cdms 
    11 import Image 
    12 from copy import copy 
    13 from pywms.render_imp import RGBARenderer 
    14 from matplotlib import cm 
    15 import genutil 
    16 from pylons import config, request, session  #config must have tmpfilebuffer and csmlstore values 
    17 import cdtime 
    18 import logging 
    19 log = logging.getLogger(__name__) 
    20 import ConfigParser 
    21  
    22 from cows.service.wxs_iface import ILayerMapper 
    23 from cows.service.wms_iface import IwmsLayer, IwmsDimension, IwmsLayerSlab 
    24  
    25 try: 
    26     from ndgUtils import ndgObject, ndgRetrieve 
    27 except: 
    28     log.warning("ndgUtils library could not be loaded, files in the eXist database won't be available, although you should still be able to access files from the csmlstore directory referenced in the ini file.") 
    29  
    30  
    31  
    32  
    33 class CSMLLayerMapper(ILayerMapper): 
    34     """ 
    35     Map keyword arguments to a collection of layers. 
    36     Supports the retrieval of sets of layers according to arbitary 
    37     keyword/value pairs. 
    38     Implements  ILayerMapper  
    39      
    40     """ 
    41     def __init__(self): 
    42         self.layermapcache={} 
    43      
    44     def _getInfo(self, feature): 
    45         ''' given a csml feature, return info about the layer/feature 
    46         @return:   title, abstract, dimensions, units, crss ''' 
    47  
    48         try: 
    49             title=feature.name.CONTENT 
    50         except: 
    51             title='' 
    52         try: 
    53             abstract=feature.description.CONTENT 
    54         except: 
    55             abstract=title 
    56          
    57         units=feature.getDomainUnits() 
    58         dimensions={} 
    59         tmpunits=copy(units) 
    60         tmpunits.reverse() 
    61         domain = feature.getDomain() 
    62         for dim in feature.getAxisLabels(): 
    63             nextdim=CSMLDimension(domain, dim, tmpunits.pop()) 
    64             if dim not in ['latitude', 'longitude']: 
    65                 dimensions[dim]=nextdim 
    66         crs=feature.getNativeCRS() 
    67         crss=[self._crscat.getCRS(crs).twoD] 
    68         if 'EPSG:4326' in crss: 
    69             crss.append('CRS:84') 
    70             crss.append('WGS84') 
    71         return title, abstract, dimensions, units, crss 
    72          
    73          
    74  
    75     def _get_csml_doc(self, fileoruri): 
    76         """ 
    77         Gets a csml document from file or exist when passed a file name or uri          
    78                 
    79         Note, access control is not implemented on file access, only on exist access. 
    80          
    81         """ 
    82         if string.find(fileoruri,'__NDG-A0__') == -1: 
    83             #it's a local file not an identifier, read from local csmlstore 
    84             file=fileoruri 
    85             csml_dir = config['csmlstore'] 
    86             path = os.path.join(csml_dir, file) 
    87             if os.path.exists(path+'.csml'): 
    88                 f = path+'.csml' 
    89             elif os.path.exists(path+'.xml'): 
    90                 f = path +'.xml' 
    91             else: 
    92                 raise ValueError("Cannot find CSML file %s" % file) 
    93              
    94             d = csml.parser.Dataset() 
    95             d.parse(f) 
    96          
    97         else: 
    98             #it's an NDG identifier, get the document from exist. 
    99             uri=fileoruri  
    100             uriN=ndgObject(uri) 
    101             conflocation=config['ndgconfig'] 
    102             cf=ConfigParser.ConfigParser() 
    103             cf.read(conflocation) 
    104             status,x=ndgRetrieve(uriN, cf, discovery=0) 
    105             d=csml.parser.Dataset() 
    106             d.parse(x)                          
    107         return d 
    108      
    109     def map(self, **kwargs): 
    110         """ 
    111         Given csml.parser.Dataset object list the names of 
    112         all layers available. 
    113          
    114         @return: A mapping of layer names to ILayer implementations. 
    115         @raise ValueError: If no layers are available for these keywords.  
    116         """ 
    117         fileoruri=kwargs['fileoruri'] 
    118         if fileoruri in self.layermapcache.keys(): 
    119             #we've accessed this layer map before, get it from the cache dictionary 
    120             return self.layermapcache[fileoruri] 
    121           
    122         ds = self._get_csml_doc(fileoruri) 
    123         layermap={} 
    124         self._crscat=csml.csmllibs.csmlcrs.CRSCatalogue() 
    125         for feature in csml.csmllibs.csmlextra.listify(ds.featureCollection.featureMembers): 
    126             title, abstract, dimensions, units, crss=self._getInfo(feature) 
    127             layermap[feature.id]=CSMLLayer(title,abstract, dimensions, units, crss, feature) 
    128         if len(layermap) > 0: 
    129             self.layermapcache[fileoruri]=layermap 
    130             return layermap 
    131         else: 
    132             raise ValueError 
    133  
    134  
    135 class CSMLLayer(IwmsLayer): 
    136     """ 
    137      representing a WMS layer.    Implements ILayer 
    138  
    139     @ivar title: The layer title.  As seen in the Capabilities document. 
    140     @ivar abstract:  Abstract as seen in the Capabilities document. 
    141     @ivar dimensions: A dictionary of IDimension objects. 
    142     @ivar units: A string describing the units. 
    143     @ivar crss: A sequence of SRS/CRSs supported by this layer. 
    144  
    145     @todo: Do we need minValue/maxValue? 
    146  
    147     """ 
    148      
    149     def __init__(self, title, abstract, dimensions, units, crss, feature): 
    150         self.featureInfoFormats=None #NotImplemented 
    151         self.title=title 
    152         self.abstract=abstract 
    153         self.dimensions=dimensions 
    154         self.units=units 
    155         self.crss=crss 
    156         self._feature=feature 
    157         self.legendSize=(30,100) 
    158         bb= self._feature.getCSMLBoundingBox().getBox() 
    159         #convert 0 - 360 to -180, 180 as per common WMS convention 
    160         if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361: 
    161             bb[0], bb[2]=-180, 180 
    162         self.wgs84BBox = bb 
    163         self.featureInfoFormats = ['text/html'] 
    164         try: 
    165             self.wgs84BBox = self.getBBox('EPSG:4326') 
    166         except: 
    167             raise ValueError("Layer must provide a bounding box in EPSG:4326 " 
    168                              "coordinates for compatibility with WMS-1.3.0") 
    169         self.featureinfofilecache={} #used for caching netcdf file in getFeatureInfo 
    170          
    171     def getBBox(self, crs): 
    172         """ 
    173         @return: A 4-typle of the bounding box in the given coordinate 
    174             reference system. 
    175         """ 
    176         #bb= self._feature.getCSMLBoundingBox().getBox() 
    177         #convert 0 - 360 to -180, 180 as per common WMS convention 
    178         #if abs(bb[2]-bb[0]) >= 359 and abs(bb[2]-bb[0]) < 361: 
    179         #    bb[0], bb[2]=-180, 180 
    180         #self.wgs84BBox = bb 
    181         return self.wgs84BBox 
    182         #raise NotImplementedError 
    183          
    184     def getSlab(self, crs, dimValues=None, renderOpts={}): 
    185         """ 
    186         Creates a slab of the layer in a particular CRS and set of 
    187         dimensions. 
    188  
    189         @param crs: The coordinate reference system. 
    190         @param dimValues: A mapping of dimension names to dimension values 
    191             as specified in the IDimension.extent 
    192         @param renderOpts: A generic mapping object for passing rendering 
    193             options 
    194         @return: An object implementing ILayerSlab 
    195         #create netcdf for whole lat/lon for given dimValues, use to init slab 
    196         """ 
    197  
    198         log.debug('getSlab(%s, %s)' % (crs, dimValues)) 
    199          
    200         #unicode conversion 
    201         for dimval in dimValues: 
    202             if dimval != 'time': 
    203                 dimValues[dimval]=float(dimValues[dimval]) 
    204             else: 
    205                 #remove any trailing Zs from time string 
    206                 if dimValues[dimval] [-1:] in ['Z', 'z']: 
    207                     dimValues[dimval]=dimValues[dimval][:-1] 
    208         if type(self._feature) == csml.parser.GridSeriesFeature: 
    209             randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc' 
    210             result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues) 
    211             #for now have to read netcdf back from  
    212             #disk (limitiation of CSML api) 
    213             netcdf=cdms.open(result[1]) 
    214             #and then delete the temporary file 
    215             os.system('rm %s'%result[1]) 
    216             bbox=self.getBBox(crs) 
    217             return CSMLLayerSlab(netcdf, self, crs, dimValues, renderOpts, bbox) 
    218         else: 
    219             raise NotImplementedError 
    220          
    221     def getCacheKey(self, crs, dimValues=None, renderOpts={}): 
    222         """ 
    223         Create a unique key for use in caching a slab. 
    224  
    225         The intention here is that most of the work should be done when 
    226         instantiating an ILayerSlab object.  These can be cached by the 
    227         server for future use.  The server will first call getCacheKey() 
    228         for the slab creation arguments and if the key is in it's cache 
    229         it will use a pre-generated ILayerSlab object. 
    230  
    231         """ 
    232         dimList = list(dimValues.items()) 
    233         dimList.sort() 
    234         return '%s:%s' % (crs, dimList) 
    235  
    236     def getFeatureInfo(self, format, crs, point, dimValues): 
    237         """ 
    238         Return a response string descibing the feature at a given 
    239         point in a given CRS. 
    240  
    241         Currently only "html" is supported as output format 
    242  
    243         @param format: One of self.featureInfoFormats.  Defines which 
    244             format the response will be in. 
    245         @param crs: One of self.crss 
    246         @param point: a tuple (x, y) in the supplied crs of the point 
    247             being selected. 
    248         @param dimValues: A mapping of dimension names to dimension values. 
    249         @return: A string containing the response. 
    250  
    251         """ 
    252          
    253         #cached netcdf is indexed by a tuple of feature id and dimvalues - i.e. typically a particular time and Z value for that feature. 
    254         #look in dictionary for cached copy, and if so use that as the file object. 
    255         dictindex=str((self._feature.id, dimValues)) 
    256         if dictindex in self.featureinfofilecache: 
    257             log.debug('calling cache') 
    258             f=self.featureinfofilecache[dictindex] 
    259         else: #else, use the csml api to subset the feature afresh 
    260             log.debug('not calling cache') 
    261             randomname= csml.csmllibs.csmlextra.getRandomID() + '.nc' 
    262             result= self._feature.subsetToGridSeries(config['tmpfilebuffer'], ncname=randomname, **dimValues) 
    263             #for now have to read netcdf back from disk (limitation of CSML api) 
    264             f=cdms.open(result[1]) 
    265             #append to cache: 
    266             self.featureinfofilecache[dictindex]=f 
    267             #and then delete the temporary file 
    268             os.system('rm %s'%result[1]) 
    269          
    270         netcdf = f(self.title)  #netcdf here is a cdms transient variable 
    271          
    272         
    273         #Now grab the netCDF object for the point specified. 
    274         #The reason for the 'cob' option is so that if the grid the data  
    275         #is defined on does not have a grid point at the point specified,  
    276         #we should  still get the nearest location 
    277          
    278         t_point = netcdf(latitude=(point[1], point[1], 'cob'), longitude=(point[0], point[0], 'cob')) 
    279         #now get the value recorded at this location 
    280         value = t_point.getValue().tolist() 
    281         log.debug(value) 
    282         log.debug(t_point.fill_value()) 
    283         #and the fill_value too 
    284         fill_value = t_point.fill_value() 
    285         #value is actually embedded in a multi dimensional list,  
    286         #so we need to extract the actual value from the list 
    287         while type(value) is list: 
    288                 value = value[0] 
    289  
    290         #now check if the value is actually the fill_value rather than  
    291         #a value recorded at the point specified 
    292         log.debug('%s %s' % (value, fill_value)) 
    293         if (2*fill_value) == value: 
    294                 value = "No value found at position: "+str(point[1])+", "+str(point[0]) 
    295         else: 
    296                 value = "Value found at position: "+str(point[1])+", "+str(point[0])+" is: "+str(value) 
    297         # finally return the value 
    298         return value 
    299  
    300  
    301 class CSMLDimension(IwmsDimension): 
    302     """ 
    303     implements IDimension 
    304     @ivar units: The units string. 
    305     @ivar extent: Sequence of extent values. 
    306  
    307     """ 
    308      
    309     def __init__(self, domain, dimname, unit): 
    310         self.units = unit 
    311         self.extent = [] 
    312         #patch to handle current limitations of multiple time dimension scanning in csml.  
    313         if string.lower(self.units)[:10] in ['days_since', 'seconds_si', 'minutes_si', 'hours_sinc','months _sin', 'years_sinc']: 
    314             if type(domain[dimname][0]) is not str   : 
    315                 tunits=self.units.replace('_', ' ') 
    316                 for val in domain[dimname]: 
    317                     csmltime= csml.csmllibs.csmltime.UDtimeToCSMLtime(cdtime.reltime(float(val), tunits).tocomp()) 
    318                     self.extent.append(csmltime) 
    319             else: 
    320                 for val in domain[dimname]: 
    321                     self.extent.append(str(val)) 
    322         else: 
    323             for val in domain[dimname]: 
    324                 self.extent.append(str(val)) 
    325         #for time axis replace units with iso string 
    326         if dimname == 'time': 
    327             self.units='ISO8601' 
    328         
    329 class CSMLLayerSlab(IwmsLayerSlab): 
    330     """ 
    331     Implements LayerSlab 
    332     Represents a particular horizontal slice of a WMS layer. 
    333  
    334     ILayerSlab objects are designed to be convenient to cache. 
    335     They should be pickleable to enable memcached support in the future. 
    336  
    337     @ivar layer: The source ILayer instance. 
    338     @ivar crs: The coordinate reference system. 
    339     @ivar dimValues: A mapping of dimension values of this view. 
    340     @ivar renderOpts: The renderOpts used to create this view. 
    341     @ivar bbox: The bounding box as a 4-tuple. 
    342     """ 
    343      
    344     def __init__(self, netcdf, layer, crs, dimValues, renderOpts, bbox): 
    345         self._netcdf=netcdf 
    346         self.layer = layer 
    347         self.crs = crs 
    348         self.dimValues = dimValues 
    349         self.renderOpts=renderOpts 
    350         self.bbox=bbox 
    351          
    352         #set colour map for ALL images from this slab 
    353         v=self._netcdf(layer.title) 
    354         tvar=v(squeeze=1) 
    355         #get the min and max values to use for the colourmapping. 
    356         #change the fill values to ensure they aren't picked up as false max/mins. 
    357         tvar.missing_value=999999999 
    358         value=tvar.getValue() 
    359         self.minval=min(min(l) for l in value) 
    360         tvar.missing_value=-0.000000001 
    361         value=tvar.getValue() 
    362         self.maxval=max(max(l) for l in value) 
    363  
    364          
    365     def getImage(self, bbox, width, height): 
    366         """ 
    367         Create an image of a sub-bbox of a given size. 
    368  
    369         @ivar bbox: A bbox 4-tuple. 
    370         @ivar width: width in pixels.`   
    371         @ivar height: height in pixels. 
    372         @return: A PIL Image object. 
    373  
    374         """ 
    375  
    376         log.debug('getImage(%s, %s, %s)' % (bbox, width, height)) 
    377          
    378  
    379         cmap=eval(config['colourmap']) # renderOpts is hook for colourmap, for now use config 
    380         grid=Grid(self.layer, self._netcdf, bbox, width, height) 
    381         #how to handle varmin,varmax? ..read array? 
    382         #minval, maxval=genutil.minmax(grid.value) 
    383         #minval=min(min(l) for l in grid.value) 
    384         #maxval=max(max(l) for l in grid.value) 
    385         minval=self.minval 
    386         maxval=self.maxval 
    387         renderer=RGBARenderer(minval, maxval)          
    388         return renderer.renderGrid(grid, bbox, width, height, cmap) 
    389      
    390 class Grid(object): 
    391     """A class encapsulating a simple regularly spaced, rectilinear 
    392     grid.  This is the only type of grid pywms is expected to 
    393     understand and adaptors should be provided to connect to 
    394     underlying implementations such as cdms or csml. 
    395  
    396     @cvar crs: Coordinate reference system 
    397  
    398     @ivar x0: coordinate of the lower left corner. 
    399     @ivar y0: coordinate of the lower left corner. 
    400     @ivar dx: The x grid spacing. 
    401     @ivar dy: The y grid spacing. 
    402     @ivar nx: The number of x grid points. 
    403     @ivar ny: The number of y grid points. 
    404     @ivar value: A masked array of the grid values. 
    405     @ivar ix: The dimension index of longidude in value 
    406     @ivar iy: The dimension index of latitude in value 
    407     @ivar long_name: The name of the field. 
    408     @ivar units: The units of the field. 
    409     """ 
    410     def __init__(self, layer, netcdf, bbox, width, height): 
    411         #we know the axes are called latitude and longitude as the CSML code has written it: 
    412         v=netcdf(layer.title) 
    413         tvar=v(latitude=(bbox[1], bbox[3]), longitude=(bbox[0],bbox[2]),squeeze=1) 
    414         order=tvar.getOrder() 
    415         #array of data 
    416         self.value=tvar.getValue() 
    417         if order == 'xy': 
    418             self.ix=0 
    419             self.iy=1 
    420         else: 
    421             self.ix=1 
    422             self.iy=0 
    423         lat = tvar.getLatitude() 
    424         lon = tvar.getLongitude() 
    425         self.x0=lon[0] 
    426         self.y0=lat[0] 
    427         self.dx=abs(lon[0]-lon[1]) 
    428         self.dy=abs(lat[0]-lat[1]) 
    429          
    430      
    431  
    432          
    433  
    434         self.nx=len(lon) 
    435         self.ny=len(lat) 
    436         self.long_name=tvar.id  #TODO, get long name from feature 
    437         self.units=tvar.units 
Note: See TracChangeset for help on using the changeset viewer.