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

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@3634
Revision 3634, 6.4 KB checked in by spascoe, 13 years ago (diff)

Test for GDAL support. Not completely working yet.

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