
.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "gallery/test_decipher_raster.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_gallery_test_decipher_raster.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_gallery_test_decipher_raster.py:

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.

.. GENERATED FROM PYTHON SOURCE LINES 9-188

.. code-block:: Python
   :lineno-start: 10


    import binascii
    import struct

    import pytest
    from sqlalchemy import Column
    from sqlalchemy import Integer
    from sqlalchemy import MetaData
    from sqlalchemy.orm import declarative_base

    from geoalchemy2 import Raster
    from geoalchemy2 import WKTElement

    # Tests imports
    from tests import test_only_with_dialects

    metadata = MetaData()
    Base = declarative_base(metadata=metadata)


    class Ocean(Base):  # type: ignore
        __tablename__ = "ocean"
        id = Column(Integer, primary_key=True)
        rast = Column(Raster)

        def __init__(self, rast):
            self.rast = rast


    def _format_e(endianness, struct_format):
        return _ENDIANNESS[endianness] + struct_format


    def wkbHeader(raw):
        # Function to decipher the WKB header
        # See http://trac.osgeo.org/postgis/browser/trunk/raster/doc/RFC2-WellKnownBinaryFormat

        header = {}

        header["endianness"] = struct.unpack("b", raw[0:1])[0]

        e = header["endianness"]
        header["version"] = struct.unpack(_format_e(e, "H"), raw[1:3])[0]
        header["nbands"] = struct.unpack(_format_e(e, "H"), raw[3:5])[0]
        header["scaleX"] = struct.unpack(_format_e(e, "d"), raw[5:13])[0]
        header["scaleY"] = struct.unpack(_format_e(e, "d"), raw[13:21])[0]
        header["ipX"] = struct.unpack(_format_e(e, "d"), raw[21:29])[0]
        header["ipY"] = struct.unpack(_format_e(e, "d"), raw[29:37])[0]
        header["skewX"] = struct.unpack(_format_e(e, "d"), raw[37:45])[0]
        header["skewY"] = struct.unpack(_format_e(e, "d"), raw[45:53])[0]
        header["srid"] = struct.unpack(_format_e(e, "i"), raw[53:57])[0]
        header["width"] = struct.unpack(_format_e(e, "H"), raw[57:59])[0]
        header["height"] = struct.unpack(_format_e(e, "H"), raw[59:61])[0]

        return header


    def read_band(data, offset, pixtype, height, width, endianness=1):
        ptype, _, psize = _PTYPE[pixtype]
        pix_data = data[offset + 1 : offset + 1 + width * height * psize]
        band = [
            [
                struct.unpack(
                    _format_e(endianness, ptype),
                    pix_data[(i * width + j) * psize : (i * width + j + 1) * psize],
                )[0]
                for j in range(width)
            ]
            for i in range(height)
        ]
        return band


    def read_band_numpy(data, offset, pixtype, height, width, endianness=1):
        import numpy as np  # noqa

        _, dtype, psize = _PTYPE[pixtype]
        dt = np.dtype(dtype)
        dt = dt.newbyteorder(_ENDIANNESS[endianness])
        band = np.frombuffer(data, dtype=dtype, count=height * width, offset=offset + 1)
        band = np.reshape(band, ((height, width)))
        return band


    _PTYPE = {
        0: ["?", "?", 1],
        1: ["B", "B", 1],
        2: ["B", "B", 1],
        3: ["b", "b", 1],
        4: ["B", "B", 1],
        5: ["h", "i2", 2],
        6: ["H", "u2", 2],
        7: ["i", "i4", 4],
        8: ["I", "u4", 4],
        10: ["f", "f4", 4],
        11: ["d", "f8", 8],
    }

    _ENDIANNESS = {
        0: ">",
        1: "<",
    }


    def wkbImage(raster_data, use_numpy=False):
        """Function to decipher the WKB raster data"""
        # Get binary data
        raw = binascii.unhexlify(raster_data)

        # Read header
        h = wkbHeader(bytes(raw))
        e = h["endianness"]

        img = []  # array to store image bands
        offset = 61  # header raw length in bytes
        band_size = h["width"] * h["height"]  # number of pixels in each band

        for _ in range(h["nbands"]):
            # Determine pixtype for this band
            pixtype = struct.unpack(_format_e(e, "b"), raw[offset : offset + 1])[0] - 64

            # Read data with either pure Python or Numpy
            if use_numpy:
                band = read_band_numpy(raw, offset, pixtype, h["height"], h["width"])
            else:
                band = read_band(raw, offset, pixtype, h["height"], h["width"])

            # Store the result
            img.append(band)
            offset = offset + 2 + band_size

        return img


    @test_only_with_dialects("postgresql")
    class TestDecipherRaster:
        @pytest.mark.parametrize(
            "pixel_type",
            [
                "1BB",
                "2BUI",
                "4BUI",
                "8BSI",
                "8BUI",
                "16BSI",
                "16BUI",
                "32BSI",
                "32BUI",
                "32BF",
                "64BF",
            ],
        )
        def test_decipher_raster(self, pixel_type, session, conn):
            """Create a raster and decipher it"""
            metadata.drop_all(conn, checkfirst=True)
            metadata.create_all(conn)

            # Create a new raster
            polygon = WKTElement("POLYGON((0 0,1 1,0 1,0 0))", srid=4326)
            o = Ocean(polygon.ST_AsRaster(5, 6, pixel_type))
            session.add(o)
            session.flush()

            # Decipher data from each raster
            image = wkbImage(o.rast.data)

            # Define expected result
            expected = [
                [0, 1, 1, 1, 1],
                [1, 1, 1, 1, 1],
                [0, 1, 1, 1, 0],
                [0, 1, 1, 0, 0],
                [0, 1, 0, 0, 0],
                [0, 0, 0, 0, 0],
            ]

            # Check results
            band = image[0]
            assert band == expected


.. _sphx_glr_download_gallery_test_decipher_raster.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: test_decipher_raster.ipynb <test_decipher_raster.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: test_decipher_raster.py <test_decipher_raster.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: test_decipher_raster.zip <test_decipher_raster.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_
