source: TI05-delivery/ows_framework/branches/ows_framework-refactor/ows_common/ows_common/test/test_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/test/test_wms_gdal.py@3956
Revision 3745, 3.1 KB checked in by spascoe, 12 years ago (diff)

More GDAL code.

Line 
1"""
2Test ows_common.service.wms_gdal.
3
4"""
5
6from ows_common.service.wms_gdal import *
7from pkg_resources import resource_filename
8from osgeo import gdal, osr
9from osgeo.gdalconst import *
10
11import unittest
12
13
14
15def geoTransform(x, y, T):
16    xg = T[0] + float(x)*T[1] + float(y)*T[2]
17    yg = T[3] + float(x)*T[4] + float(y)*T[5]
18
19    return (xg, yg)
20
21class StaticDataSource(IGDALDataSource):
22    def __init__(self, filename):
23        self._ds = gdal.Open(filename, GA_ReadOnly)
24        if self._ds is None:
25            raise ValueError('Dataset open failed')
26
27        self.title = 'Static data source %s' % filename
28        self.abstract = None
29        self.dimensions = {}
30        self.units = None
31        # Deduce CRS from the WKT
32        sr = osr.SpatialReference()
33        sr.ImportFromWkt(self._ds.GetProjection())
34
35        if sr.GetAuthorityName('PROJCS'):
36            self.crs = '%s:%s' % (sr.GetAuthorityName('PROJCS'),
37                                  sr.GetAuthorityCode('PROJCS'))
38        else:
39            self.crs = '%s:%s' % (sr.GetAuthorityName('GEOGCS'),
40                                  sr.GetAuthorityCode('GEOGCS'))
41
42    def getWKT(self):
43        return self._ds.GetProjection()
44
45    def getBBox(self):
46        # Calculate the bounding box from the geoTransform.
47        T = self._ds.GetGeoTransform()
48        w, h = self._ds.RasterXSize, self._ds.RasterYSize
49
50        # the GeoTransform is relative to the top-left corner of pixel 0,0
51        # therefore to encompass all of the pixels we add 1 to w and h.
52        llx, lly = geoTransform(0, h+1, T)
53        urx, ury = geoTransform(w+1, 0, T)
54        bbox = (llx, lly, urx, ury)
55
56        return bbox
57
58    def getDataset(self, dimValues=None, renderOpts={}):
59        return self._ds
60
61class TestGDAL(unittest.TestSuite):
62    def setUp(self):
63        filename = resource_filename('ows_common.test', 'data/obs_gdal.png')
64        self.globalDS = StaticDataSource(filename)
65        fn2 = resource_filename('ows_common.test', 'data/ukcip_rainfall.png')
66        self.ukDS = StaticDataSource(fn2)
67
68    def setupWarp(self):
69        layer = GDALLayer(self.ukDS)
70        sr = osr.SpatialReference()
71        sr.ImportFromEPSG(4326)
72        layer.warpCRS['EPSG:4326'] = sr.ExportToWkt()
73
74        return layer
75
76    def test1(self):
77        layer = GDALLayer(self.globalDS)
78        assert layer.crss == ['EPSG:4326']
79
80    def test2(self):
81        layer = GDALLayer(self.globalDS)
82        slab = layer.getSlab('EPSG:4326')
83        assert slab.bbox == (-180.0, -90.0, 180.0, 90.0)
84
85    def test3(self):
86        layer = GDALLayer(self.globalDS)
87        slab = layer.getSlab('EPSG:4326')
88        img = slab.getImage((-10.0, 45.0, 5.0, 60.0), 100, 100)
89        img.save('out1.png')
90
91    def test4(self):
92        layer = self.setupWarp()
93        print layer.crss
94        assert layer.crss == ['EPSG:27700', 'EPSG:4326']
95
96    def test5(self):
97        layer = self.setupWarp()
98        slab = layer.getSlab('EPSG:4326')
99
100        # Create a copy of the warped image
101        ds = slab._data
102        import pdb; pdb.set_trace()
103
104        img = slab.getImage((-3.0, 54.0, 3.0, 60.0), 100, 100)
105        img.save('out2.png')
Note: See TracBrowser for help on using the repository browser.