Fixed #23804 -- Added RasterField for PostGIS.

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

14 files changed, 444 insertions, 110 deletions

diff --git a/django/contrib/gis/db/backends/base/features.py b/django/contrib/gis/db/backends/base/features.py
index ebd26af431..4fd4fd18a7 100644
--- a/django/contrib/gis/db/backends/base/features.py
+++ b/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
diff --git a/django/contrib/gis/db/backends/postgis/const.py b/django/contrib/gis/db/backends/postgis/const.py
new file mode 100644
index 0000000000..0d6d809ec3
--- /dev/null
+++ b/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
+}
diff --git a/django/contrib/gis/db/backends/postgis/features.py b/django/contrib/gis/db/backends/postgis/features.py
index 181c859641..3fad9e70c0 100644
--- a/django/contrib/gis/db/backends/postgis/features.py
+++ b/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
diff --git a/django/contrib/gis/db/backends/postgis/introspection.py b/django/contrib/gis/db/backends/postgis/introspection.py
index 5e08b2f9ce..7f231a578a 100644
--- a/django/contrib/gis/db/backends/postgis/introspection.py
+++ b/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
diff --git a/django/contrib/gis/db/backends/postgis/operations.py b/django/contrib/gis/db/backends/postgis/operations.py
index 3e7aa5cdde..31406e6506 100644
--- a/django/contrib/gis/db/backends/postgis/operations.py
+++ b/django/contrib/gis/db/backends/postgis/operations.py
@@ -4,6 +4,9 @@ from django.conf import settings
 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.db.utils import ProgrammingError
 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 @@ class PostGISOperations(BaseSpatialOperations, DatabaseOperations):
 
     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 @@ class PostGISOperations(BaseSpatialOperations, DatabaseOperations):
         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 @@ class PostGISOperations(BaseSpatialOperations, DatabaseOperations):
 
     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)
diff --git a/django/contrib/gis/db/backends/postgis/pgraster.py b/django/contrib/gis/db/backends/postgis/pgraster.py
new file mode 100644
index 0000000000..555fc17f1e
--- /dev/null
+++ b/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 it 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()
diff --git a/django/contrib/gis/db/backends/postgis/schema.py b/django/contrib/gis/db/backends/postgis/schema.py
index be2b12bb66..8c2cb38608 100644
--- a/django/contrib/gis/db/backends/postgis/schema.py
+++ b/django/contrib/gis/db/backends/postgis/schema.py
@@ -4,6 +4,7 @@ from django.db.backends.postgresql_psycopg2.schema import DatabaseSchemaEditor
 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,8 +17,8 @@ class PostGISSchemaEditor(DatabaseSchemaEditor):
         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)
@@ -25,8 +26,13 @@ class PostGISSchemaEditor(DatabaseSchemaEditor):
         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 @@ class PostGISSchemaEditor(DatabaseSchemaEditor):
                 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,
                 }
diff --git a/django/contrib/gis/db/models/__init__.py b/django/contrib/gis/db/models/__init__.py
index 835e907526..bbba20dca8 100644
--- a/django/contrib/gis/db/models/__init__.py
+++ b/django/contrib/gis/db/models/__init__.py
@@ -11,4 +11,4 @@ from django.contrib.gis.db.models.manager import GeoManager  # NOQA
 from django.contrib.gis.db.models.fields import (  # NOQA
     GeometryField, PointField, LineStringField, PolygonField,
     MultiPointField, MultiLineStringField, MultiPolygonField,
-    GeometryCollectionField)
+    GeometryCollectionField, RasterField)
diff --git a/django/contrib/gis/db/models/fields.py b/django/contrib/gis/db/models/fields.py
index 6d5c2b2b6a..5cdb81e440 100644
--- a/django/contrib/gis/db/models/fields.py
+++ b/django/contrib/gis/db/models/fields.py
@@ -1,7 +1,9 @@
 from django.contrib.gis import forms
 from django.contrib.gis.db.models.lookups import gis_lookups
-from django.contrib.gis.db.models.proxy import GeometryProxy
+from django.contrib.gis.db.models.proxy import SpatialProxy
+from django.contrib.gis.gdal.raster.source import GDALRaster
 from django.contrib.gis.geometry.backend import Geometry, GeometryException
+from django.core.exceptions import ImproperlyConfigured
 from django.db.models.expressions import Expression
 from django.db.models.fields import Field
 from django.utils import six
@@ -65,22 +67,21 @@ class GeoSelectFormatMixin(object):
         return sel_fmt % sql, params
 
 
-class GeometryField(GeoSelectFormatMixin, Field):
-    "The base GIS field -- maps to the OpenGIS Specification Geometry type."
-
-    # The OpenGIS Geometry name.
-    geom_type = 'GEOMETRY'
-    form_class = forms.GeometryField
+class BaseSpatialField(Field):
+    """
+    The Base GIS Field.
 
+    It's used as a base class for GeometryField and RasterField. Defines
+    properties that are common to all GIS fields such as the characteristics
+    of the spatial reference system of the field.
+    """
+    description = _("The base GIS field.")
     # Geodetic units.
     geodetic_units = ('decimal degree', 'degree')
 
-    description = _("The base GIS field -- maps to the OpenGIS Specification Geometry type.")
-
-    def __init__(self, verbose_name=None, srid=4326, spatial_index=True, dim=2,
-                 geography=False, **kwargs):
+    def __init__(self, verbose_name=None, srid=4326, spatial_index=True, **kwargs):
         """
-        The initialization function for geometry fields.  Takes the following
+        The initialization function for base spatial fields. Takes the following
         as keyword arguments:
 
         srid:
@@ -91,18 +92,6 @@ class GeometryField(GeoSelectFormatMixin, Field):
          Indicates whether to create a spatial index.  Defaults to True.
          Set this instead of 'db_index' for geographic fields since index
          creation is different for geometry columns.
-
-        dim:
-         The number of dimensions for this geometry.  Defaults to 2.
-
-        extent:
-         Customize the extent, in a 4-tuple of WGS 84 coordinates, for the
-         geometry field entry in the `USER_SDO_GEOM_METADATA` table.  Defaults
-         to (-180.0, -90.0, 180.0, 90.0).
-
-        tolerance:
-         Define the tolerance, in meters, to use for the geometry field
-         entry in the `USER_SDO_GEOM_METADATA` table.  Defaults to 0.05.
         """
 
         # Setting the index flag with the value of the `spatial_index` keyword.
@@ -112,38 +101,26 @@ class GeometryField(GeoSelectFormatMixin, Field):
         # easily available in the field instance for distance queries.
         self.srid = srid
 
-        # Setting the dimension of the geometry field.
-        self.dim = dim
-
         # Setting the verbose_name keyword argument with the positional
         # first parameter, so this works like normal fields.
         kwargs['verbose_name'] = verbose_name
 
-        # Is this a geography rather than a geometry column?
-        self.geography = geography
-
-        # Oracle-specific private attributes for creating the entry in
-        # `USER_SDO_GEOM_METADATA`
-        self._extent = kwargs.pop('extent', (-180.0, -90.0, 180.0, 90.0))
-        self._tolerance = kwargs.pop('tolerance', 0.05)
-
-        super(GeometryField, self).__init__(**kwargs)
+        super(BaseSpatialField, self).__init__(**kwargs)
 
     def deconstruct(self):
-        name, path, args, kwargs = super(GeometryField, self).deconstruct()
-        # Always include SRID for less fragility; include others if they're
-        # not the default values.
+        name, path, args, kwargs = super(BaseSpatialField, self).deconstruct()
+        # Always include SRID for less fragility; include spatial index if it's
+        # not the default value.
         kwargs['srid'] = self.srid
-        if self.dim != 2:
-            kwargs['dim'] = self.dim
         if self.spatial_index is not True:
             kwargs['spatial_index'] = self.spatial_index
-        if self.geography is not False:
-            kwargs['geography'] = self.geography
         return name, path, args, kwargs
 
+    def db_type(self, connection):
+        return connection.ops.geo_db_type(self)
+
     # The following functions are used to get the units, their name, and
-    # the spheroid corresponding to the SRID of the GeometryField.
+    # the spheroid corresponding to the SRID of the BaseSpatialField.
     def _get_srid_info(self, connection):
         # Get attributes from `get_srid_info`.
         self._units, self._units_name, self._spheroid = get_srid_info(self.srid, connection)
@@ -163,7 +140,6 @@ class GeometryField(GeoSelectFormatMixin, Field):
             self._get_srid_info(connection)
         return self._units_name
 
-    # ### Routines specific to GeometryField ###
     def geodetic(self, connection):
         """
         Returns true if this field's SRID corresponds with a coordinate
@@ -174,6 +150,64 @@ class GeometryField(GeoSelectFormatMixin, Field):
         # test if srid is 4326 (WGS84), even if this is over-simplification.
         return units_name.lower() in self.geodetic_units if units_name else self.srid == 4326
 
+    def get_placeholder(self, value, compiler, connection):
+        """
+        Returns the placeholder for the spatial column for the
+        given value.
+        """
+        return connection.ops.get_geom_placeholder(self, value, compiler)
+
+
+class GeometryField(GeoSelectFormatMixin, BaseSpatialField):
+    """
+    The base Geometry field -- maps to the OpenGIS Specification Geometry type.
+    """
+    description = _("The base Geometry field -- maps to the OpenGIS Specification Geometry type.")
+    form_class = forms.GeometryField
+    # The OpenGIS Geometry name.
+    geom_type = 'GEOMETRY'
+
+    def __init__(self, verbose_name=None, dim=2, geography=False, **kwargs):
+        """
+        The initialization function for geometry fields. In addition to the
+        parameters from BaseSpatialField, it takes the following as keyword
+        arguments:
+
+        dim:
+         The number of dimensions for this geometry.  Defaults to 2.
+
+        extent:
+         Customize the extent, in a 4-tuple of WGS 84 coordinates, for the
+         geometry field entry in the `USER_SDO_GEOM_METADATA` table.  Defaults
+         to (-180.0, -90.0, 180.0, 90.0).
+
+        tolerance:
+         Define the tolerance, in meters, to use for the geometry field
+         entry in the `USER_SDO_GEOM_METADATA` table.  Defaults to 0.05.
+        """
+        # Setting the dimension of the geometry field.
+        self.dim = dim
+
+        # Is this a geography rather than a geometry column?
+        self.geography = geography
+
+        # Oracle-specific private attributes for creating the entry in
+        # `USER_SDO_GEOM_METADATA`
+        self._extent = kwargs.pop('extent', (-180.0, -90.0, 180.0, 90.0))
+        self._tolerance = kwargs.pop('tolerance', 0.05)
+
+        super(GeometryField, self).__init__(verbose_name=verbose_name, **kwargs)
+
+    def deconstruct(self):
+        name, path, args, kwargs = super(GeometryField, self).deconstruct()
+        # Include kwargs if they're not the default values.
+        if self.dim != 2:
+            kwargs['dim'] = self.dim
+        if self.geography is not False:
+            kwargs['geography'] = self.geography
+        return name, path, args, kwargs
+
+    # ### Routines specific to GeometryField ###
     def get_distance(self, value, lookup_type, connection):
         """
         Returns a distance number in units of the field.  For example, if
@@ -244,10 +278,7 @@ class GeometryField(GeoSelectFormatMixin, Field):
         super(GeometryField, self).contribute_to_class(cls, name, **kwargs)
 
         # Setup for lazy-instantiated Geometry object.
-        setattr(cls, self.attname, GeometryProxy(Geometry, self))
-
-    def db_type(self, connection):
-        return connection.ops.geo_db_type(self)
+        setattr(cls, self.attname, SpatialProxy(Geometry, self))
 
     def formfield(self, **kwargs):
         defaults = {'form_class': self.form_class,
@@ -309,13 +340,6 @@ class GeometryField(GeoSelectFormatMixin, Field):
         else:
             return connection.ops.Adapter(self.get_prep_value(value))
 
-    def get_placeholder(self, value, compiler, connection):
-        """
-        Returns the placeholder for the geometry column for the
-        given value.
-        """
-        return connection.ops.get_geom_placeholder(self, value, compiler)
-
 
 for klass in gis_lookups.values():
     GeometryField.register_lookup(klass)
@@ -371,3 +395,39 @@ class ExtentField(GeoSelectFormatMixin, Field):
 
     def get_internal_type(self):
         return "ExtentField"
+
+
+class RasterField(BaseSpatialField):
+    """
+    Raster field for GeoDjango -- evaluates into GDALRaster objects.
+    """
+
+    description = _("Raster Field")
+    geom_type = 'RASTER'
+
+    def _check_connection(self, connection):
+        # Make sure raster fields are used only on backends with raster support.
+        if not connection.features.gis_enabled or not connection.features.supports_raster:
+            raise ImproperlyConfigured('Raster fields require backends with raster support.')
+
+    def db_type(self, connection):
+        self._check_connection(connection)
+        return super(RasterField, self).db_type(connection)
+
+    def from_db_value(self, value, expression, connection, context):
+        return connection.ops.parse_raster(value)
+
+    def get_db_prep_value(self, value, connection, prepared=False):
+        self._check_connection(connection)
+        # Prepare raster for writing to database.
+        if not prepared:
+            value = connection.ops.deconstruct_raster(value)
+        return super(RasterField, self).get_db_prep_value(value, connection, prepared)
+
+    def contribute_to_class(self, cls, name, **kwargs):
+        super(RasterField, self).contribute_to_class(cls, name, **kwargs)
+        # Setup for lazy-instantiated Raster object. For large querysets, the
+        # instantiation of all GDALRasters can potentially be expensive. This
+        # delays the instantiation of the objects to the moment of evaluation
+        # of the raster attribute.
+        setattr(cls, self.attname, SpatialProxy(GDALRaster, self))
diff --git a/django/contrib/gis/db/models/proxy.py b/django/contrib/gis/db/models/proxy.py
index 8f2b04c452..494d4a949d 100644
--- a/django/contrib/gis/db/models/proxy.py
+++ b/django/contrib/gis/db/models/proxy.py
@@ -1,66 +1,73 @@
 """
-The GeometryProxy object, allows for lazy-geometries.  The proxy uses
-Python descriptors for instantiating and setting Geometry objects
-corresponding to geographic model fields.
+The SpatialProxy object allows for lazy-geometries and lazy-rasters. The proxy
+uses Python descriptors for instantiating and setting Geometry or Raster
+objects corresponding to geographic model fields.
 
 Thanks to Robert Coup for providing this functionality (see #4322).
 """
 from django.utils import six
 
 
-class GeometryProxy(object):
+class SpatialProxy(object):
     def __init__(self, klass, field):
         """
-        Proxy initializes on the given Geometry class (not an instance) and
-        the GeometryField.
+        Proxy initializes on the given Geometry or Raster class (not an instance)
+        and the corresponding field.
         """
         self._field = field
         self._klass = klass
 
     def __get__(self, obj, type=None):
         """
-        This accessor retrieves the geometry, initializing it using the geometry
-        class specified during initialization and the HEXEWKB value of the field.
-        Currently, only GEOS or OGR geometries are supported.
+        This accessor retrieves the geometry or raster, initializing it using
+        the corresponding class specified during initialization and the value
+        of the field. Currently, GEOS or OGR geometries as well as GDALRasters
+        are supported.
         """
         if obj is None:
             # Accessed on a class, not an instance
             return self
 
         # Getting the value of the field.
-        geom_value = obj.__dict__[self._field.attname]
+        geo_value = obj.__dict__[self._field.attname]
 
-        if isinstance(geom_value, self._klass):
-            geom = geom_value
-        elif (geom_value is None) or (geom_value == ''):
-            geom = None
+        if isinstance(geo_value, self._klass):
+            geo_obj = geo_value
+        elif (geo_value is None) or (geo_value == ''):
+            geo_obj = None
         else:
-            # Otherwise, a Geometry object is built using the field's contents,
-            # and the model's corresponding attribute is set.
-            geom = self._klass(geom_value)
-            setattr(obj, self._field.attname, geom)
-        return geom
+            # Otherwise, a geometry or raster object is built using the field's
+            # contents, and the model's corresponding attribute is set.
+            geo_obj = self._klass(geo_value)
+            setattr(obj, self._field.attname, geo_obj)
+        return geo_obj
 
     def __set__(self, obj, value):
         """
-        This accessor sets the proxied geometry with the geometry class
-        specified during initialization.  Values of None, HEXEWKB, or WKT may
-        be used to set the geometry as well.
+        This accessor sets the proxied geometry or raster with the
+        corresponding class specified during initialization.
+
+        To set geometries, values of None, HEXEWKB, or WKT may be used.
+        To set rasters, JSON or dict values may be used.
         """
-        # The OGC Geometry type of the field.
+        # The geographic type of the field.
         gtype = self._field.geom_type
 
-        # The geometry type must match that of the field -- unless the
-        # general GeometryField is used.
-        if isinstance(value, self._klass) and (str(value.geom_type).upper() == gtype or gtype == 'GEOMETRY'):
-            # Assigning the SRID to the geometry.
+        if gtype == 'RASTER' and (value is None or isinstance(value, six.string_types + (dict, self._klass))):
+            # For raster fields, assure input is None or a string, dict, or
+            # raster instance.
+            pass
+        elif isinstance(value, self._klass) and (str(value.geom_type).upper() == gtype or gtype == 'GEOMETRY'):
+            # The geometry type must match that of the field -- unless the
+            # general GeometryField is used.
             if value.srid is None:
+                # Assigning the field SRID if the geometry has no SRID.
                 value.srid = self._field.srid
         elif value is None or isinstance(value, six.string_types + (six.memoryview,)):
-            # Set with None, WKT, HEX, or WKB
+            # Set geometries with None, WKT, HEX, or WKB
             pass
         else:
-            raise TypeError('Cannot set %s GeometryProxy (%s) with value of type: %s' % (
+            raise TypeError('Cannot set %s SpatialProxy (%s) with value of type: %s' % (
                 obj.__class__.__name__, gtype, type(value)))
 
         # Setting the objects dictionary with the value, and returning.
diff --git a/django/contrib/gis/gdal/raster/band.py b/django/contrib/gis/gdal/raster/band.py
index b5aa5aa4df..f1eb50e506 100644
--- a/django/contrib/gis/gdal/raster/band.py
+++ b/django/contrib/gis/gdal/raster/band.py
@@ -6,7 +6,7 @@ from django.contrib.gis.shortcuts import numpy
 from django.utils import six
 from django.utils.encoding import force_text
 
-from .const import GDAL_PIXEL_TYPES, GDAL_TO_CTYPES
+from .const import GDAL_INTEGER_TYPES, GDAL_PIXEL_TYPES, GDAL_TO_CTYPES
 
 
 class GDALBand(GDALBase):
@@ -64,9 +64,15 @@ class GDALBand(GDALBase):
         """
         Returns the nodata value for this band, or None if it isn't set.
         """
+        # Get value and nodata exists flag
         nodata_exists = c_int()
         value = capi.get_band_nodata_value(self._ptr, nodata_exists)
-        return value if nodata_exists else None
+        if not nodata_exists:
+            value = None
+        # If the pixeltype is an integer, convert to int
+        elif self.datatype() in GDAL_INTEGER_TYPES:
+            value = int(value)
+        return value
 
     @nodata_value.setter
     def nodata_value(self, value):
diff --git a/django/contrib/gis/gdal/raster/const.py b/django/contrib/gis/gdal/raster/const.py
index bdf858afaf..719d1bf0b2 100644
--- a/django/contrib/gis/gdal/raster/const.py
+++ b/django/contrib/gis/gdal/raster/const.py
@@ -21,6 +21,9 @@ GDAL_PIXEL_TYPES = {
     11: 'GDT_CFloat64',  # Complex Float64
 }
 
+# A list of gdal datatypes that are integers.
+GDAL_INTEGER_TYPES = [1, 2, 3, 4, 5]
+
 # Lookup values to convert GDAL pixel type indices into ctypes objects.
 # The GDAL band-io works with ctypes arrays to hold data to be written
 # or to hold the space for data to be read into. The lookup below helps
diff --git a/django/contrib/gis/gdal/raster/source.py b/django/contrib/gis/gdal/raster/source.py
index 6f36a10823..34f941f40d 100644
--- a/django/contrib/gis/gdal/raster/source.py
+++ b/django/contrib/gis/gdal/raster/source.py
@@ -111,7 +111,7 @@ class GDALRaster(GDALBase):
                 if 'nodata_value' in band_input:
                     self.bands[i].nodata_value = band_input['nodata_value']
 
-            # Set SRID, default to 0 (this assures SRS is always instanciated)
+            # Set SRID
             self.srs = ds_input.get('srid')
 
             # Set additional properties if provided
diff --git a/django/db/backends/postgresql_psycopg2/introspection.py b/django/db/backends/postgresql_psycopg2/introspection.py
index 32c69bc510..9b3e9074b2 100644
--- a/django/db/backends/postgresql_psycopg2/introspection.py
+++ b/django/db/backends/postgresql_psycopg2/introspection.py
@@ -34,6 +34,16 @@ class DatabaseIntrospection(BaseDatabaseIntrospection):
 
     ignored_tables = []
 
+    _get_indexes_query = """
+        SELECT 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
+        WHERE c.oid = idx.indrelid
+            AND idx.indexrelid = c2.oid
+            AND attr.attrelid = c.oid
+            AND attr.attnum = idx.indkey[0]
+            AND c.relname = %s"""
+
     def get_field_type(self, data_type, description):
         field_type = super(DatabaseIntrospection, self).get_field_type(data_type, description)
         if field_type == 'IntegerField' and description.default and 'nextval' in description.default:
@@ -108,15 +118,7 @@ class DatabaseIntrospection(BaseDatabaseIntrospection):
     def get_indexes(self, cursor, table_name):
         # This query retrieves each index on the given table, including the
         # first associated field name
-        cursor.execute("""
-            SELECT 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
-            WHERE c.oid = idx.indrelid
-                AND idx.indexrelid = c2.oid
-                AND attr.attrelid = c.oid
-                AND attr.attnum = idx.indkey[0]
-                AND c.relname = %s""", [table_name])
+        cursor.execute(self._get_indexes_query, [table_name])
         indexes = {}
         for row in cursor.fetchall():
             # row[1] (idx.indkey) is stored in the DB as an array. It comes out as