RasterField for PostGIS backends

Thanks to Tim Graham and Claude Paroz for the reviews and patches.

django/contrib/gis/db/backends/base/features.py

@@ -37,6 +37,9 @@ class BaseSpatialFeatures(object):
     supports_distances_lookups = True
     supports_left_right_lookups = False
 
+    # Does the database have raster support?
+    supports_raster = False
+
     @property
     def supports_bbcontains_lookup(self):
         return 'bbcontains' in self.connection.ops.gis_operators

django/contrib/gis/db/backends/postgis/const.py

@@ -0,0 +1,43 @@
+"""
+PostGIS to GDAL conversion constant definitions
+"""
+# Lookup to convert pixel type values from GDAL to PostGIS
+GDAL_TO_POSTGIS = [None, 4, 6, 5, 8, 7, 10, 11, None, None, None, None]
+
+# Lookup to convert pixel type values from PostGIS to GDAL
+POSTGIS_TO_GDAL = [1, 1, 1, 3, 1, 3, 2, 5, 4, None, 6, 7, None, None]
+
+# Struct pack structure for raster header, the raster header has the
+# following structure:
+#
+# Endianness, PostGIS raster version, number of bands, scale, origin,
+# skew, srid, width, and height.
+#
+# Scale, origin, and skew have x and y values. PostGIS currently uses
+# a fixed endianness (1) and there is only one version (0).
+POSTGIS_HEADER_STRUCTURE = 'B H H d d d d d d i H H'
+
+# Lookup values to convert GDAL pixel types to struct characters. This is
+# used to pack and unpack the pixel values of PostGIS raster bands.
+GDAL_TO_STRUCT = [
+    None, 'B', 'H', 'h', 'L', 'l', 'f', 'd',
+    None, None, None, None,
+]
+
+# Size of the packed value in bytes for different numerical types.
+# This is needed to cut chunks of band data out of PostGIS raster strings
+# when decomposing them into GDALRasters.
+# See https://docs.python.org/3/library/struct.html#format-characters
+STRUCT_SIZE = {
+    'b': 1,  # Signed char
+    'B': 1,  # Unsigned char
+    '?': 1,  # _Bool
+    'h': 2,  # Short
+    'H': 2,  # Unsigned short
+    'i': 4,  # Integer
+    'I': 4,  # Unsigned Integer
+    'l': 4,  # Long
+    'L': 4,  # Unsigned Long
+    'f': 4,  # Float
+    'd': 8,  # Double
+}

django/contrib/gis/db/backends/postgis/features.py

@@ -7,3 +7,4 @@ class DatabaseFeatures(BaseSpatialFeatures, Psycopg2DatabaseFeatures):
     supports_3d_storage = True
     supports_3d_functions = True
     supports_left_right_lookups = True
+    supports_raster = True

django/contrib/gis/db/backends/postgis/introspection.py

@@ -20,6 +20,26 @@ class PostGISIntrospection(DatabaseIntrospection):
         'raster_overviews',
     ]
 
+    # Overridden from parent to include raster indices in retrieval.
+    # Raster indices have pg_index.indkey value 0 because they are an
+    # expression over the raster column through the ST_ConvexHull function.
+    # So the default query has to be adapted to include raster indices.
+    _get_indexes_query = """
+        SELECT DISTINCT attr.attname, idx.indkey, idx.indisunique, idx.indisprimary
+        FROM pg_catalog.pg_class c, pg_catalog.pg_class c2,
+            pg_catalog.pg_index idx, pg_catalog.pg_attribute attr
+        LEFT JOIN pg_catalog.pg_type t ON t.oid = attr.atttypid
+        WHERE
+            c.oid = idx.indrelid
+            AND idx.indexrelid = c2.oid
+            AND attr.attrelid = c.oid
+            AND (
+                attr.attnum = idx.indkey[0] OR
+                (t.typname LIKE 'raster' AND idx.indkey = '0')
+            )
+            AND attr.attnum > 0
+            AND c.relname = %s"""
+
     def get_postgis_types(self):
         """
         Returns a dictionary with keys that are the PostgreSQL object

django/contrib/gis/db/backends/postgis/operations.py

@@ -4,6 +4,9 @@
 from django.contrib.gis.db.backends.base.operations import \
     BaseSpatialOperations
 from django.contrib.gis.db.backends.postgis.adapter import PostGISAdapter
+from django.contrib.gis.db.backends.postgis.pgraster import (
+    from_pgraster, to_pgraster,
+)
 from django.contrib.gis.db.backends.utils import SpatialOperator
 from django.contrib.gis.geometry.backend import Geometry
 from django.contrib.gis.measure import Distance
@@ -14,6 +17,7 @@
 from django.utils.functional import cached_property
 
 from .models import PostGISGeometryColumns, PostGISSpatialRefSys
+from .pgraster import get_pgraster_srid
 
 
 class PostGISOperator(SpatialOperator):
@@ -205,12 +209,11 @@ def convert_geom(self, hex, geo_field):
 
     def geo_db_type(self, f):
         """
-        Return the database field type for the given geometry field.
-        Typically this is `None` because geometry columns are added via
-        the `AddGeometryColumn` stored procedure, unless the field
-        has been specified to be of geography type instead.
+        Return the database field type for the given spatial field.
         """
-        if f.geography:
+        if f.geom_type == 'RASTER':
+            return 'raster'
+        elif f.geography:
             if f.srid != 4326:
                 raise NotImplementedError('PostGIS only supports geography columns with an SRID of 4326.')
 
@@ -272,10 +275,21 @@ def get_geom_placeholder(self, f, value, compiler):
         SRID of the field.  Specifically, this routine will substitute in the
         ST_Transform() function call.
         """
-        if value is None or value.srid == f.srid:
+        # Get the srid for this object
+        if value is None:
+            value_srid = None
+        elif f.geom_type == 'RASTER':
+            value_srid = get_pgraster_srid(value)
+        else:
+            value_srid = value.srid
+
+        # Adding Transform() to the SQL placeholder if the value srid
+        # is not equal to the field srid.
+        if value_srid is None or value_srid == f.srid:
             placeholder = '%s'
+        elif f.geom_type == 'RASTER':
+            placeholder = '%s((%%s)::raster, %s)' % (self.transform, f.srid)
         else:
-            # Adding Transform() to the SQL placeholder.
             placeholder = '%s(%%s, %s)' % (self.transform, f.srid)
 
         if hasattr(value, 'as_sql'):
@@ -359,3 +373,11 @@ def geometry_columns(self):
 
     def spatial_ref_sys(self):
         return PostGISSpatialRefSys
+
+    # Methods to convert between PostGIS rasters and dicts that are
+    # readable by GDALRaster.
+    def parse_raster(self, value):
+        return from_pgraster(value)
+
+    def deconstruct_raster(self, value):
+        return to_pgraster(value)

django/contrib/gis/db/backends/postgis/pgraster.py

@@ -0,0 +1,161 @@
+import binascii
+import struct
+
+from django.forms import ValidationError
+
+from .const import (
+    GDAL_TO_POSTGIS, GDAL_TO_STRUCT, POSTGIS_HEADER_STRUCTURE, POSTGIS_TO_GDAL,
+    STRUCT_SIZE,
+)
+
+
+def pack(structure, data):
+    """
+    Pack data into hex string with little endian format.
+    """
+    return binascii.hexlify(struct.pack('<' + structure, *data)).upper()
+
+
+def unpack(structure, data):
+    """
+    Unpack little endian hexlified binary string into a list.
+    """
+    return struct.unpack('<' + structure, binascii.unhexlify(data))
+
+
+def chunk(data, index):
+    """
+    Split a string into two parts at the input index.
+    """
+    return data[:index], data[index:]
+
+
+def get_pgraster_srid(data):
+    """
+    Extract the SRID from a PostGIS raster string.
+    """
+    if data is None:
+        return
+    # The positional arguments here extract the hex-encoded srid from the
+    # header of the PostGIS raster string. This can be understood through
+    # the POSTGIS_HEADER_STRUCTURE constant definition in the const module.
+    return unpack('i', data[106:114])[0]
+
+
+def from_pgraster(data):
+    """
+    Convert a PostGIS HEX String into a dictionary.
+    """
+    if data is None:
+        return
+
+    # Split raster header from data
+    header, data = chunk(data, 122)
+    header = unpack(POSTGIS_HEADER_STRUCTURE, header)
+
+    # Parse band data
+    bands = []
+    pixeltypes = []
+    while data:
+        # Get pixel type for this band
+        pixeltype, data = chunk(data, 2)
+        pixeltype = unpack('B', pixeltype)[0]
+
+        # Subtract nodata byte from band nodata value if exists
+        has_nodata = pixeltype >= 64
+        if has_nodata:
+            pixeltype -= 64
+
+        # Convert datatype from PostGIS to GDAL & get pack type and size
+        pixeltype = POSTGIS_TO_GDAL[pixeltype]
+        pack_type = GDAL_TO_STRUCT[pixeltype]
+        pack_size = 2 * STRUCT_SIZE[pack_type]
+
+        # Parse band nodata value. The nodata value is part of the
+        # PGRaster string even if the nodata flag is True, so it always
+        # has to be chunked off the data string.
+        nodata, data = chunk(data, pack_size)
+        nodata = unpack(pack_type, nodata)[0]
+
+        # Chunk and unpack band data (pack size times nr of pixels)
+        band, data = chunk(data, pack_size * header[10] * header[11])
+        band_result = {'data': binascii.unhexlify(band)}
+
+        # If the nodata flag is True, set the nodata value.
+        if has_nodata:
+            band_result['nodata_value'] = nodata
+
+        # Append band data to band list
+        bands.append(band_result)
+
+        # Store pixeltype of this band in pixeltypes array
+        pixeltypes.append(pixeltype)
+
+    # Check that all bands have the same pixeltype.
+    # This is required by GDAL. PostGIS rasters could have different pixeltypes
+    # for bands of the same raster.
+    if len(set(pixeltypes)) != 1:
+        raise ValidationError("Band pixeltypes are not all equal.")
+
+    return {
+        'srid': int(header[9]),
+        'width': header[10], 'height': header[11],
+        'datatype': pixeltypes[0],
+        'origin': (header[5], header[6]),
+        'scale': (header[3], header[4]),
+        'skew': (header[7], header[8]),
+        'bands': bands,
+    }
+
+
+def to_pgraster(rast):
+    """
+    Convert a GDALRaster into PostGIS Raster format.
+    """
+    # Return if the raster is null
+    if rast is None or rast == '':
+        return
+
+    # Prepare the raster header data as a tuple. The first two numbers are
+    # the endianness and the PostGIS Raster Version, both are fixed by
+    # PostGIS at the moment.
+    rasterheader = (
+        1, 0, len(rast.bands), rast.scale.x, rast.scale.y,
+        rast.origin.x, rast.origin.y, rast.skew.x, rast.skew.y,
+        rast.srs.srid, rast.width, rast.height
+    )
+
+    # Hexlify raster header
+    result = pack(POSTGIS_HEADER_STRUCTURE, rasterheader)
+
+    for band in rast.bands:
+        # The PostGIS raster band header has exactly two elements, a 8BUI byte
+        # and the nodata value.
+        #
+        # The 8BUI stores both the PostGIS pixel data type and a nodata flag.
+        # It is composed as the datatype integer plus 64 as a flag for existing
+        # nodata values:
+        # 8BUI_VALUE = PG_PIXEL_TYPE (0-11) + FLAG (0 or 64)
+        #
+        # For example, if the byte value is 71, then the datatype is
+        # 71-64 = 7 (32BSI) and the nodata value is True.
+        structure = 'B' + GDAL_TO_STRUCT[band.datatype()]
+
+        # Get band pixel type in PostGIS notation
+        pixeltype = GDAL_TO_POSTGIS[band.datatype()]
+
+        # Set the nodata flag
+        if band.nodata_value is not None:
+            pixeltype += 64
+
+        # Pack band header
+        bandheader = pack(structure, (pixeltype, band.nodata_value or 0))
+
+        # Hexlify band data
+        band_data_hex = binascii.hexlify(band.data(as_memoryview=True)).upper()
+
+        # Add packed header and band data to result
+        result += bandheader + band_data_hex
+
+    # Cast raster to string before passing it to the DB
+    return result.decode()

django/contrib/gis/db/backends/postgis/schema.py

@@ -4,6 +4,7 @@
 class PostGISSchemaEditor(DatabaseSchemaEditor):
     geom_index_type = 'GIST'
     geom_index_ops_nd = 'GIST_GEOMETRY_OPS_ND'
+    rast_index_wrapper = 'ST_ConvexHull(%s)'
 
     sql_add_spatial_index = "CREATE INDEX %(index)s ON %(table)s USING %(index_type)s (%(column)s %(ops)s)"
     sql_clear_geometry_columns = "DELETE FROM geometry_columns WHERE f_table_name = %(table)s"
@@ -16,17 +17,22 @@ def geo_quote_name(self, name):
         return self.connection.ops.geo_quote_name(name)
 
     def column_sql(self, model, field, include_default=False):
-        from django.contrib.gis.db.models.fields import GeometryField
-        if not isinstance(field, GeometryField):
+        from django.contrib.gis.db.models.fields import BaseSpatialField
+        if not isinstance(field, BaseSpatialField):
             return super(PostGISSchemaEditor, self).column_sql(model, field, include_default)
 
         column_sql = super(PostGISSchemaEditor, self).column_sql(model, field, include_default)
 
         if field.spatial_index:
             # Spatial indexes created the same way for both Geometry and
             # Geography columns.
-
-            if field.geography:
+            field_column = self.quote_name(field.column)
+            if field.geom_type == 'RASTER':
+                # For raster fields, wrap index creation SQL statement with ST_ConvexHull.
+                # Indexes on raster columns are based on the convex hull of the raster.
+                field_column = self.rast_index_wrapper % field_column
+                index_ops = ''
+            elif field.geography:
                 index_ops = ''
             else:
                 # Use either "nd" ops  which are fast on multidimensional cases
@@ -39,7 +45,7 @@ def column_sql(self, model, field, include_default=False):
                 self.sql_add_spatial_index % {
                     "index": self.quote_name('%s_%s_id' % (model._meta.db_table, field.column)),
                     "table": self.quote_name(model._meta.db_table),
-                    "column": self.quote_name(field.column),
+                    "column": field_column,
                     "index_type": self.geom_index_type,
                     "ops": index_ops,
                 }

django/contrib/gis/db/models/__init__.py

@@ -11,4 +11,4 @@
 from django.contrib.gis.db.models.fields import (  # NOQA
     GeometryField, PointField, LineStringField, PolygonField,
     MultiPointField, MultiLineStringField, MultiPolygonField,
-    GeometryCollectionField)
+    GeometryCollectionField, RasterField)