Fixed #25588 -- Added spatial lookups to RasterField.

Thanks Tim Graham for the review.

11 files changed, 744 insertions, 159 deletions

diff --git a/django/contrib/gis/db/backends/postgis/adapter.py b/django/contrib/gis/db/backends/postgis/adapter.py
index cb0d466a80..2d2d9578ec 100644
--- a/django/contrib/gis/db/backends/postgis/adapter.py
+++ b/django/contrib/gis/db/backends/postgis/adapter.py
@@ -6,16 +6,27 @@ from __future__ import unicode_literals
 from psycopg2 import Binary
 from psycopg2.extensions import ISQLQuote
 
+from django.contrib.gis.db.backends.postgis.pgraster import to_pgraster
+from django.contrib.gis.geometry.backend import Geometry
+
 
 class PostGISAdapter(object):
-    def __init__(self, geom, geography=False):
-        "Initializes on the geometry."
+    def __init__(self, obj, geography=False):
+        """
+        Initialize on the spatial object.
+        """
+        self.is_geometry = isinstance(obj, Geometry)
+
         # Getting the WKB (in string form, to allow easy pickling of
-        # the adaptor) and the SRID from the geometry.
-        self.ewkb = bytes(geom.ewkb)
-        self.srid = geom.srid
+        # the adaptor) and the SRID from the geometry or raster.
+        if self.is_geometry:
+            self.ewkb = bytes(obj.ewkb)
+            self._adapter = Binary(self.ewkb)
+        else:
+            self.ewkb = to_pgraster(obj)
+
+        self.srid = obj.srid
         self.geography = geography
-        self._adapter = Binary(self.ewkb)
 
     def __conform__(self, proto):
         # Does the given protocol conform to what Psycopg2 expects?
@@ -40,12 +51,19 @@ class PostGISAdapter(object):
         This method allows escaping the binary in the style required by the
         server's `standard_conforming_string` setting.
         """
-        self._adapter.prepare(conn)
+        if self.is_geometry:
+            self._adapter.prepare(conn)
 
     def getquoted(self):
-        "Returns a properly quoted string for use in PostgreSQL/PostGIS."
-        # psycopg will figure out whether to use E'\\000' or '\000'
-        return str('%s(%s)' % (
-            'ST_GeogFromWKB' if self.geography else 'ST_GeomFromEWKB',
-            self._adapter.getquoted().decode())
-        )
+        """
+        Return a properly quoted string for use in PostgreSQL/PostGIS.
+        """
+        if self.is_geometry:
+            # Psycopg will figure out whether to use E'\\000' or '\000'.
+            return str('%s(%s)' % (
+                'ST_GeogFromWKB' if self.geography else 'ST_GeomFromEWKB',
+                self._adapter.getquoted().decode())
+            )
+        else:
+            # For rasters, add explicit type cast to WKB string.
+            return "'%s'::raster" % self.ewkb
diff --git a/django/contrib/gis/db/backends/postgis/operations.py b/django/contrib/gis/db/backends/postgis/operations.py
index 89ddc11127..8216668ca5 100644
--- a/django/contrib/gis/db/backends/postgis/operations.py
+++ b/django/contrib/gis/db/backends/postgis/operations.py
@@ -4,30 +4,83 @@ from django.conf import settings
 from django.contrib.gis.db.backends.base.operations import \
     BaseSpatialOperations
 from django.contrib.gis.db.backends.utils import SpatialOperator
+from django.contrib.gis.gdal import GDALRaster
 from django.contrib.gis.geometry.backend import Geometry
 from django.contrib.gis.measure import Distance
 from django.core.exceptions import ImproperlyConfigured
 from django.db.backends.postgresql.operations import DatabaseOperations
 from django.db.utils import ProgrammingError
+from django.utils import six
 from django.utils.functional import cached_property
 
 from .adapter import PostGISAdapter
 from .models import PostGISGeometryColumns, PostGISSpatialRefSys
 from .pgraster import from_pgraster, get_pgraster_srid, to_pgraster
 
+# Identifier to mark raster lookups as bilateral.
+BILATERAL = 'bilateral'
+
 
 class PostGISOperator(SpatialOperator):
-    def __init__(self, geography=False, **kwargs):
-        # Only a subset of the operators and functions are available
-        # for the geography type.
+    def __init__(self, geography=False, raster=False, **kwargs):
+        # Only a subset of the operators and functions are available for the
+        # geography type.
         self.geography = geography
+        # Only a subset of the operators and functions are available for the
+        # raster type. Lookups that don't suport raster will be converted to
+        # polygons. If the raster argument is set to BILATERAL, then the
+        # operator cannot handle mixed geom-raster lookups.
+        self.raster = raster
         super(PostGISOperator, self).__init__(**kwargs)
 
-    def as_sql(self, connection, lookup, *args):
+    def as_sql(self, connection, lookup, template_params, *args):
         if lookup.lhs.output_field.geography and not self.geography:
             raise ValueError('PostGIS geography does not support the "%s" '
                              'function/operator.' % (self.func or self.op,))
-        return super(PostGISOperator, self).as_sql(connection, lookup, *args)
+
+        template_params = self.check_raster(lookup, template_params)
+        return super(PostGISOperator, self).as_sql(connection, lookup, template_params, *args)
+
+    def check_raster(self, lookup, template_params):
+        # Get rhs value.
+        if isinstance(lookup.rhs, (tuple, list)):
+            rhs_val = lookup.rhs[0]
+            spheroid = lookup.rhs[-1] == 'spheroid'
+        else:
+            rhs_val = lookup.rhs
+            spheroid = False
+
+        # Check which input is a raster.
+        lhs_is_raster = lookup.lhs.field.geom_type == 'RASTER'
+        rhs_is_raster = isinstance(rhs_val, GDALRaster)
+
+        # Look for band indices and inject them if provided.
+        if lookup.band_lhs is not None and lhs_is_raster:
+            if not self.func:
+                raise ValueError('Band indices are not allowed for this operator, it works on bbox only.')
+            template_params['lhs'] = '%s, %s' % (template_params['lhs'], lookup.band_lhs)
+
+        if lookup.band_rhs is not None and rhs_is_raster:
+            if not self.func:
+                raise ValueError('Band indices are not allowed for this operator, it works on bbox only.')
+            template_params['rhs'] = '%s, %s' % (template_params['rhs'], lookup.band_rhs)
+
+        # Convert rasters to polygons if necessary.
+        if not self.raster or spheroid:
+            # Operators without raster support.
+            if lhs_is_raster:
+                template_params['lhs'] = 'ST_Polygon(%s)' % template_params['lhs']
+            if rhs_is_raster:
+                template_params['rhs'] = 'ST_Polygon(%s)' % template_params['rhs']
+        elif self.raster == BILATERAL:
+            # Operators with raster support but don't support mixed (rast-geom)
+            # lookups.
+            if lhs_is_raster and not rhs_is_raster:
+                template_params['lhs'] = 'ST_Polygon(%s)' % template_params['lhs']
+            elif rhs_is_raster and not lhs_is_raster:
+                template_params['rhs'] = 'ST_Polygon(%s)' % template_params['rhs']
+
+        return template_params
 
 
 class PostGISDistanceOperator(PostGISOperator):
@@ -35,6 +88,7 @@ class PostGISDistanceOperator(PostGISOperator):
 
     def as_sql(self, connection, lookup, template_params, sql_params):
         if not lookup.lhs.output_field.geography and lookup.lhs.output_field.geodetic(connection):
+            template_params = self.check_raster(lookup, template_params)
             sql_template = self.sql_template
             if len(lookup.rhs) == 3 and lookup.rhs[-1] == 'spheroid':
                 template_params.update({'op': self.op, 'func': 'ST_Distance_Spheroid'})
@@ -58,33 +112,33 @@ class PostGISOperations(BaseSpatialOperations, DatabaseOperations):
     Adapter = PostGISAdapter
 
     gis_operators = {
-        'bbcontains': PostGISOperator(op='~'),
-        'bboverlaps': PostGISOperator(op='&&', geography=True),
-        'contained': PostGISOperator(op='@'),
-        'contains': PostGISOperator(func='ST_Contains'),
-        'overlaps_left': PostGISOperator(op='&<'),
-        'overlaps_right': PostGISOperator(op='&>'),
+        'bbcontains': PostGISOperator(op='~', raster=True),
+        'bboverlaps': PostGISOperator(op='&&', geography=True, raster=True),
+        'contained': PostGISOperator(op='@', raster=True),
+        'overlaps_left': PostGISOperator(op='&<', raster=BILATERAL),
+        'overlaps_right': PostGISOperator(op='&>', raster=BILATERAL),
         'overlaps_below': PostGISOperator(op='&<|'),
         'overlaps_above': PostGISOperator(op='|&>'),
         'left': PostGISOperator(op='<<'),
         'right': PostGISOperator(op='>>'),
         'strictly_below': PostGISOperator(op='<<|'),
         'strictly_above': PostGISOperator(op='|>>'),
-        'same_as': PostGISOperator(op='~='),
-        'exact': PostGISOperator(op='~='),  # alias of same_as
-        'contains_properly': PostGISOperator(func='ST_ContainsProperly'),
-        'coveredby': PostGISOperator(func='ST_CoveredBy', geography=True),
-        'covers': PostGISOperator(func='ST_Covers', geography=True),
+        'same_as': PostGISOperator(op='~=', raster=BILATERAL),
+        'exact': PostGISOperator(op='~=', raster=BILATERAL),  # alias of same_as
+        'contains': PostGISOperator(func='ST_Contains', raster=BILATERAL),
+        'contains_properly': PostGISOperator(func='ST_ContainsProperly', raster=BILATERAL),
+        'coveredby': PostGISOperator(func='ST_CoveredBy', geography=True, raster=BILATERAL),
+        'covers': PostGISOperator(func='ST_Covers', geography=True, raster=BILATERAL),
         'crosses': PostGISOperator(func='ST_Crosses'),
-        'disjoint': PostGISOperator(func='ST_Disjoint'),
+        'disjoint': PostGISOperator(func='ST_Disjoint', raster=BILATERAL),
         'equals': PostGISOperator(func='ST_Equals'),
-        'intersects': PostGISOperator(func='ST_Intersects', geography=True),
+        'intersects': PostGISOperator(func='ST_Intersects', geography=True, raster=BILATERAL),
         'isvalid': PostGISOperator(func='ST_IsValid'),
-        'overlaps': PostGISOperator(func='ST_Overlaps'),
+        'overlaps': PostGISOperator(func='ST_Overlaps', raster=BILATERAL),
         'relate': PostGISOperator(func='ST_Relate'),
-        'touches': PostGISOperator(func='ST_Touches'),
-        'within': PostGISOperator(func='ST_Within'),
-        'dwithin': PostGISOperator(func='ST_DWithin', geography=True),
+        'touches': PostGISOperator(func='ST_Touches', raster=BILATERAL),
+        'within': PostGISOperator(func='ST_Within', raster=BILATERAL),
+        'dwithin': PostGISOperator(func='ST_DWithin', geography=True, raster=BILATERAL),
         'distance_gt': PostGISDistanceOperator(func='ST_Distance', op='>', geography=True),
         'distance_gte': PostGISDistanceOperator(func='ST_Distance', op='>=', geography=True),
         'distance_lt': PostGISDistanceOperator(func='ST_Distance', op='<', geography=True),
@@ -272,14 +326,14 @@ class PostGISOperations(BaseSpatialOperations, DatabaseOperations):
 
     def get_geom_placeholder(self, f, value, compiler):
         """
-        Provides a proper substitution value for Geometries that are not in the
-        SRID of the field.  Specifically, this routine will substitute in the
-        ST_Transform() function call.
+        Provide a proper substitution value for Geometries or rasters that are
+        not in the SRID of the field. Specifically, this routine will
+        substitute in the ST_Transform() function call.
         """
         # Get the srid for this object
         if value is None:
             value_srid = None
-        elif f.geom_type == 'RASTER':
+        elif f.geom_type == 'RASTER' and isinstance(value, six.string_types):
             value_srid = get_pgraster_srid(value)
         else:
             value_srid = value.srid
@@ -288,7 +342,7 @@ class PostGISOperations(BaseSpatialOperations, DatabaseOperations):
         # is not equal to the field srid.
         if value_srid is None or value_srid == f.srid:
             placeholder = '%s'
-        elif f.geom_type == 'RASTER':
+        elif f.geom_type == 'RASTER' and isinstance(value, six.string_types):
             placeholder = '%s((%%s)::raster, %s)' % (self.transform, f.srid)
         else:
             placeholder = '%s(%%s, %s)' % (self.transform, f.srid)
diff --git a/django/contrib/gis/db/models/fields.py b/django/contrib/gis/db/models/fields.py
index 8452a11c5f..91a12623c7 100644
--- a/django/contrib/gis/db/models/fields.py
+++ b/django/contrib/gis/db/models/fields.py
@@ -1,7 +1,10 @@
 from django.contrib.gis import forms
-from django.contrib.gis.db.models.lookups import gis_lookups
+from django.contrib.gis.db.models.lookups import (
+    RasterBandTransform, gis_lookups,
+)
 from django.contrib.gis.db.models.proxy import SpatialProxy
 from django.contrib.gis.gdal import HAS_GDAL
+from django.contrib.gis.gdal.error import GDALException
 from django.contrib.gis.geometry.backend import Geometry, GeometryException
 from django.core.exceptions import ImproperlyConfigured
 from django.db.models.expressions import Expression
@@ -157,6 +160,82 @@ class BaseSpatialField(Field):
         """
         return connection.ops.get_geom_placeholder(self, value, compiler)
 
+    def get_srid(self, obj):
+        """
+        Return the default SRID for the given geometry or raster, taking into
+        account the SRID set for the field. For example, if the input geometry
+        or raster doesn't have an SRID, then the SRID of the field will be
+        returned.
+        """
+        srid = obj.srid  # SRID of given geometry.
+        if srid is None or self.srid == -1 or (srid == -1 and self.srid != -1):
+            return self.srid
+        else:
+            return srid
+
+    def get_db_prep_save(self, value, connection):
+        """
+        Prepare the value for saving in the database.
+        """
+        if not value:
+            return None
+        else:
+            return connection.ops.Adapter(self.get_prep_value(value))
+
+    def get_prep_value(self, value):
+        """
+        Spatial lookup values are either a parameter that is (or may be
+        converted to) a geometry or raster, or a sequence of lookup values
+        that begins with a geometry or raster. This routine sets up the
+        geometry or raster value properly and preserves any other lookup
+        parameters.
+        """
+        from django.contrib.gis.gdal import GDALRaster
+
+        value = super(BaseSpatialField, self).get_prep_value(value)
+        # For IsValid lookups, boolean values are allowed.
+        if isinstance(value, (Expression, bool)):
+            return value
+        elif isinstance(value, (tuple, list)):
+            obj = value[0]
+            seq_value = True
+        else:
+            obj = value
+            seq_value = False
+
+        # When the input is not a geometry or raster, attempt to construct one
+        # from the given string input.
+        if isinstance(obj, (Geometry, GDALRaster)):
+            pass
+        elif isinstance(obj, (bytes, six.string_types)) or hasattr(obj, '__geo_interface__'):
+            try:
+                obj = Geometry(obj)
+            except (GeometryException, GDALException):
+                try:
+                    obj = GDALRaster(obj)
+                except GDALException:
+                    raise ValueError("Couldn't create spatial object from lookup value '%s'." % obj)
+        elif isinstance(obj, dict):
+            try:
+                obj = GDALRaster(obj)
+            except GDALException:
+                raise ValueError("Couldn't create spatial object from lookup value '%s'." % obj)
+        else:
+            raise ValueError('Cannot use object with type %s for a spatial lookup parameter.' % type(obj).__name__)
+
+        # Assigning the SRID value.
+        obj.srid = self.get_srid(obj)
+
+        if seq_value:
+            lookup_val = [obj]
+            lookup_val.extend(value[1:])
+            return tuple(lookup_val)
+        else:
+            return obj
+
+for klass in gis_lookups.values():
+    BaseSpatialField.register_lookup(klass)
+
 
 class GeometryField(GeoSelectFormatMixin, BaseSpatialField):
     """
@@ -224,6 +303,8 @@ class GeometryField(GeoSelectFormatMixin, BaseSpatialField):
         value properly, and preserve any other lookup parameters before
         returning to the caller.
         """
+        from django.contrib.gis.gdal import GDALRaster
+
         value = super(GeometryField, self).get_prep_value(value)
         if isinstance(value, (Expression, bool)):
             return value
@@ -236,7 +317,7 @@ class GeometryField(GeoSelectFormatMixin, BaseSpatialField):
 
         # When the input is not a GEOS geometry, attempt to construct one
         # from the given string input.
-        if isinstance(geom, Geometry):
+        if isinstance(geom, (Geometry, GDALRaster)):
             pass
         elif isinstance(geom, (bytes, six.string_types)) or hasattr(geom, '__geo_interface__'):
             try:
@@ -265,18 +346,6 @@ class GeometryField(GeoSelectFormatMixin, BaseSpatialField):
                 value.srid = self.srid
         return value
 
-    def get_srid(self, geom):
-        """
-        Returns the default SRID for the given geometry, taking into account
-        the SRID set for the field.  For example, if the input geometry
-        has no SRID, then that of the field will be returned.
-        """
-        gsrid = geom.srid  # SRID of given geometry.
-        if gsrid is None or self.srid == -1 or (gsrid == -1 and self.srid != -1):
-            return self.srid
-        else:
-            return gsrid
-
     # ### Routines overloaded from Field ###
     def contribute_to_class(self, cls, name, **kwargs):
         super(GeometryField, self).contribute_to_class(cls, name, **kwargs)
@@ -316,17 +385,6 @@ class GeometryField(GeoSelectFormatMixin, BaseSpatialField):
             params = [connection.ops.Adapter(value)]
         return params
 
-    def get_db_prep_save(self, value, connection):
-        "Prepares the value for saving in the database."
-        if not value:
-            return None
-        else:
-            return connection.ops.Adapter(self.get_prep_value(value))
-
-
-for klass in gis_lookups.values():
-    GeometryField.register_lookup(klass)
-
 
 # The OpenGIS Geometry Type Fields
 class PointField(GeometryField):
@@ -387,6 +445,7 @@ class RasterField(BaseSpatialField):
 
     description = _("Raster Field")
     geom_type = 'RASTER'
+    geography = False
 
     def __init__(self, *args, **kwargs):
         if not HAS_GDAL:
@@ -421,3 +480,15 @@ class RasterField(BaseSpatialField):
         # delays the instantiation of the objects to the moment of evaluation
         # of the raster attribute.
         setattr(cls, self.attname, SpatialProxy(GDALRaster, self))
+
+    def get_transform(self, name):
+        try:
+            band_index = int(name)
+            return type(
+                'SpecificRasterBandTransform',
+                (RasterBandTransform, ),
+                {'band_index': band_index}
+            )
+        except ValueError:
+            pass
+        return super(RasterField, self).get_transform(name)
diff --git a/django/contrib/gis/db/models/lookups.py b/django/contrib/gis/db/models/lookups.py
index 158fb5cfa1..a5bf428aa4 100644
--- a/django/contrib/gis/db/models/lookups.py
+++ b/django/contrib/gis/db/models/lookups.py
@@ -5,16 +5,23 @@ import re
 from django.core.exceptions import FieldDoesNotExist
 from django.db.models.constants import LOOKUP_SEP
 from django.db.models.expressions import Col, Expression
-from django.db.models.lookups import BuiltinLookup, Lookup
+from django.db.models.lookups import BuiltinLookup, Lookup, Transform
 from django.utils import six
 
 gis_lookups = {}
 
 
+class RasterBandTransform(Transform):
+    def as_sql(self, compiler, connection):
+        return compiler.compile(self.lhs)
+
+
 class GISLookup(Lookup):
     sql_template = None
     transform_func = None
     distance = False
+    band_rhs = None
+    band_lhs = None
 
     def __init__(self, *args, **kwargs):
         super(GISLookup, self).__init__(*args, **kwargs)
@@ -28,10 +35,10 @@ class GISLookup(Lookup):
         'point, 'the_geom', or a related lookup on a geographic field like
         'address__point'.
 
-        If a GeometryField exists according to the given lookup on the model
-        options, it will be returned.  Otherwise returns None.
+        If a BaseSpatialField exists according to the given lookup on the model
+        options, it will be returned. Otherwise return None.
         """
-        from django.contrib.gis.db.models.fields import GeometryField
+        from django.contrib.gis.db.models.fields import BaseSpatialField
         # This takes into account the situation where the lookup is a
         # lookup to a related geographic field, e.g., 'address__point'.
         field_list = lookup.split(LOOKUP_SEP)
@@ -55,11 +62,34 @@ class GISLookup(Lookup):
             return False
 
         # Finally, make sure we got a Geographic field and return.
-        if isinstance(geo_fld, GeometryField):
+        if isinstance(geo_fld, BaseSpatialField):
             return geo_fld
         else:
             return False
 
+    def process_band_indices(self, only_lhs=False):
+        """
+        Extract the lhs band index from the band transform class and the rhs
+        band index from the input tuple.
+        """
+        # PostGIS band indices are 1-based, so the band index needs to be
+        # increased to be consistent with the GDALRaster band indices.
+        if only_lhs:
+            self.band_rhs = 1
+            self.band_lhs = self.lhs.band_index + 1
+            return
+
+        if isinstance(self.lhs, RasterBandTransform):
+            self.band_lhs = self.lhs.band_index + 1
+        else:
+            self.band_lhs = 1
+
+        self.band_rhs = self.rhs[1]
+        if len(self.rhs) == 1:
+            self.rhs = self.rhs[0]
+        else:
+            self.rhs = (self.rhs[0], ) + self.rhs[2:]
+
     def get_db_prep_lookup(self, value, connection):
         # get_db_prep_lookup is called by process_rhs from super class
         if isinstance(value, (tuple, list)):
@@ -70,10 +100,9 @@ class GISLookup(Lookup):
         return ('%s', params)
 
     def process_rhs(self, compiler, connection):
-        rhs, rhs_params = super(GISLookup, self).process_rhs(compiler, connection)
         if hasattr(self.rhs, '_as_sql'):
             # If rhs is some QuerySet, don't touch it
-            return rhs, rhs_params
+            return super(GISLookup, self).process_rhs(compiler, connection)
 
         geom = self.rhs
         if isinstance(self.rhs, Col):
@@ -85,9 +114,19 @@ class GISLookup(Lookup):
                 raise ValueError('No geographic field found in expression.')
             self.rhs.srid = geo_fld.srid
         elif isinstance(self.rhs, Expression):
-            raise ValueError('Complex expressions not supported for GeometryField')
+            raise ValueError('Complex expressions not supported for spatial fields.')
         elif isinstance(self.rhs, (list, tuple)):
             geom = self.rhs[0]
+            # Check if a band index was passed in the query argument.
+            if ((len(self.rhs) == 2 and not self.lookup_name == 'relate') or
+                    (len(self.rhs) == 3 and self.lookup_name == 'relate')):
+                self.process_band_indices()
+            elif len(self.rhs) > 2:
+                raise ValueError('Tuple too long for lookup %s.' % self.lookup_name)
+        elif isinstance(self.lhs, RasterBandTransform):
+            self.process_band_indices(only_lhs=True)
+
+        rhs, rhs_params = super(GISLookup, self).process_rhs(compiler, connection)
         rhs = connection.ops.get_geom_placeholder(self.lhs.output_field, geom, compiler)
         return rhs, rhs_params
 
@@ -274,6 +313,8 @@ class IsValidLookup(BuiltinLookup):
     lookup_name = 'isvalid'
 
     def as_sql(self, compiler, connection):
+        if self.lhs.field.geom_type == 'RASTER':
+            raise ValueError('The isvalid lookup is only available on geometry fields.')
         gis_op = connection.ops.gis_operators[self.lookup_name]
         sql, params = self.process_lhs(compiler, connection)
         sql = '%(func)s(%(lhs)s)' % {'func': gis_op.func, 'lhs': sql}
@@ -323,9 +364,17 @@ class DistanceLookupBase(GISLookup):
     sql_template = '%(func)s(%(lhs)s, %(rhs)s) %(op)s %(value)s'
 
     def process_rhs(self, compiler, connection):
-        if not isinstance(self.rhs, (tuple, list)) or not 2 <= len(self.rhs) <= 3:
-            raise ValueError("2 or 3-element tuple required for '%s' lookup." % self.lookup_name)
+        if not isinstance(self.rhs, (tuple, list)) or not 2 <= len(self.rhs) <= 4:
+            raise ValueError("2, 3, or 4-element tuple required for '%s' lookup." % self.lookup_name)
+        elif len(self.rhs) == 4 and not self.rhs[3] == 'spheroid':
+            raise ValueError("For 4-element tuples the last argument must be the 'speroid' directive.")
+
+        # Check if the second parameter is a band index.
+        if len(self.rhs) > 2 and not self.rhs[2] == 'spheroid':
+            self.process_band_indices()
+
         params = [connection.ops.Adapter(self.rhs[0])]
+
         # Getting the distance parameter in the units of the field.
         dist_param = self.rhs[1]
         if hasattr(dist_param, 'resolve_expression'):
diff --git a/docs/ref/contrib/gis/db-api.txt b/docs/ref/contrib/gis/db-api.txt
index ef26f8f75e..4403fff392 100644
--- a/docs/ref/contrib/gis/db-api.txt
+++ b/docs/ref/contrib/gis/db-api.txt
@@ -53,7 +53,12 @@ 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.
+lookups are available for raster fields, but spatial database functions and
+aggregates aren't implemented for raster fields.
+
+.. versionchanged:: 1.10
+
+    ``RasterField`` now supports spatial lookups.
 
 Creating and Saving Models with Geometry Fields
 ===============================================
@@ -136,11 +141,20 @@ Spatial Lookups
 
 GeoDjango's lookup types may be used with any manager method like
 ``filter()``, ``exclude()``, etc.  However, the lookup types unique to
-GeoDjango are only available on geometry fields.
+GeoDjango are only available on spatial fields.
+
 Filters on 'normal' fields (e.g. :class:`~django.db.models.CharField`)
-may be chained with those on geographic fields.  Thus, geographic queries
-take the following general form (assuming  the ``Zipcode`` model used in the
-:doc:`model-api`)::
+may be chained with those on geographic fields. Geographic lookups accept
+geometry and raster input on both sides and input types can be mixed freely.
+
+The general structure of geographic lookups is described below. A complete
+reference can be found in the :ref:`spatial lookup reference<spatial-lookups>`.
+
+Geometry Lookups
+----------------
+
+Geographic queries with geometries take the following general form (assuming
+the ``Zipcode`` model used in the :doc:`model-api`)::
 
     >>> qs = Zipcode.objects.filter(<field>__<lookup_type>=<parameter>)
     >>> qs = Zipcode.objects.exclude(...)
@@ -148,14 +162,60 @@ take the following general form (assuming  the ``Zipcode`` model used in the
 For example::
 
     >>> qs = Zipcode.objects.filter(poly__contains=pnt)
+    >>> qs = Elevation.objects.filter(poly__contains=rst)
 
 In this case, ``poly`` is the geographic field, :lookup:`contains <gis-contains>`
-is the spatial lookup type, and ``pnt`` is the parameter (which may be a
+is the spatial lookup type, ``pnt`` is the parameter (which may be a
 :class:`~django.contrib.gis.geos.GEOSGeometry` object or a string of
-GeoJSON , WKT, or HEXEWKB).
+GeoJSON , WKT, or HEXEWKB), and ``rst`` is a
+:class:`~django.contrib.gis.gdal.GDALRaster` object.
+
+.. _spatial-lookup-raster:
+
+Raster Lookups
+--------------
+
+.. versionadded:: 1.10
 
-A complete reference can be found in the :ref:`spatial lookup reference
-<spatial-lookups>`.
+The raster lookup syntax is similar to the syntax for geometries. The only
+difference is that a band index can specified as additional input. If no band
+index is specified, the first band is used by default (index ``0``). In that
+case the syntax is identical to the syntax for geometry lookups.
+
+To specify the band index, an additional parameter can be specified on both
+sides of the lookup. On the left hand side, the double underscore syntax is
+used to pass a band index. On the right hand side, a tuple of the raster and
+band index can be specified.
+
+This results in the following general form for lookups involving rasters
+(assuming the ``Elevation`` model used in the :doc:`model-api`)::
+
+    >>> qs = Elevation.objects.filter(<field>__<lookup_type>=<parameter>)
+    >>> qs = Elevation.objects.filter(<field>__<band_index>__<lookup_type>=<parameter>)
+    >>> qs = Elevation.objects.filter(<field>__<lookup_type>=(<raster_input, <band_index>)
+
+Fore example::
+
+    >>> qs = Elevation.objects.filter(rast__contains=geom)
+    >>> qs = Elevation.objects.filter(rast__contains=rst)
+    >>> qs = Elevation.objects.filter(rast__1__contains=geom)
+    >>> qs = Elevation.objects.filter(rast__contains=(rst, 1))
+    >>> qs = Elevation.objects.filter(rast__1__contains=(rst, 1))
+
+On the left hand side of the example, ``rast`` is the geographic raster field
+and :lookup:`contains <gis-contains>` is the spatial lookup type. On the right
+hand side, ``geom`` is a geometry input and ``rst`` is a
+:class:`~django.contrib.gis.gdal.GDALRaster` object. The band index defaults to
+``0`` in the first two queries and is set to ``1`` on the others.
+
+While all spatial lookups can be used with raster objects on both sides, not all
+underlying operators natively accept raster input. For cases where the operator
+expects geometry input, the raster is automatically converted to a geometry.
+It's important to keep this in mind when interpreting the lookup results.
+
+The type of raster support is listed for all lookups in the :ref:`compatibility
+table <spatial-lookup-compatibility>`. Lookups involving rasters are currently
+only available for the PostGIS backend.
 
 .. _distance-queries:
 
@@ -176,7 +236,7 @@ in the :doc:`model-api` documentation for more details.
 Distance Lookups
 ----------------
 
-*Availability*: PostGIS, Oracle, SpatiaLite
+*Availability*: PostGIS, Oracle, SpatiaLite, PGRaster (Native)
 
 The following distance lookups are available:
 
@@ -193,7 +253,7 @@ The following distance lookups are available:
 
 Distance lookups take a tuple parameter comprising:
 
-#. A geometry to base calculations from; and
+#. A geometry or raster to base calculations from; and
 #. A number or :class:`~django.contrib.gis.measure.Distance` object containing the distance.
 
 If a :class:`~django.contrib.gis.measure.Distance` object is used,
@@ -241,6 +301,16 @@ Then distance queries may be performed as follows::
     >>> qs = SouthTexasCity.objects.filter(point__distance_gte=(pnt, D(mi=20)))
     >>> qs = SouthTexasCity.objects.filter(point__distance_gte=(pnt, D(chain=100)))
 
+Raster queries work the same way by simply replacing the geometry field
+``point`` with a raster field, or the ``pnt`` object with a raster object, or
+both. To specify the band index of a raster input on the right hand side, a
+3-tuple can be passed to the lookup as follows::
+
+    >>> qs = SouthTexasCity.objects.filter(point__distance_gte=(rst, 2, D(km=7)))
+
+Where the band with index 2 (the third band) of the raster ``rst`` would be
+used for the lookup.
+
 __ https://github.com/django/django/blob/master/tests/gis_tests/distapp/models.py
 
 .. _compatibility-table:
@@ -254,43 +324,46 @@ Spatial Lookups
 ---------------
 
 The following table provides a summary of what spatial lookups are available
-for each spatial database backend.
-
-=================================  =========  ========  ============ ==========
-Lookup Type                        PostGIS    Oracle    MySQL [#]_   SpatiaLite
-=================================  =========  ========  ============ ==========
-:lookup:`bbcontains`               X                    X            X
-:lookup:`bboverlaps`               X                    X            X
-:lookup:`contained`                X                    X            X
-:lookup:`contains <gis-contains>`  X          X         X            X
-:lookup:`contains_properly`        X
-:lookup:`coveredby`                X          X
-:lookup:`covers`                   X          X
-:lookup:`crosses`                  X                                 X
-:lookup:`disjoint`                 X          X         X            X
-:lookup:`distance_gt`              X          X                      X
-:lookup:`distance_gte`             X          X                      X
-:lookup:`distance_lt`              X          X                      X
-:lookup:`distance_lte`             X          X                      X
-:lookup:`dwithin`                  X          X
-:lookup:`equals`                   X          X         X            X
-:lookup:`exact`                    X          X         X            X
-:lookup:`intersects`               X          X         X            X
+for each spatial database backend. The PostGIS Raster (PGRaster) lookups are
+divided into the three categories described in the :ref:`raster lookup details
+<spatial-lookup-raster>`: native support ``N``, bilateral native support ``B``,
+and geometry conversion support ``C``.
+
+=================================  =========  ========  ============ ========== ========
+Lookup Type                        PostGIS    Oracle    MySQL [#]_   SpatiaLite PGRaster
+=================================  =========  ========  ============ ========== ========
+:lookup:`bbcontains`               X                    X            X          N
+:lookup:`bboverlaps`               X                    X            X          N
+:lookup:`contained`                X                    X            X          N
+:lookup:`contains <gis-contains>`  X          X         X            X          B
+:lookup:`contains_properly`        X                                            B
+:lookup:`coveredby`                X          X                                 B
+:lookup:`covers`                   X          X                                 B
+:lookup:`crosses`                  X                                 X          C
+:lookup:`disjoint`                 X          X         X            X          B
+:lookup:`distance_gt`              X          X                      X          N
+:lookup:`distance_gte`             X          X                      X          N
+:lookup:`distance_lt`              X          X                      X          N
+:lookup:`distance_lte`             X          X                      X          N
+:lookup:`dwithin`                  X          X                                 B
+:lookup:`equals`                   X          X         X            X          C
+:lookup:`exact`                    X          X         X            X          B
+:lookup:`intersects`               X          X         X            X          B
 :lookup:`isvalid`                  X
-:lookup:`overlaps`                 X          X         X            X
-:lookup:`relate`                   X          X                      X
-:lookup:`same_as`                  X          X         X            X
-:lookup:`touches`                  X          X         X            X
-:lookup:`within`                   X          X         X            X
-:lookup:`left`                     X
-:lookup:`right`                    X
-:lookup:`overlaps_left`            X
-:lookup:`overlaps_right`           X
-:lookup:`overlaps_above`           X
-:lookup:`overlaps_below`           X
-:lookup:`strictly_above`           X
-:lookup:`strictly_below`           X
-=================================  =========  ========  ============ ==========
+:lookup:`overlaps`                 X          X         X            X          B
+:lookup:`relate`                   X          X                      X          C
+:lookup:`same_as`                  X          X         X            X          B
+:lookup:`touches`                  X          X         X            X          B
+:lookup:`within`                   X          X         X            X          B
+:lookup:`left`                     X                                            C
+:lookup:`right`                    X                                            C
+:lookup:`overlaps_left`            X                                            B
+:lookup:`overlaps_right`           X                                            B
+:lookup:`overlaps_above`           X                                            C
+:lookup:`overlaps_below`           X                                            C
+:lookup:`strictly_above`           X                                            C
+:lookup:`strictly_below`           X                                            C
+=================================  =========  ========  ============ ========== ========
 
 .. _database-functions-compatibility:
 
diff --git a/docs/ref/contrib/gis/geoquerysets.txt b/docs/ref/contrib/gis/geoquerysets.txt
index 6b12c07a57..1c9d13df85 100644
--- a/docs/ref/contrib/gis/geoquerysets.txt
+++ b/docs/ref/contrib/gis/geoquerysets.txt
@@ -11,22 +11,70 @@ GeoQuerySet API Reference
 Spatial Lookups
 ===============
 
-The spatial lookups in this section are available for :class:`GeometryField`.
+The spatial lookups in this section are available for :class:`GeometryField`
+and :class:`RasterField`.
 
 For an introduction, see the :ref:`spatial lookups introduction
 <spatial-lookups-intro>`.  For an overview of what lookups are
 compatible with a particular spatial backend, refer to the
 :ref:`spatial lookup compatibility table <spatial-lookup-compatibility>`.
 
+.. versionchanged:: 1.10
+
+    Spatial lookups now support raster input.
+
+Lookups with rasters
+--------------------
+
+All examples in the reference below are given for geometry fields and inputs,
+but the lookups can be used the same way with rasters on both sides. Whenever
+a lookup doesn't support raster input, the input is automatically
+converted to a geometry where necessary using the `ST_Polygon
+<http://postgis.net/docs/RT_ST_Polygon.html>`_ function. See also the
+:ref:`introduction to raster lookups <spatial-lookup-raster>`.
+
+The database operators used by the lookups can be divided into three categories:
+
+- Native raster support ``N``: the operator accepts rasters natively on both
+  sides of the lookup, and raster input can be mixed with geometry inputs.
+
+- Bilateral raster support ``B``: the operator supports rasters only if both
+  sides of the lookup receive raster inputs. Raster data is automatically
+  converted to geometries for mixed lookups.
+
+- Geometry conversion support ``C``. The lookup does not have native raster
+  support, all raster data is automatically converted to geometries.
+
+The examples below show the SQL equivalent for the lookups in the different
+types of raster support. The same pattern applies to all spatial lookups.
+
+==== ============================== =======================================================
+Case Lookup                         SQL Equivalent
+==== ============================== =======================================================
+N, B ``rast__contains=rst``         ``ST_Contains(rast, rst)``
+N, B ``rast__1__contains=(rst, 2)`` ``ST_Contains(rast, 1, rst, 2)``
+B, C ``rast__contains=geom``        ``ST_Contains(ST_Polygon(rast), geom)``
+B, C ``rast__1__contains=geom``     ``ST_Contains(ST_Polygon(rast, 1), geom)``
+B, C ``poly__contains=rst``         ``ST_Contains(poly, ST_Polygon(rst))``
+B, C ``poly__contains=(rst, 1)``    ``ST_Contains(poly, ST_Polygon(rst, 1))``
+C    ``rast__crosses=rst``          ``ST_Crosses(ST_Polygon(rast), ST_Polygon(rst))``
+C    ``rast__1__crosses=(rst, 2)``  ``ST_Crosses(ST_Polygon(rast, 1), ST_Polygon(rst, 2))``
+C    ``rast__crosses=geom``         ``ST_Crosses(ST_Polygon(rast), geom)``
+C    ``poly__crosses=rst``          ``ST_Crosses(poly, ST_Polygon(rst))``
+==== ============================== =======================================================
+
+Spatial lookups with rasters are only supported for PostGIS backends
+(denominated as PGRaster in this section).
+
 .. fieldlookup:: bbcontains
 
 ``bbcontains``
 --------------
 
-*Availability*: PostGIS, MySQL, SpatiaLite
+*Availability*: PostGIS, MySQL, SpatiaLite, PGRaster (Native)
 
-Tests if the geometry field's bounding box completely contains the lookup
-geometry's bounding box.
+Tests if the geometry or raster field's bounding box completely contains the
+lookup geometry's bounding box.
 
 Example::
 
@@ -45,7 +93,7 @@ SpatiaLite  ``MbrContains(poly, geom)``
 ``bboverlaps``
 --------------
 
-*Availability*: PostGIS, MySQL, SpatiaLite
+*Availability*: PostGIS, MySQL, SpatiaLite, PGRaster (Native)
 
 Tests if the geometry field's bounding box overlaps the lookup geometry's
 bounding box.
@@ -67,7 +115,7 @@ SpatiaLite  ``MbrOverlaps(poly, geom)``
 ``contained``
 -------------
 
-*Availability*: PostGIS, MySQL, SpatiaLite
+*Availability*: PostGIS, MySQL, SpatiaLite, PGRaster (Native)
 
 Tests if the geometry field's bounding box is completely contained by the
 lookup geometry's bounding box.
@@ -89,7 +137,7 @@ SpatiaLite  ``MbrWithin(poly, geom)``
 ``contains``
 ------------
 
-*Availability*: PostGIS, Oracle, MySQL, SpatiaLite
+*Availability*: PostGIS, Oracle, MySQL, SpatiaLite, PGRaster (Bilateral)
 
 Tests if the geometry field spatially contains the lookup geometry.
 
@@ -111,7 +159,7 @@ SpatiaLite  ``Contains(poly, geom)``
 ``contains_properly``
 ---------------------
 
-*Availability*: PostGIS
+*Availability*: PostGIS, PGRaster (Bilateral)
 
 Returns true if the lookup geometry intersects the interior of the
 geometry field, but not the boundary (or exterior). [#fncontainsproperly]_
@@ -131,7 +179,7 @@ PostGIS     ``ST_ContainsProperly(poly, geom)``
 ``coveredby``
 -------------
 
-*Availability*: PostGIS, Oracle
+*Availability*: PostGIS, Oracle, PGRaster (Bilateral)
 
 Tests if no point in the geometry field is outside the lookup geometry.
 [#fncovers]_
@@ -152,7 +200,7 @@ Oracle      ``SDO_COVEREDBY(poly, geom)``
 ``covers``
 ----------
 
-*Availability*: PostGIS, Oracle
+*Availability*: PostGIS, Oracle, PGRaster (Bilateral)
 
 Tests if no point in the lookup geometry is outside the geometry field.
 [#fncovers]_
@@ -173,7 +221,7 @@ Oracle      ``SDO_COVERS(poly, geom)``
 ``crosses``
 -----------
 
-*Availability*: PostGIS, SpatiaLite
+*Availability*: PostGIS, SpatiaLite, PGRaster (Conversion)
 
 Tests if the geometry field spatially crosses the lookup geometry.
 
@@ -193,7 +241,7 @@ SpatiaLite  ``Crosses(poly, geom)``
 ``disjoint``
 ------------
 
-*Availability*: PostGIS, Oracle, MySQL, SpatiaLite
+*Availability*: PostGIS, Oracle, MySQL, SpatiaLite, PGRaster (Bilateral)
 
 Tests if the geometry field is spatially disjoint from the lookup geometry.
 
@@ -215,7 +263,7 @@ SpatiaLite  ``Disjoint(poly, geom)``
 ``equals``
 ----------
 
-*Availability*: PostGIS, Oracle, MySQL, SpatiaLite
+*Availability*: PostGIS, Oracle, MySQL, SpatiaLite, PGRaster (Conversion)
 
 .. fieldlookup:: exact
 .. fieldlookup:: same_as
@@ -223,14 +271,14 @@ SpatiaLite  ``Disjoint(poly, geom)``
 ``exact``, ``same_as``
 ----------------------
 
-*Availability*: PostGIS, Oracle, MySQL, SpatiaLite
+*Availability*: PostGIS, Oracle, MySQL, SpatiaLite, PGRaster (Bilateral)
 
 .. fieldlookup:: intersects
 
 ``intersects``
 --------------
 
-*Availability*: PostGIS, Oracle, MySQL, SpatiaLite
+*Availability*: PostGIS, Oracle, MySQL, SpatiaLite, PGRaster (Bilateral)
 
 Tests if the geometry field spatially intersects the lookup geometry.
 
@@ -271,14 +319,14 @@ PostGIS equivalent::
 ``overlaps``
 ------------
 
-*Availability*: PostGIS, Oracle, MySQL, SpatiaLite
+*Availability*: PostGIS, Oracle, MySQL, SpatiaLite, PGRaster (Bilateral)
 
 .. fieldlookup:: relate
 
 ``relate``
 ----------
 
-*Availability*: PostGIS, Oracle, SpatiaLite
+*Availability*: PostGIS, Oracle, SpatiaLite, PGRaster (Conversion)
 
 Tests if the geometry field is spatially related to the lookup geometry by
 the values given in the given pattern.  This lookup requires a tuple parameter,
@@ -293,7 +341,7 @@ The intersection pattern matrix may only use the following characters:
 ``1``, ``2``, ``T``, ``F``, or ``*``.  This lookup type allows users to "fine tune"
 a specific geometric relationship consistent with the DE-9IM model. [#fnde9im]_
 
-Example::
+Geometry example::
 
     # A tuple lookup parameter is used to specify the geometry and
     # the intersection pattern (the pattern here is for 'contains').
@@ -307,6 +355,16 @@ SpatiaLite SQL equivalent::
 
     SELECT ... WHERE Relate(poly, geom, 'T*T***FF*')
 
+Raster example::
+
+    Zipcode.objects.filter(poly__relate=(rast, 1, 'T*T***FF*'))
+    Zipcode.objects.filter(rast__2__relate=(rast, 1, 'T*T***FF*'))
+
+PostGIS SQL equivalent::
+
+    SELECT ... WHERE ST_Relate(poly, ST_Polygon(rast, 1), 'T*T***FF*')
+    SELECT ... WHERE ST_Relate(ST_Polygon(rast, 2), ST_Polygon(rast, 1), 'T*T***FF*')
+
 Oracle
 ~~~~~~
 
@@ -352,7 +410,7 @@ SpatiaLite  ``Touches(poly, geom)``
 ``within``
 ----------
 
-*Availability*: PostGIS, Oracle, MySQL, SpatiaLite
+*Availability*: PostGIS, Oracle, MySQL, SpatiaLite, PGRaster (Bilateral)
 
 Tests if the geometry field is spatially within the lookup geometry.
 
@@ -374,7 +432,7 @@ SpatiaLite  ``Within(poly, geom)``
 ``left``
 --------
 
-*Availability*: PostGIS
+*Availability*: PostGIS, PGRaster (Conversion)
 
 Tests if the geometry field's bounding box is strictly to the left of the
 lookup geometry's bounding box.
@@ -392,7 +450,7 @@ PostGIS equivalent::
 ``right``
 ---------
 
-*Availability*: PostGIS
+*Availability*: PostGIS, PGRaster (Conversion)
 
 Tests if the geometry field's bounding box is strictly to the right of the
 lookup geometry's bounding box.
@@ -410,7 +468,7 @@ PostGIS equivalent::
 ``overlaps_left``
 -----------------
 
-*Availability*: PostGIS
+*Availability*: PostGIS, PGRaster (Bilateral)
 
 Tests if the geometry field's bounding box overlaps or is to the left of the lookup
 geometry's bounding box.
@@ -429,7 +487,7 @@ PostGIS equivalent::
 ``overlaps_right``
 ------------------
 
-*Availability*: PostGIS
+*Availability*: PostGIS, PGRaster (Bilateral)
 
 Tests if the geometry field's bounding box overlaps or is to the right of the lookup
 geometry's bounding box.
@@ -447,7 +505,7 @@ PostGIS equivalent::
 ``overlaps_above``
 ------------------
 
-*Availability*: PostGIS
+*Availability*: PostGIS, PGRaster (Conversion)
 
 Tests if the geometry field's bounding box overlaps or is above the lookup
 geometry's bounding box.
@@ -465,7 +523,7 @@ PostGIS equivalent::
 ``overlaps_below``
 ------------------
 
-*Availability*: PostGIS
+*Availability*: PostGIS, PGRaster (Conversion)
 
 Tests if the geometry field's bounding box overlaps or is below the lookup
 geometry's bounding box.
@@ -483,7 +541,7 @@ PostGIS equivalent::
 ``strictly_above``
 ------------------
 
-*Availability*: PostGIS
+*Availability*: PostGIS, PGRaster (Conversion)
 
 Tests if the geometry field's bounding box is strictly above the lookup
 geometry's bounding box.
@@ -501,7 +559,7 @@ PostGIS equivalent::
 ``strictly_below``
 ------------------
 
-*Availability*: PostGIS
+*Availability*: PostGIS, PGRaster (Conversion)
 
 Tests if the geometry field's bounding box is strictly below the lookup
 geometry's bounding box.
@@ -520,27 +578,31 @@ PostGIS equivalent::
 Distance Lookups
 ================
 
-*Availability*: PostGIS, Oracle, SpatiaLite
+*Availability*: PostGIS, Oracle, SpatiaLite, PGRaster (Native)
 
 For an overview on performing distance queries, please refer to
 the :ref:`distance queries introduction <distance-queries>`.
 
 Distance lookups take the following form::
 
-    <field>__<distance lookup>=(<geometry>, <distance value>[, 'spheroid'])
+    <field>__<distance lookup>=(<geometry/raster>, <distance value>[, 'spheroid'])
+    <field>__<distance lookup>=(<raster>, <band_index>, <distance value>[, 'spheroid'])
+    <field>__<band_index>__<distance lookup>=(<raster>, <band_index>, <distance value>[, 'spheroid'])
 
 The value passed into a distance lookup is a tuple; the first two
 values are mandatory, and are the geometry to calculate distances to,
 and a distance value (either a number in units of the field, a
 :class:`~django.contrib.gis.measure.Distance` object, or a `query expression
-<ref/models/expressions>`).
+<ref/models/expressions>`). To pass a band index to the lookup, use a 3-tuple
+where the second entry is the band index.
 
 With PostGIS, on every distance lookup but :lookup:`dwithin`, an optional
-third element, ``'spheroid'``, may be included to tell GeoDjango
-to use the more accurate spheroid distance calculation functions on
-fields with a geodetic coordinate system (e.g., ``ST_Distance_Spheroid``
-would be used instead of ``ST_Distance_Sphere``). The simpler ``ST_Distance``
-function is used with projected coordinate systems.
+element, ``'spheroid'``, may be included to tell GeoDjango to use the more
+accurate spheroid distance calculation functions on fields with a geodetic
+coordinate system (e.g., ``ST_Distance_Spheroid`` would be used instead of
+``ST_Distance_Sphere``). The simpler ``ST_Distance`` function is used with
+projected coordinate systems. Rasters are converted to geometries for spheroid
+based lookups.
 
 .. versionadded:: 1.10
 
diff --git a/docs/releases/1.10.txt b/docs/releases/1.10.txt
index e750fb1d47..dd874e2ae8 100644
--- a/docs/releases/1.10.txt
+++ b/docs/releases/1.10.txt
@@ -160,6 +160,9 @@ Minor features
   :lookup:`isvalid` lookup, all for PostGIS. This allows filtering and
   repairing invalid geometries on the database side.
 
+* Added raster support for all :doc:`spatial lookups
+  </ref/contrib/gis/geoquerysets>`.
+
 :mod:`django.contrib.messages`
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/docs/spelling_wordlist b/docs/spelling_wordlist
index 5c3daa5caf..c9fb05f5ea 100644
--- a/docs/spelling_wordlist
+++ b/docs/spelling_wordlist
@@ -558,6 +558,7 @@ pessimization
 Petri
 Peucker
 pgAdmin
+PGRaster
 phishing
 php
 picklable
diff --git a/tests/gis_tests/gis_migrations/test_operations.py b/tests/gis_tests/gis_migrations/test_operations.py
index e2f5e106d7..aa6b13bdda 100644
--- a/tests/gis_tests/gis_migrations/test_operations.py
+++ b/tests/gis_tests/gis_migrations/test_operations.py
@@ -82,9 +82,9 @@ class OperationTests(TransactionTestCase):
         operation = migration_class(*args)
         new_state = project_state.clone()
         operation.state_forwards('gis', new_state)
+        self.current_state = new_state
         with connection.schema_editor() as editor:
             operation.database_forwards('gis', editor, project_state, new_state)
-        self.current_state = new_state
 
     def test_add_geom_field(self):
         """
diff --git a/tests/gis_tests/rasterapp/models.py b/tests/gis_tests/rasterapp/models.py
index 4f172d4470..81f039b8d6 100644
--- a/tests/gis_tests/rasterapp/models.py
+++ b/tests/gis_tests/rasterapp/models.py
@@ -3,6 +3,18 @@ from ..models import models
 
 class RasterModel(models.Model):
     rast = models.RasterField('A Verbose Raster Name', null=True, srid=4326, spatial_index=True, blank=True)
+    rastprojected = models.RasterField('A Projected Raster Table', srid=3086, null=True)
+    geom = models.PointField(null=True)
+
+    class Meta:
+        required_db_features = ['supports_raster']
+
+    def __str__(self):
+        return str(self.id)
+
+
+class RasterRelatedModel(models.Model):
+    rastermodel = models.ForeignKey(RasterModel, models.CASCADE)
 
     class Meta:
         required_db_features = ['supports_raster']
diff --git a/tests/gis_tests/rasterapp/test_rasterfield.py b/tests/gis_tests/rasterapp/test_rasterfield.py
index 4dbe019602..82bed69def 100644
--- a/tests/gis_tests/rasterapp/test_rasterfield.py
+++ b/tests/gis_tests/rasterapp/test_rasterfield.py
@@ -1,20 +1,48 @@
 import json
 
+from django.contrib.gis.db.models.lookups import (
+    DistanceLookupBase, gis_lookups,
+)
+from django.contrib.gis.gdal import HAS_GDAL
+from django.contrib.gis.geos import GEOSGeometry
+from django.contrib.gis.measure import D
 from django.contrib.gis.shortcuts import numpy
 from django.core.exceptions import ImproperlyConfigured
+from django.db.models import Q
 from django.test import (
     TestCase, TransactionTestCase, mock, skipUnlessDBFeature,
 )
 
 from ..data.rasters.textrasters import JSON_RASTER
 from ..models import models
-from .models import RasterModel
+from .models import RasterModel, RasterRelatedModel
+
+if HAS_GDAL:
+    from django.contrib.gis.gdal import GDALRaster
 
 
 @skipUnlessDBFeature('supports_raster')
 class RasterFieldTest(TransactionTestCase):
     available_apps = ['gis_tests.rasterapp']
 
+    def setUp(self):
+        rast = GDALRaster({
+            "srid": 4326,
+            "origin": [0, 0],
+            "scale": [-1, 1],
+            "skew": [0, 0],
+            "width": 5,
+            "height": 5,
+            "nr_of_bands": 2,
+            "bands": [{"data": range(25)}, {"data": range(25, 50)}],
+        })
+        model_instance = RasterModel.objects.create(
+            rast=rast,
+            rastprojected=rast,
+            geom="POINT (-95.37040 29.70486)",
+        )
+        RasterRelatedModel.objects.create(rastermodel=model_instance)
+
     def test_field_null_value(self):
         """
         Test creating a model where the RasterField has a null value.
@@ -89,6 +117,220 @@ class RasterFieldTest(TransactionTestCase):
             'A Verbose Raster Name'
         )
 
+    def test_all_gis_lookups_with_rasters(self):
+        """
+        Evaluate all possible lookups for all input combinations (i.e.
+        raster-raster, raster-geom, geom-raster) and for projected and
+        unprojected coordinate systems. This test just checks that the lookup
+        can be called, but doesn't check if the result makes logical sense.
+        """
+        from django.contrib.gis.db.backends.postgis.operations import PostGISOperations
+
+        # Create test raster and geom.
+        rast = GDALRaster(json.loads(JSON_RASTER))
+        stx_pnt = GEOSGeometry('POINT (-95.370401017314293 29.704867409475465)', 4326)
+        stx_pnt.transform(3086)
+
+        # Loop through all the GIS lookups.
+        for name, lookup in gis_lookups.items():
+            # Construct lookup filter strings.
+            combo_keys = [
+                field + name for field in [
+                    'rast__', 'rast__', 'rastprojected__0__', 'rast__',
+                    'rastprojected__', 'geom__', 'rast__',
+                ]
+            ]
+            if issubclass(lookup, DistanceLookupBase):
+                # Set lookup values for distance lookups.
+                combo_values = [
+                    (rast, 50, 'spheroid'),
+                    (rast, 0, 50, 'spheroid'),
+                    (rast, 0, D(km=1)),
+                    (stx_pnt, 0, 500),
+                    (stx_pnt, D(km=1000)),
+                    (rast, 500),
+                    (json.loads(JSON_RASTER), 500),
+                ]
+            elif name == 'relate':
+                # Set lookup values for the relate lookup.
+                combo_values = [
+                    (rast, 'T*T***FF*'),
+                    (rast, 0, 'T*T***FF*'),
+                    (rast, 0, 'T*T***FF*'),
+                    (stx_pnt, 0, 'T*T***FF*'),
+                    (stx_pnt, 'T*T***FF*'),
+                    (rast, 'T*T***FF*'),
+                    (json.loads(JSON_RASTER), 'T*T***FF*'),
+                ]
+            elif name == 'isvalid':
+                # The isvalid lookup doesn't make sense for rasters.
+                continue
+            elif PostGISOperations.gis_operators[name].func:
+                # Set lookup values for all function based operators.
+                combo_values = [
+                    rast, (rast, 0), (rast, 0), (stx_pnt, 0), stx_pnt,
+                    rast, rast, json.loads(JSON_RASTER)
+                ]
+            else:
+                # Override band lookup for these, as it's not supported.
+                combo_keys[2] = 'rastprojected__' + name
+                # Set lookup values for all other operators.
+                combo_values = [rast, rast, rast, stx_pnt, stx_pnt, rast, rast, json.loads(JSON_RASTER)]
+
+            # Create query filter combinations.
+            combos = [{x[0]: x[1]} for x in zip(combo_keys, combo_values)]
+
+            for combo in combos:
+                # Apply this query filter.
+                qs = RasterModel.objects.filter(**combo)
+
+                # Evaluate normal filter qs.
+                self.assertTrue(qs.count() in [0, 1])
+
+            # Evaluate on conditional Q expressions.
+            qs = RasterModel.objects.filter(Q(**combos[0]) & Q(**combos[1]))
+            self.assertTrue(qs.count() in [0, 1])
+
+    def test_dwithin_gis_lookup_ouptut_with_rasters(self):
+        """
+        Check the logical functionality of the dwithin lookup for different
+        input parameters.
+        """
+        # Create test raster and geom.
+        rast = GDALRaster(json.loads(JSON_RASTER))
+        stx_pnt = GEOSGeometry('POINT (-95.370401017314293 29.704867409475465)', 4326)
+        stx_pnt.transform(3086)
+
+        # Filter raster with different lookup raster formats.
+        qs = RasterModel.objects.filter(rastprojected__dwithin=(rast, D(km=1)))
+        self.assertEqual(qs.count(), 1)
+
+        qs = RasterModel.objects.filter(rastprojected__dwithin=(json.loads(JSON_RASTER), D(km=1)))
+        self.assertEqual(qs.count(), 1)
+
+        qs = RasterModel.objects.filter(rastprojected__dwithin=(JSON_RASTER, D(km=1)))
+        self.assertEqual(qs.count(), 1)
+
+        # Filter in an unprojected coordinate system.
+        qs = RasterModel.objects.filter(rast__dwithin=(rast, 40))
+        self.assertEqual(qs.count(), 1)
+
+        # Filter with band index transform.
+        qs = RasterModel.objects.filter(rast__1__dwithin=(rast, 1, 40))
+        self.assertEqual(qs.count(), 1)
+        qs = RasterModel.objects.filter(rast__1__dwithin=(rast, 40))
+        self.assertEqual(qs.count(), 1)
+        qs = RasterModel.objects.filter(rast__dwithin=(rast, 1, 40))
+        self.assertEqual(qs.count(), 1)
+
+        # Filter raster by geom.
+        qs = RasterModel.objects.filter(rast__dwithin=(stx_pnt, 500))
+        self.assertEqual(qs.count(), 1)
+
+        qs = RasterModel.objects.filter(rastprojected__dwithin=(stx_pnt, D(km=10000)))
+        self.assertEqual(qs.count(), 1)
+
+        qs = RasterModel.objects.filter(rast__dwithin=(stx_pnt, 5))
+        self.assertEqual(qs.count(), 0)
+
+        qs = RasterModel.objects.filter(rastprojected__dwithin=(stx_pnt, D(km=100)))
+        self.assertEqual(qs.count(), 0)
+
+        # Filter geom by raster.
+        qs = RasterModel.objects.filter(geom__dwithin=(rast, 500))
+        self.assertEqual(qs.count(), 1)
+
+        # Filter through related model.
+        qs = RasterRelatedModel.objects.filter(rastermodel__rast__dwithin=(rast, 40))
+        self.assertEqual(qs.count(), 1)
+
+        # Filter through related model with band index transform
+        qs = RasterRelatedModel.objects.filter(rastermodel__rast__1__dwithin=(rast, 40))
+        self.assertEqual(qs.count(), 1)
+
+        # Filter through conditional statements.
+        qs = RasterModel.objects.filter(Q(rast__dwithin=(rast, 40)) & Q(rastprojected__dwithin=(stx_pnt, D(km=10000))))
+        self.assertEqual(qs.count(), 1)
+
+        # Filter through different lookup.
+        qs = RasterModel.objects.filter(rastprojected__bbcontains=rast)
+        self.assertEqual(qs.count(), 1)
+
+    def test_lookup_input_tuple_too_long(self):
+        rast = GDALRaster(json.loads(JSON_RASTER))
+        qs = RasterModel.objects.filter(rast__bbcontains=(rast, 1, 2))
+        msg = 'Tuple too long for lookup bbcontains.'
+        with self.assertRaisesMessage(ValueError, msg):
+            qs.count()
+
+    def test_lookup_input_band_not_allowed(self):
+        rast = GDALRaster(json.loads(JSON_RASTER))
+        qs = RasterModel.objects.filter(rast__bbcontains=(rast, 1))
+        msg = 'Band indices are not allowed for this operator, it works on bbox only.'
+        with self.assertRaisesMessage(ValueError, msg):
+            qs.count()
+
+    def test_isvalid_lookup_with_raster_error(self):
+        qs = RasterModel.objects.filter(rast__isvalid=True)
+        msg = 'The isvalid lookup is only available on geometry fields.'
+        with self.assertRaisesMessage(ValueError, msg):
+            qs.count()
+
+    def test_result_of_gis_lookup_with_rasters(self):
+        # Point is in the interior
+        qs = RasterModel.objects.filter(rast__contains=GEOSGeometry('POINT (-0.5 0.5)', 4326))
+        self.assertEqual(qs.count(), 1)
+        # Point is in the exterior
+        qs = RasterModel.objects.filter(rast__contains=GEOSGeometry('POINT (0.5 0.5)', 4326))
+        self.assertEqual(qs.count(), 0)
+        # A point on the boundary is not contained properly
+        qs = RasterModel.objects.filter(rast__contains_properly=GEOSGeometry('POINT (0 0)', 4326))
+        self.assertEqual(qs.count(), 0)
+        # Raster is located left of the point
+        qs = RasterModel.objects.filter(rast__left=GEOSGeometry('POINT (1 0)', 4326))
+        self.assertEqual(qs.count(), 1)
+
+    def test_lookup_with_raster_bbox(self):
+        rast = GDALRaster(json.loads(JSON_RASTER))
+        # Shift raster upwards
+        rast.origin.y = 2
+        # The raster in the model is not strictly below
+        qs = RasterModel.objects.filter(rast__strictly_below=rast)
+        self.assertEqual(qs.count(), 0)
+        # Shift raster further upwards
+        rast.origin.y = 6
+        # The raster in the model is strictly below
+        qs = RasterModel.objects.filter(rast__strictly_below=rast)
+        self.assertEqual(qs.count(), 1)
+
+    def test_lookup_with_polygonized_raster(self):
+        rast = GDALRaster(json.loads(JSON_RASTER))
+        # Move raster to overlap with the model point on the left side
+        rast.origin.x = -95.37040 + 1
+        rast.origin.y = 29.70486
+        # Raster overlaps with point in model
+        qs = RasterModel.objects.filter(geom__intersects=rast)
+        self.assertEqual(qs.count(), 1)
+        # Change left side of raster to be nodata values
+        rast.bands[0].data(data=[0, 0, 0, 1, 1], shape=(5, 1))
+        rast.bands[0].nodata_value = 0
+        qs = RasterModel.objects.filter(geom__intersects=rast)
+        # Raster does not overlap anymore after polygonization
+        # where the nodata zone is not included.
+        self.assertEqual(qs.count(), 0)
+
+    def test_lookup_value_error(self):
+        # Test with invalid dict lookup parameter
+        obj = dict()
+        msg = "Couldn't create spatial object from lookup value '%s'." % obj
+        with self.assertRaisesMessage(ValueError, msg):
+            RasterModel.objects.filter(geom__intersects=obj)
+        # Test with invalid string lookup parameter
+        obj = '00000'
+        msg = "Couldn't create spatial object from lookup value '%s'." % obj
+        with self.assertRaisesMessage(ValueError, msg):
+            RasterModel.objects.filter(geom__intersects=obj)
+
 
 @mock.patch('django.contrib.gis.db.models.fields.HAS_GDAL', False)
 class RasterFieldWithoutGDALTest(TestCase):