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

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@3745
Revision 3745, 6.8 KB checked in by spascoe, 12 years ago (diff)

More GDAL code.

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