source: TI05-delivery/ows_framework/branches/ows_framework-refactor/ows_common/ows_common/service/wms_gdal.py @ 3621

Subversion URL: http://proj.badc.rl.ac.uk/svn/ndg/TI05-delivery/ows_framework/branches/ows_framework-refactor/ows_common/ows_common/service/wms_gdal.py@3621
Revision 3621, 5.8 KB checked in by spascoe, 14 years ago (diff)

Work in progress on GDAL warping support.

Line 
1"""
2An implementation of ows_common.service.wms_iface that uses GDAL to support
3warping between multiple coordinate reference systems.  This implementation
4relies on a further interface, IGDALDataSource, to provide the data.
5
6"""
7
8from ows_common.service.wms_iface import *
9from ows_common.bbox_util import geoToPixel
10
11import logging
12log = logging.getLogger(__name__)
13
14class IGDALDataSource(object):
15    """
16    This interface is very similar to ILayer except that it returns GDAL
17    datasets rather than PIL images.  It also doesn't try to handle multiple
18    CRSs as this is handled by GDALLayer. 
19
20    @ivar title: The layer title.  As seen in the Capabilities document.
21    @ivar abstract:  Abstract as seen in the Capabilities document.
22    @ivar dimensions: A mapping of dimension names to IDimension objects.
23    @ivar units: A string describing the units.
24    @ivar crs: The CRS that GDAL datasets will be returned in by
25        self.getDataset()
26
27    @todo: Legend plotting needs support but should probably be done in a
28        seperate interface.
29
30    """
31    dimensions = NotImplemented
32
33    def getWKT(self):
34        """
35        Because mapping between CRS codes and WKT format can be flakey in GDAL
36        this function allows the problem to be solved on a case-by-case basis.
37       
38        @return: the description of self.crs in GDAL well known text format.
39
40        """
41       
42    def getBBox(self):
43        """
44        @return: the bounding box (llx, lly, urx, ury) in self.crs.
45
46        """
47
48    def getDataset(self, dimValues=None, renderOpts={}):
49        """
50        Create the equivilent of ILayerSlab as a GDAL dataset.  The dataset
51        could have 1,3 or 4 bands representing PIL modes 'L', 'RGB' or 'RGBA'.
52       
53        @param dimValues: A mapping of dimension names to dimension values
54            as specified in the IDimension.extent.
55        @param renderOpts: A generic mapping object for passing rendering
56            options.
57        @return: A GDAL Dataset object for this horizontal slice.
58           
59        """
60
61
62class GDALLayer(ILayer):
63    """
64    This implementation of ILayer can warp images from a source CRS to
65    various other CRSs.
66
67    @ivar sourceCRS: The CRS of the data source.
68    @ivar warpCRS: A mapping of CRS itentifiers to WKT descriptions of
69       CRSs that are supported for this ILayer via warping.
70
71    """
72
73    def __init__(self, dataSource):
74        """
75        @param dataSource: A IGDALDataSource implementation.
76
77        """
78        self._ds = dataSource
79        self.warpCRS = {}
80        self.sourceCRS = dataSource.crs
81
82        self.title = dataSource.title
83        self.abstract = dataSource.abstract
84        self.dimensions = dataSource.dimensions
85        self.units = dataSource.units
86        #!NOTE: self.crss is implemented as property
87
88    def _getCRSs(self):
89        return [self._ds.crs] + self._warpCRS.keys()
90    crss = property(_getCRSs)
91
92    def getBBox(self, crs):
93        src_bb = self._ds.getBBox()
94        if self.crs == self._ds.crs:
95            return src_bb
96
97        sr_src = osr.SpatialReference(self._ds.getWKT())
98        sr_dst = osr.SpatialReference(self.warpCRS[crs])
99
100        ct = osr.CoordinateTransformation(sr_srs, sr_dst)
101
102        llx, lly = ct.TransformPoint(*src_bb[:2])[:2]
103        urx, ury = ct.TransformPoint(*src_bb[2:])[:2]
104        return (llx, lly, urx, ury)
105
106    def getSlab(self, crs, dimValues=None, renderOpts={}):
107        return GDALLayerSlab(self, crs, dimValues=dimValues,
108                             renderOpts=renderOpts)
109
110    def getCacheKey(self, crs, dimValues=None, renderOpts={}):
111        """
112        A fairly sane cache key generation algorithm.
113
114        """
115        if dimValues is None:
116            x = None
117        else:
118            x = dimValues.items()
119        x.sort()
120        y = renderOpts.items(); y.sort()
121
122        return str((x, y))
123
124class GDALLayerSlab(ILayerSlab):
125    def __init__(self, layer, crs, dimValues=None, renderOpts={}):
126        self.layer = layer
127        self.crs = crs
128        self.dimValues = dimValues
129        self.rendOpts = renderOpts
130        self.bbox = layer.getBBOX(crs)
131       
132        if crs == layer.sourceCRS:
133            self._data = layer._ds.getDataset()
134        else:
135            self._data = warpDataset(layer._ds, layer.warpCRS[crs])
136
137    def getImage(self, bbox, width, height):       
138        # Calculate the pixel coordinates of bbox within self.bbox
139        xoff, yoff = geoToPixel(bbox[0], bbox[1], self.bbox, width, height)
140        xsize, ysize = geoToPixel(bbox[2], bbox[3], self.bbox, width, height)
141
142        return datasetToImage(self._data, xoff, yoff, xsize, ysize)
143
144
145#-----------------------------------------------------------------------------
146# Utility functions
147
148def datasetToImage(ds, xoff, yoff, xsize, ysize):
149    """
150    Convert a GDAL dataset into a PIL image with cropping.
151
152    """
153    ds = self._ds
154    bandImages = []
155    for iband in range(1, ds.RasterCount+1):
156        band = ds.GetRasterBand(iband)
157        bandImages.append(Image.fromstring('L', (x, y),
158                                           band.ReadRaster(xoff, yoff,
159                                                           xsize, ysize)))
160    return Image.merge('RGBA', bandImages)
161
162def warpDataset(dataSource, wkt, driverName='MEM', datasetName=''):
163    """
164    Warp a GDAL dataset from one CRS to another.
165
166    @param dataset: An object implementing IGDALDataSource.
167    @param wkt: The Well Known Text string of the destination CRS
168    @param driverName: The GDAL driver to use for the new dataset.
169        This driver must support the Create() method.
170    @param datasetName: The name to give the dataset (i.e. the filename if
171        a file-based driver)
172    @param return: A GDAL dataset in the new CRS.
173
174    """
175    ds = dataSource.getDataset()
176    dr = gdal.GetDriverByName(driverName)
177    # ...
178    #dsOut = dr.Create(datasetName, ds.
179   
180
181
Note: See TracBrowser for help on using the repository browser.