RasterField for PostGIS backends

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

@@ -35,6 +35,9 @@ class BaseSpatialFeatures(object):
     supports_distances_lookups = True
     supports_left_right_lookups = False
 
+    # Does the database have raster suport?
+    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 endiannes (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/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
@@ -222,6 +225,9 @@ def geo_db_type(self, f):
         the `AddGeometryColumn` stored procedure, unless the field
         has been specified to be of geography type instead.
         """
+        if f.geom_type == 'RASTER':
+            return 'raster'
+
         if f.geography:
             if f.srid != 4326:
                 raise NotImplementedError('PostGIS only supports geography columns with an SRID of 4326.')
@@ -373,3 +379,9 @@ def geometry_columns(self):
 
     def spatial_ref_sys(self):
         return PostGISSpatialRefSys
+
+    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,147 @@
+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):
+    """
+    Packs data into hex string with little endian format.
+    """
+    return binascii.hexlify(struct.pack('<' + structure, *data)).upper()
+
+
+def unpack(structure, data):
+    """
+    Unpacks little endian hexlified binary string into python list.
+    """
+    return struct.unpack('<' + structure, binascii.unhexlify(data))
+
+
+def chunk(data, index):
+    """
+    Splits a string into two parts at the input index.
+    """
+    return data[:index], data[index:]
+
+
+def band_to_hex(band):
+    """
+    Returns a GDALBand's pixel values as PGRaster Band hex string.
+    """
+    return binascii.hexlify(band.data(as_memoryview=True)).upper()
+
+
+def from_pgraster(data):
+    """
+    Converts a PostGIS HEX String into a python dictionary.
+    """
+    # Abort if input data is null
+    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, even if it is ignored by the nodata flag
+        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 has_nodata:
+            band_result['nodata_value'] = nodata
+
+        bands.append(band_result)
+        pixeltypes.append(pixeltype)
+
+    # Check that all bands have the same pixeltype, this is required by GDAL
+    if len(set(pixeltypes)) != 1:
+        raise ValidationError("Band pixeltypes are not all equal.")
+
+    # Process raster header
+    return {
+        'srid': int(header[8]),
+        '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):
+    """
+    Converts a GDALRaster into PostGIS Raster format.
+    """
+    # Abort if raster is null
+    if rast is None:
+        return
+
+    # Prepare the raster header data as 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 WKB band header has 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. As a formula one could write:
+        # 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))
+
+        # Add packed header and band data to result string
+        result += bandheader + band_to_hex(band)
+
+    # Cast raster to string before passing it to the DB
+    result = result.decode()
+
+    return result

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

@@ -5,6 +5,7 @@ class PostGISSchemaEditor(DatabaseSchemaEditor):
     geom_index_type = 'GIST'
     geom_index_ops = 'GIST_GEOMETRY_OPS'
     geom_index_ops_nd = 'GIST_GEOMETRY_OPS_ND'
+    rast_index_wrap = 'ST_ConvexHull(%s)'
 
     sql_add_geometry_column = "SELECT AddGeometryColumn(%(table)s, %(column)s, %(srid)s, %(geom_type)s, %(dim)s)"
     sql_drop_geometry_column = "SELECT DropGeometryColumn(%(table)s, %(column)s)"
@@ -20,12 +21,12 @@ 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 GeometryField, RasterField
+        if not isinstance(field, (GeometryField, RasterField)):
             return super(PostGISSchemaEditor, self).column_sql(model, field, include_default)
 
-        if field.geography or self.connection.ops.geometry:
-            # Geography and Geometry (PostGIS 2.0+) columns are
+        if isinstance(field, RasterField) or field.geography or self.connection.ops.geometry:
+            # Geography, Geometry and Raster (PostGIS 2.0+) columns are
             # created normally.
             column_sql = super(PostGISSchemaEditor, self).column_sql(model, field, include_default)
         else:
@@ -57,7 +58,11 @@ def column_sql(self, model, field, include_default=False):
             # we use GIST_GEOMETRY_OPS, on 2.0 we use either "nd" ops
             # which are fast on multidimensional cases, or just plain
             # gist index for the 2d case.
-            if field.geography:
+            column = self.quote_name(field.column)
+            if isinstance(field, RasterField):
+                column = self.rast_index_wrap % field.column
+                index_ops = ''
+            elif field.geography:
                 index_ops = ''
             elif self.connection.ops.geometry:
                 if field.dim > 2:
@@ -66,11 +71,12 @@ def column_sql(self, model, field, include_default=False):
                     index_ops = ''
             else:
                 index_ops = self.geom_index_ops
+
             self.geometry_sql.append(
                 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": column,
                     "index_type": self.geom_index_type,
                     "ops": index_ops,
                 }