Note
Click here to download the full example code
Decipher RasterΒΆ
The RasterElement objects store the Raster data in WKB form. When using rasters it is usually better to convert them into TIFF, PNG, JPEG or whatever. Nevertheless, it is possible to decipher the WKB to get a 2D list of values. This example uses SQLAlchemy ORM queries.
10 import binascii
11 import struct
12
13 import pytest
14 from sqlalchemy import Column
15 from sqlalchemy import create_engine
16 from sqlalchemy import Integer
17 from sqlalchemy import MetaData
18 from sqlalchemy.ext.declarative import declarative_base
19 from sqlalchemy.orm import sessionmaker
20
21 from geoalchemy2 import Raster, WKTElement
22
23
24 engine = create_engine('postgresql://gis:gis@localhost/gis', echo=False)
25 metadata = MetaData(engine)
26 Base = declarative_base(metadata=metadata)
27
28 session = sessionmaker(bind=engine)()
29
30
31 class Ocean(Base):
32 __tablename__ = 'ocean'
33 id = Column(Integer, primary_key=True)
34 rast = Column(Raster)
35
36 def __init__(self, rast):
37 self.rast = rast
38
39
40 def _format_e(endianess, struct_format):
41 return _ENDIANESS[endianess] + struct_format
42
43
44 def wkbHeader(raw):
45 # Function to decipher the WKB header
46 # See http://trac.osgeo.org/postgis/browser/trunk/raster/doc/RFC2-WellKnownBinaryFormat
47
48 header = {}
49
50 header['endianess'] = struct.unpack('b', raw[0:1])[0]
51
52 e = header['endianess']
53 header['version'] = struct.unpack(_format_e(e, 'H'), raw[1:3])[0]
54 header['nbands'] = struct.unpack(_format_e(e, 'H'), raw[3:5])[0]
55 header['scaleX'] = struct.unpack(_format_e(e, 'd'), raw[5:13])[0]
56 header['scaleY'] = struct.unpack(_format_e(e, 'd'), raw[13:21])[0]
57 header['ipX'] = struct.unpack(_format_e(e, 'd'), raw[21:29])[0]
58 header['ipY'] = struct.unpack(_format_e(e, 'd'), raw[29:37])[0]
59 header['skewX'] = struct.unpack(_format_e(e, 'd'), raw[37:45])[0]
60 header['skewY'] = struct.unpack(_format_e(e, 'd'), raw[45:53])[0]
61 header['srid'] = struct.unpack(_format_e(e, 'i'), raw[53:57])[0]
62 header['width'] = struct.unpack(_format_e(e, 'H'), raw[57:59])[0]
63 header['height'] = struct.unpack(_format_e(e, 'H'), raw[59:61])[0]
64
65 return header
66
67
68 def read_band(data, offset, pixtype, height, width, endianess=1):
69 ptype, _, psize = _PTYPE[pixtype]
70 pix_data = data[offset + 1: offset + 1 + width * height * psize]
71 band = [
72 [
73 struct.unpack(_format_e(endianess, ptype), pix_data[
74 (i * width + j) * psize: (i * width + j + 1) * psize
75 ])[0]
76 for j in range(width)
77 ]
78 for i in range(height)
79 ]
80 return band
81
82
83 def read_band_numpy(data, offset, pixtype, height, width, endianess=1):
84 import numpy as np # noqa
85 _, dtype, psize = _PTYPE[pixtype]
86 dt = np.dtype(dtype)
87 dt = dt.newbyteorder(_ENDIANESS[endianess])
88 band = np.frombuffer(data, dtype=dtype,
89 count=height * width, offset=offset + 1)
90 band = (np.reshape(band, ((height, width))))
91 return band
92
93
94 _PTYPE = {
95 0: ['?', '?', 1],
96 1: ['B', 'B', 1],
97 2: ['B', 'B', 1],
98 3: ['b', 'b', 1],
99 4: ['B', 'B', 1],
100 5: ['h', 'i2', 2],
101 6: ['H', 'u2', 2],
102 7: ['i', 'i4', 4],
103 8: ['I', 'u4', 4],
104 10: ['f', 'f4', 4],
105 11: ['d', 'f8', 8],
106 }
107
108 _ENDIANESS = {
109 0: '>',
110 1: '<',
111 }
112
113
114 def wkbImage(raster_data, use_numpy=False):
115 """Function to decipher the WKB raster data"""
116
117 # Get binary data
118 raw = binascii.unhexlify(raster_data)
119
120 # Read header
121 h = wkbHeader(bytes(raw))
122 e = h["endianess"]
123
124 img = [] # array to store image bands
125 offset = 61 # header raw length in bytes
126 band_size = h['width'] * h['height'] # number of pixels in each band
127
128 for i in range(h['nbands']):
129 # Determine pixtype for this band
130 pixtype = struct.unpack(_format_e(e, 'b'), raw[offset: offset + 1])[0] - 64
131
132 # Read data with either pure Python or Numpy
133 if use_numpy:
134 band = read_band_numpy(
135 raw, offset, pixtype, h['height'], h['width'])
136 else:
137 band = read_band(
138 raw, offset, pixtype, h['height'], h['width'])
139
140 # Store the result
141 img.append(band)
142 offset = offset + 2 + band_size
143
144 return img
145
146
147 class TestDecipherRaster():
148
149 def setup(self):
150 metadata.drop_all(checkfirst=True)
151 metadata.create_all()
152
153 def teardown(self):
154 session.rollback()
155 metadata.drop_all()
156
157 @pytest.mark.parametrize("pixel_type", [
158 '1BB',
159 '2BUI',
160 '4BUI',
161 '8BSI',
162 '8BUI',
163 '16BSI',
164 '16BUI',
165 '32BSI',
166 '32BUI',
167 '32BF',
168 '64BF'
169 ])
170 def test_decipher_raster(self, pixel_type):
171 """Create a raster and decipher it"""
172
173 # Create a new raster
174 polygon = WKTElement('POLYGON((0 0,1 1,0 1,0 0))', srid=4326)
175 o = Ocean(polygon.ST_AsRaster(5, 6, pixel_type))
176 session.add(o)
177 session.flush()
178
179 # Decipher data from each raster
180 image = wkbImage(o.rast.data)
181
182 # Define expected result
183 expected = [
184 [0, 1, 1, 1, 1],
185 [1, 1, 1, 1, 1],
186 [0, 1, 1, 1, 0],
187 [0, 1, 1, 0, 0],
188 [0, 1, 0, 0, 0],
189 [0, 0, 0, 0, 0]
190 ]
191
192 # Check results
193 band = image[0]
194 assert band == expected
Total running time of the script: ( 0 minutes 0.000 seconds)