Fixed #23804 -- Added RasterField for PostGIS.

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

27 files changed, 817 insertions, 238 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
diff --git a/docs/ref/contrib/gis/db-api.txt b/docs/ref/contrib/gis/db-api.txt
index b5d4c4e3ce..cecb02de6c 100644
--- a/docs/ref/contrib/gis/db-api.txt
+++ b/docs/ref/contrib/gis/db-api.txt
@@ -49,8 +49,14 @@ on a different spatial backend.
     lookups and the integrity of your data -- MyISAM tables do
     not support transactions or foreign key constraints.
 
-Creating and Saving Geographic Models
-=====================================
+Raster Support
+--------------
+
+``RasterField`` is currently only implemented for the PostGIS backend. Spatial
+queries (such as lookups and distance) are not yet available for raster fields.
+
+Creating and Saving Models with Geometry Fields
+===============================================
 
 Here is an example of how to create a geometry object (assuming the ``Zipcode``
 model)::
@@ -87,6 +93,42 @@ create a ``GEOSGeometry`` instance from the input.
 For more information creating :class:`~django.contrib.gis.geos.GEOSGeometry`
 objects, refer to the :ref:`GEOS tutorial <geos-tutorial>`.
 
+.. _creating-and-saving-raster-models:
+
+Creating and Saving Models with Raster Fields
+=============================================
+
+.. versionadded:: 1.9
+
+When creating raster models, the raster field will implicitly convert the input
+into a :class:`~django.contrib.gis.gdal.GDALRaster` using lazy-evaluation.
+The raster field will therefore accept any input that is accepted by the
+:class:`~django.contrib.gis.gdal.GDALRaster` constructor.
+
+Here is an example of how to create a raster object from a raster file
+``volcano.tif`` (assuming the ``Elevation`` model)::
+
+    >>> from elevation.models import Elevation
+    >>> dem = Elevation(name='Volcano', rast='/path/to/raster/volcano.tif')
+    >>> dem.save()
+
+:class:`~django.contrib.gis.gdal.GDALRaster` objects may also be used to save
+raster models::
+
+    >>> from django.contrib.gis.gdal import GDALRaster
+    >>> rast = GDALRaster({'width': 10, 'height': 10, 'name': 'Canyon', 'srid': 4326,
+    ...                    'scale': [0.1, -0.1]'bands': [{"data": range(100)}]}
+    >>> dem = Elevation(name='Canyon', rast=rast)
+    >>> dem.save()
+
+Note that this equivalent to::
+
+    >>> dem = Elevation.objects.create(
+    ...     name='Canyon',
+    ...     rast={'width': 10, 'height': 10, 'name': 'Canyon', 'srid': 4326,
+    ...           'scale': [0.1, -0.1]'bands': [{"data": range(100)}]}
+    ... )
+
 .. _spatial-lookups-intro:
 
 Spatial Lookups
@@ -122,6 +164,7 @@ Distance Queries
 
 Introduction
 ------------
+
 Distance calculations with spatial data is tricky because, unfortunately,
 the Earth is not flat.  Some distance queries with fields in a geographic
 coordinate system may have to be expressed differently because of
@@ -132,6 +175,7 @@ in the :doc:`model-api` documentation for more details.
 
 Distance Lookups
 ----------------
+
 *Availability*: PostGIS, Oracle, SpatiaLite
 
 The following distance lookups are available:
diff --git a/docs/ref/contrib/gis/gdal.txt b/docs/ref/contrib/gis/gdal.txt
index 49e26f7a82..25f3ecb153 100644
--- a/docs/ref/contrib/gis/gdal.txt
+++ b/docs/ref/contrib/gis/gdal.txt
@@ -1123,7 +1123,7 @@ blue.
     can be created from different input sources (using the sample data from the
     GeoDjango tests, see also the :ref:`gdal_sample_data` section)::
 
-        >>> from django.contrib.gis.gdal.raster.source import GDALRaster
+        >>> from django.contrib.gis.gdal import GDALRaster
         >>> rst = GDALRaster('/path/to/your/raster.tif', write=False)
         >>> rst.name
         '/path/to/your/raster.tif'
diff --git a/docs/ref/contrib/gis/model-api.txt b/docs/ref/contrib/gis/model-api.txt
index 92f2efc8c8..8f386ef284 100644
--- a/docs/ref/contrib/gis/model-api.txt
+++ b/docs/ref/contrib/gis/model-api.txt
@@ -6,8 +6,8 @@ GeoDjango Model API
    :synopsis: GeoDjango model and field API.
 
 This document explores the details of the GeoDjango Model API.  Throughout this
-section, we'll be using the following geographic model of a `ZIP code`__ as our
-example::
+section, we'll be using the following geographic model of a `ZIP code`__ and
+of a `Digital Elevation Model`__ as our examples::
 
     from django.contrib.gis.db import models
 
@@ -16,13 +16,19 @@ example::
         poly = models.PolygonField()
         objects = models.GeoManager()
 
+    class Elevation(models.Model):
+        name = models.CharField(max_length=100)
+        rast = models.RasterField()
+
 __ http://en.wikipedia.org/wiki/ZIP_code
+__ https://en.wikipedia.org/wiki/Digital_elevation_model
 
-Geometry Field Types
-====================
+Spatial Field Types
+===================
 
-Each of the following geometry field types correspond with the
-OpenGIS Simple Features specification [#fnogc]_.
+Spatial fields consist of a series of geometry field types and one raster field
+type. Each of the geometry field types correspond to the OpenGIS Simple
+Features specification [#fnogc]_. There is no such standard for raster data.
 
 ``GeometryField``
 -----------------
@@ -64,19 +70,31 @@ OpenGIS Simple Features specification [#fnogc]_.
 
 .. class:: GeometryCollectionField
 
-.. _geometry-field-options:
+``RasterField``
+---------------
 
-Geometry Field Options
-======================
+.. versionadded:: 1.9
+
+.. class:: RasterField
+
+``RasterField`` is currently only implemented for the PostGIS backend.
+
+Spatial Field Options
+=====================
+
+.. versionchanged:: 1.9
+
+    The geometry field options ``srid`` and ``spatial_index`` are now shared by
+    ``GeometryField`` and ``RasterField`` through the ``BaseSpatialField``.
 
 In addition to the regular :ref:`common-model-field-options` available for
-Django model fields, geometry fields have the following additional options.
+Django model fields, spatial fields have the following additional options.
 All are optional.
 
 ``srid``
 --------
 
-.. attribute:: GeometryField.srid
+.. attribute:: BaseSpatialField.srid
 
 Sets the SRID [#fnogcsrid]_ (Spatial Reference System Identity) of the geometry field to
 the given value. Defaults to 4326 (also known as `WGS84`__, units are in degrees
@@ -144,7 +162,7 @@ __ http://web.archive.org/web/20080302095452/http://welcome.warnercnr.colostate.
 ``spatial_index``
 -----------------
 
-.. attribute:: GeometryField.spatial_index
+.. attribute:: BaseSpatialField.spatial_index
 
 Defaults to ``True``.  Creates a spatial index for the given geometry
 field.
@@ -157,6 +175,14 @@ field.
     a variant of the R-Tree, while regular database indexes typically
     use B-Trees.
 
+.. _geometry-field-options:
+
+Geometry Field Options
+======================
+
+There are additional options available for Geometry fields. All the following
+options are optional.
+
 ``dim``
 -------
 
@@ -223,7 +249,14 @@ determining `when to use geography data type over geometry data type
 In order to conduct geographic queries, each geographic model requires
 a ``GeoManager`` model manager.  This manager allows for the proper SQL
 construction for geographic queries; thus, without it, all geographic filters
-will fail.  It should also be noted that ``GeoManager`` is required even if the
+will fail.
+
+.. note::
+
+    Geographic filtering support is limited to geometry fields. ``RasterField``
+    does not currently allow spatial querying.
+
+It should also be noted that ``GeoManager`` is required even if the
 model does not have a geographic field itself, e.g., in the case of a
 ``ForeignKey`` relation to a model with a geographic field.  For example,
 if we had an ``Address`` model with a ``ForeignKey`` to our ``Zipcode``
diff --git a/docs/ref/contrib/gis/tutorial.txt b/docs/ref/contrib/gis/tutorial.txt
index 71e5304082..9abcd857ff 100644
--- a/docs/ref/contrib/gis/tutorial.txt
+++ b/docs/ref/contrib/gis/tutorial.txt
@@ -10,10 +10,10 @@ world-class geographic Web framework.  GeoDjango strives to make it as simple
 as possible to create geographic Web applications, like location-based services.
 Its features include:
 
-* Django model fields for `OGC`_ geometries.
+* Django model fields for `OGC`_ geometries and raster data.
 * Extensions to Django's ORM for querying and manipulating spatial data.
-* Loosely-coupled, high-level Python interfaces for GIS geometry operations and
-  data formats.
+* Loosely-coupled, high-level Python interfaces for GIS geometry and raster
+  operations and data manipulation in different formats.
 * Editing geometry fields from the admin.
 
 This tutorial assumes familiarity with Django; thus, if you're brand new to
diff --git a/docs/releases/1.9.txt b/docs/releases/1.9.txt
index 73ed556d31..9641ab2307 100644
--- a/docs/releases/1.9.txt
+++ b/docs/releases/1.9.txt
@@ -168,6 +168,11 @@ Minor features
   Setters for raster properties such as projection or pixel values have
   been added.
 
+* For PostGIS users, the new :class:`~django.contrib.gis.db.models.RasterField`
+  allows :ref:`storing GDALRaster objects <creating-and-saving-raster-models>`.
+  It supports automatic spatial index creation and reprojection when saving a
+  model. It does not yet support spatial querying.
+
 :mod:`django.contrib.messages`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
diff --git a/tests/gis_tests/data/rasters/textrasters.py b/tests/gis_tests/data/rasters/textrasters.py
index b12ce62afc..01b211b3ce 100644
--- a/tests/gis_tests/data/rasters/textrasters.py
+++ b/tests/gis_tests/data/rasters/textrasters.py
@@ -1,7 +1,11 @@
+"""
+Text-based test rasters
+"""
+
 JSON_RASTER = """{
     "srid": 4326,
     "origin": [0, 0],
-    "scale": [1, 1],
+    "scale": [-1, 1],
     "skew": [0, 0],
     "width": 5,
     "height": 5,
diff --git a/tests/gis_tests/gis_migrations/migrations/0001_initial.py b/tests/gis_tests/gis_migrations/migrations/0001_initial.py
index bf2966ee79..6615fdf2c4 100644
--- a/tests/gis_tests/gis_migrations/migrations/0001_initial.py
+++ b/tests/gis_tests/gis_migrations/migrations/0001_initial.py
@@ -1,53 +1,70 @@
-from django.db import migrations, models
+from django.db import connection, migrations, models
 
 from ...models import models as gis_models
 
+ops = [
+    migrations.CreateModel(
+        name='Neighborhood',
+        fields=[
+            ('id', models.AutoField(verbose_name='ID', serialize=False, auto_created=True, primary_key=True)),
+            ('name', models.CharField(max_length=100, unique=True)),
+            ('geom', gis_models.MultiPolygonField(srid=4326)),
+        ],
+        options={
+            'required_db_features': ['gis_enabled'],
+        },
+        bases=(models.Model,),
+    ),
+    migrations.CreateModel(
+        name='Household',
+        fields=[
+            ('id', models.AutoField(verbose_name='ID', serialize=False, auto_created=True, primary_key=True)),
+            ('neighborhood', models.ForeignKey(to='gis_migrations.Neighborhood', to_field='id', null=True)),
+            ('address', models.CharField(max_length=100)),
+            ('zip_code', models.IntegerField(null=True, blank=True)),
+            ('geom', gis_models.PointField(srid=4326, geography=True)),
+        ],
+        options={
+            'required_db_features': ['gis_enabled'],
+        },
+        bases=(models.Model,),
+    ),
+    migrations.CreateModel(
+        name='Family',
+        fields=[
+            ('id', models.AutoField(verbose_name='ID', serialize=False, auto_created=True, primary_key=True)),
+            ('name', models.CharField(max_length=100, unique=True)),
+        ],
+        options={
+        },
+        bases=(models.Model,),
+    ),
+    migrations.AddField(
+        model_name='household',
+        name='family',
+        field=models.ForeignKey(blank=True, to='gis_migrations.Family', null=True),
+        preserve_default=True,
+    )
+]
 
-class Migration(migrations.Migration):
-    """
-    Used for gis.specific migration tests.
-    """
-    operations = [
-        migrations.CreateModel(
-            name='Neighborhood',
-            fields=[
-                ('id', models.AutoField(verbose_name='ID', serialize=False, auto_created=True, primary_key=True)),
-                ('name', models.CharField(max_length=100, unique=True)),
-                ('geom', gis_models.MultiPolygonField(srid=4326)),
-            ],
-            options={
-                'required_db_features': ['gis_enabled'],
-            },
-            bases=(models.Model,),
-        ),
-        migrations.CreateModel(
-            name='Household',
-            fields=[
-                ('id', models.AutoField(verbose_name='ID', serialize=False, auto_created=True, primary_key=True)),
-                ('neighborhood', models.ForeignKey(to='gis_migrations.Neighborhood', to_field='id', null=True)),
-                ('address', models.CharField(max_length=100)),
-                ('zip_code', models.IntegerField(null=True, blank=True)),
-                ('geom', gis_models.PointField(srid=4326, geography=True)),
-            ],
-            options={
-                'required_db_features': ['gis_enabled'],
-            },
-            bases=(models.Model,),
-        ),
+if connection.features.gis_enabled and connection.features.supports_raster:
+    ops += [
         migrations.CreateModel(
-            name='Family',
+            name='Heatmap',
             fields=[
                 ('id', models.AutoField(verbose_name='ID', serialize=False, auto_created=True, primary_key=True)),
                 ('name', models.CharField(max_length=100, unique=True)),
+                ('rast', gis_models.fields.RasterField(srid=4326)),
             ],
             options={
             },
             bases=(models.Model,),
         ),
-        migrations.AddField(
-            model_name='household',
-            name='family',
-            field=models.ForeignKey(blank=True, to='gis_migrations.Family', null=True),
-            preserve_default=True,
-        ),
     ]
+
+
+class Migration(migrations.Migration):
+    """
+    Used for gis-specific migration tests.
+    """
+    operations = ops
diff --git a/tests/gis_tests/gis_migrations/test_commands.py b/tests/gis_tests/gis_migrations/test_commands.py
index b4fe597bca..77bbdf3349 100644
--- a/tests/gis_tests/gis_migrations/test_commands.py
+++ b/tests/gis_tests/gis_migrations/test_commands.py
@@ -39,12 +39,16 @@ class MigrateTests(TransactionTestCase):
         self.assertTableExists("gis_migrations_neighborhood")
         self.assertTableExists("gis_migrations_household")
         self.assertTableExists("gis_migrations_family")
+        if connection.features.supports_raster:
+            self.assertTableExists("gis_migrations_heatmap")
         # Unmigrate everything
         call_command("migrate", "gis_migrations", "zero", verbosity=0)
         # Make sure it's all gone
         self.assertTableNotExists("gis_migrations_neighborhood")
         self.assertTableNotExists("gis_migrations_household")
         self.assertTableNotExists("gis_migrations_family")
+        if connection.features.supports_raster:
+            self.assertTableNotExists("gis_migrations_heatmap")
         # Even geometry columns metadata
         try:
             GeoColumn = connection.ops.geometry_columns()
diff --git a/tests/gis_tests/gis_migrations/test_operations.py b/tests/gis_tests/gis_migrations/test_operations.py
index 7af8da74e5..680742f41d 100644
--- a/tests/gis_tests/gis_migrations/test_operations.py
+++ b/tests/gis_tests/gis_migrations/test_operations.py
@@ -1,14 +1,17 @@
 from __future__ import unicode_literals
 
+from django.contrib.gis.db.models import fields
+from django.core.exceptions import ImproperlyConfigured
 from django.db import connection, migrations, models
 from django.db.migrations.migration import Migration
 from django.db.migrations.state import ProjectState
-from django.test import TransactionTestCase, skipUnlessDBFeature
+from django.test import (
+    TransactionTestCase, skipIfDBFeature, skipUnlessDBFeature,
+)
 
 from ..utils import mysql
 
 if connection.features.gis_enabled:
-    from django.contrib.gis.db.models import fields
     try:
         GeometryColumns = connection.ops.geometry_columns()
         HAS_GEOMETRY_COLUMNS = True
@@ -16,13 +19,14 @@ if connection.features.gis_enabled:
         HAS_GEOMETRY_COLUMNS = False
 
 
-@skipUnlessDBFeature("gis_enabled")
+@skipUnlessDBFeature('gis_enabled')
 class OperationTests(TransactionTestCase):
-    available_apps = ["gis_tests.gis_migrations"]
+    available_apps = ['gis_tests.gis_migrations']
 
     def tearDown(self):
         # Delete table after testing
-        self.apply_operations('gis', self.current_state, [migrations.DeleteModel("Neighborhood")])
+        if hasattr(self, 'current_state'):
+            self.apply_operations('gis', self.current_state, [migrations.DeleteModel('Neighborhood')])
         super(OperationTests, self).tearDown()
 
     def get_table_description(self, table):
@@ -41,19 +45,19 @@ class OperationTests(TransactionTestCase):
         with connection.schema_editor() as editor:
             return migration.apply(project_state, editor)
 
-    def set_up_test_model(self):
-        operations = [migrations.CreateModel(
-            "Neighborhood",
-            [
-                ("id", models.AutoField(primary_key=True)),
-                ('name', models.CharField(max_length=100, unique=True)),
-                ('geom', fields.MultiPolygonField(srid=4326)),
-            ],
-        )]
+    def set_up_test_model(self, force_raster_creation=False):
+        test_fields = [
+            ('id', models.AutoField(primary_key=True)),
+            ('name', models.CharField(max_length=100, unique=True)),
+            ('geom', fields.MultiPolygonField(srid=4326))
+        ]
+        if connection.features.supports_raster or force_raster_creation:
+            test_fields += [('rast', fields.RasterField(srid=4326))]
+        operations = [migrations.CreateModel('Neighborhood', test_fields)]
         return self.apply_operations('gis', ProjectState(), operations)
 
     def assertGeometryColumnsCount(self, expected_count):
-        table_name = "gis_neighborhood"
+        table_name = 'gis_neighborhood'
         if connection.features.uppercases_column_names:
             table_name = table_name.upper()
         self.assertEqual(
@@ -63,91 +67,137 @@ class OperationTests(TransactionTestCase):
             expected_count
         )
 
-    def test_add_gis_field(self):
-        """
-        Tests the AddField operation with a GIS-enabled column.
-        """
+    def assertSpatialIndexExists(self, table, column):
+        with connection.cursor() as cursor:
+            indexes = connection.introspection.get_indexes(cursor, table)
+        self.assertIn(column, indexes)
+
+    def alter_gis_model(self, migration_class, model_name, field_name,
+            blank=False, field_class=None):
         project_state = self.set_up_test_model()
         self.current_state = project_state
-        operation = migrations.AddField(
-            "Neighborhood",
-            "path",
-            fields.LineStringField(srid=4326),
-        )
+        args = [model_name, field_name]
+        if field_class:
+            args.append(field_class(srid=4326, blank=blank))
+        operation = migration_class(*args)
         new_state = project_state.clone()
-        operation.state_forwards("gis", new_state)
+        operation.state_forwards('gis', new_state)
         with connection.schema_editor() as editor:
-            operation.database_forwards("gis", editor, project_state, new_state)
+            operation.database_forwards('gis', editor, project_state, new_state)
         self.current_state = new_state
-        self.assertColumnExists("gis_neighborhood", "path")
+
+    def test_add_geom_field(self):
+        """
+        Test the AddField operation with a geometry-enabled column.
+        """
+        self.alter_gis_model(migrations.AddField, 'Neighborhood',
+            'path', False, fields.LineStringField)
+        self.assertColumnExists('gis_neighborhood', 'path')
 
         # Test GeometryColumns when available
         if HAS_GEOMETRY_COLUMNS:
             self.assertGeometryColumnsCount(2)
 
+        # Test spatial indices when available
         if self.has_spatial_indexes:
-            with connection.cursor() as cursor:
-                indexes = connection.introspection.get_indexes(cursor, "gis_neighborhood")
-            self.assertIn('path', indexes)
+            self.assertSpatialIndexExists('gis_neighborhood', 'path')
 
-    def test_add_blank_gis_field(self):
+    @skipUnlessDBFeature('supports_raster')
+    def test_add_raster_field(self):
+        """
+        Test the AddField operation with a raster-enabled column.
+        """
+        self.alter_gis_model(migrations.AddField, 'Neighborhood',
+            'heatmap', False, fields.RasterField)
+        self.assertColumnExists('gis_neighborhood', 'heatmap')
+
+        # Test spatial indices when available
+        if self.has_spatial_indexes:
+            self.assertSpatialIndexExists('gis_neighborhood', 'heatmap')
+
+    @skipIfDBFeature('supports_raster')
+    def test_create_raster_model_on_db_without_raster_support(self):
+        """
+        Test creating a model with a raster field on a db without raster support.
+        """
+        msg = 'Raster fields require backends with raster support.'
+        with self.assertRaisesMessage(ImproperlyConfigured, msg):
+            self.set_up_test_model(True)
+
+    @skipIfDBFeature('supports_raster')
+    def test_add_raster_field_on_db_without_raster_support(self):
+        """
+        Test adding a raster field on a db without raster support.
+        """
+        msg = 'Raster fields require backends with raster support.'
+        with self.assertRaisesMessage(ImproperlyConfigured, msg):
+            self.alter_gis_model(
+                migrations.AddField, 'Neighborhood', 'heatmap',
+                False, fields.RasterField
+            )
+
+    def test_add_blank_geom_field(self):
         """
         Should be able to add a GeometryField with blank=True.
         """
-        project_state = self.set_up_test_model()
-        self.current_state = project_state
-        operation = migrations.AddField(
-            "Neighborhood",
-            "path",
-            fields.LineStringField(blank=True, srid=4326),
-        )
-        new_state = project_state.clone()
-        operation.state_forwards("gis", new_state)
-        with connection.schema_editor() as editor:
-            operation.database_forwards("gis", editor, project_state, new_state)
-        self.current_state = new_state
-        self.assertColumnExists("gis_neighborhood", "path")
+        self.alter_gis_model(migrations.AddField, 'Neighborhood',
+            'path', True, fields.LineStringField)
+        self.assertColumnExists('gis_neighborhood', 'path')
 
         # Test GeometryColumns when available
         if HAS_GEOMETRY_COLUMNS:
             self.assertGeometryColumnsCount(2)
 
+        # Test spatial indices when available
         if self.has_spatial_indexes:
-            with connection.cursor() as cursor:
-                indexes = connection.introspection.get_indexes(cursor, "gis_neighborhood")
-            self.assertIn('path', indexes)
+            self.assertSpatialIndexExists('gis_neighborhood', 'path')
 
-    def test_remove_gis_field(self):
+    @skipUnlessDBFeature('supports_raster')
+    def test_add_blank_raster_field(self):
         """
-        Tests the RemoveField operation with a GIS-enabled column.
+        Should be able to add a RasterField with blank=True.
         """
-        project_state = self.set_up_test_model()
-        self.current_state = project_state
-        operation = migrations.RemoveField("Neighborhood", "geom")
-        new_state = project_state.clone()
-        operation.state_forwards("gis", new_state)
-        with connection.schema_editor() as editor:
-            operation.database_forwards("gis", editor, project_state, new_state)
-        self.current_state = new_state
-        self.assertColumnNotExists("gis_neighborhood", "geom")
+        self.alter_gis_model(migrations.AddField, 'Neighborhood',
+            'heatmap', True, fields.RasterField)
+        self.assertColumnExists('gis_neighborhood', 'heatmap')
+
+        # Test spatial indices when available
+        if self.has_spatial_indexes:
+            self.assertSpatialIndexExists('gis_neighborhood', 'heatmap')
+
+    def test_remove_geom_field(self):
+        """
+        Test the RemoveField operation with a geometry-enabled column.
+        """
+        self.alter_gis_model(migrations.RemoveField, 'Neighborhood', 'geom')
+        self.assertColumnNotExists('gis_neighborhood', 'geom')
 
         # Test GeometryColumns when available
         if HAS_GEOMETRY_COLUMNS:
             self.assertGeometryColumnsCount(0)
 
+    @skipUnlessDBFeature('supports_raster')
+    def test_remove_raster_field(self):
+        """
+        Test the RemoveField operation with a raster-enabled column.
+        """
+        self.alter_gis_model(migrations.RemoveField, 'Neighborhood', 'rast')
+        self.assertColumnNotExists('gis_neighborhood', 'rast')
+
     def test_create_model_spatial_index(self):
         self.current_state = self.set_up_test_model()
 
         if not self.has_spatial_indexes:
-            self.skipTest("No support for Spatial indexes")
+            self.skipTest('No support for Spatial indexes')
 
-        with connection.cursor() as cursor:
-            indexes = connection.introspection.get_indexes(cursor, "gis_neighborhood")
-        self.assertIn('geom', indexes)
+        self.assertSpatialIndexExists('gis_neighborhood', 'geom')
+
+        if connection.features.supports_raster:
+            self.assertSpatialIndexExists('gis_neighborhood', 'rast')
 
     @property
     def has_spatial_indexes(self):
         if mysql:
             with connection.cursor() as cursor:
-                return connection.introspection.supports_spatial_index(cursor, "gis_neighborhood")
+                return connection.introspection.supports_spatial_index(cursor, 'gis_neighborhood')
         return True
diff --git a/tests/gis_tests/models.py b/tests/gis_tests/models.py
index cfb7f96656..f5c694e688 100644
--- a/tests/gis_tests/models.py
+++ b/tests/gis_tests/models.py
@@ -1,14 +1,18 @@
 from django.core.exceptions import ImproperlyConfigured
+from django.db import models
+
+
+class DummyField(models.Field):
+    def __init__(self, dim=None, srid=None, geography=None, spatial_index=True, *args, **kwargs):
+        super(DummyField, self).__init__(*args, **kwargs)
 
 try:
     from django.contrib.gis.db import models
+    try:
+        models.RasterField()
+    except ImproperlyConfigured:
+        models.RasterField = DummyField
 except ImproperlyConfigured:
-    from django.db import models
-
-    class DummyField(models.Field):
-        def __init__(self, dim=None, srid=None, geography=None, *args, **kwargs):
-            super(DummyField, self).__init__(*args, **kwargs)
-
     models.GeoManager = models.Manager
     models.GeometryField = DummyField
     models.LineStringField = DummyField
@@ -16,3 +20,4 @@ except ImproperlyConfigured:
     models.MultiPolygonField = DummyField
     models.PointField = DummyField
     models.PolygonField = DummyField
+    models.RasterField = DummyField
diff --git a/tests/gis_tests/rasterapp/__init__.py b/tests/gis_tests/rasterapp/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
--- /dev/null
+++ b/tests/gis_tests/rasterapp/__init__.py
diff --git a/tests/gis_tests/rasterapp/models.py b/tests/gis_tests/rasterapp/models.py
new file mode 100644
index 0000000000..d7cc4d15e2
--- /dev/null
+++ b/tests/gis_tests/rasterapp/models.py
@@ -0,0 +1,11 @@
+from ..models import models
+
+
+class RasterModel(models.Model):
+    rast = models.RasterField(null=True, srid=4326, spatial_index=True, blank=True)
+
+    class Meta:
+        required_db_features = ['supports_raster']
+
+    def __str__(self):
+        return str(self.id)
diff --git a/tests/gis_tests/rasterapp/test_rasterfield.py b/tests/gis_tests/rasterapp/test_rasterfield.py
new file mode 100644
index 0000000000..7b4413ae0d
--- /dev/null
+++ b/tests/gis_tests/rasterapp/test_rasterfield.py
@@ -0,0 +1,72 @@
+import json
+
+from django.contrib.gis.shortcuts import numpy
+from django.test import TransactionTestCase, skipUnlessDBFeature
+
+from ..data.rasters.textrasters import JSON_RASTER
+from .models import RasterModel
+
+
+@skipUnlessDBFeature('supports_raster')
+class RasterFieldTest(TransactionTestCase):
+    available_apps = ['gis_tests.rasterapp']
+
+    def test_field_null_value(self):
+        """
+        Test creating a model where the RasterField has a null value.
+        """
+        r = RasterModel.objects.create(rast=None)
+        r.refresh_from_db()
+        self.assertIsNone(r.rast)
+
+    def test_model_creation(self):
+        """
+        Test RasterField through a test model.
+        """
+        # Create model instance from JSON raster
+        r = RasterModel.objects.create(rast=JSON_RASTER)
+        r.refresh_from_db()
+        # Test raster metadata properties
+        self.assertEqual((5, 5), (r.rast.width, r.rast.height))
+        self.assertEqual([0.0, -1.0, 0.0, 0.0, 0.0, 1.0], r.rast.geotransform)
+        self.assertIsNone(r.rast.bands[0].nodata_value)
+        # Compare srs
+        self.assertEqual(r.rast.srs.srid, 4326)
+        # Compare pixel values
+        band = r.rast.bands[0].data()
+        # If numpy, convert result to list
+        if numpy:
+            band = band.flatten().tolist()
+        # Loop through rows in band data and assert single
+        # value is as expected.
+        self.assertEqual(
+            [
+                0.0, 1.0, 2.0, 3.0, 4.0,
+                5.0, 6.0, 7.0, 8.0, 9.0,
+                10.0, 11.0, 12.0, 13.0, 14.0,
+                15.0, 16.0, 17.0, 18.0, 19.0,
+                20.0, 21.0, 22.0, 23.0, 24.0
+            ],
+            band
+        )
+
+    def test_implicit_raster_transformation(self):
+        """
+        Test automatic transformation of rasters with srid different from the
+        field srid.
+        """
+        # Parse json raster
+        rast = json.loads(JSON_RASTER)
+        # Update srid to another value
+        rast['srid'] = 3086
+        # Save model and get it from db
+        r = RasterModel.objects.create(rast=rast)
+        r.refresh_from_db()
+        # Confirm raster has been transformed to the default srid
+        self.assertEqual(r.rast.srs.srid, 4326)
+        # Confirm geotransform is in lat/lon
+        self.assertEqual(
+            r.rast.geotransform,
+            [-87.9298551266551, 9.459646421449934e-06, 0.0,
+             23.94249275457565, 0.0, -9.459646421449934e-06]
+        )