From bbfad84dd980a97174c3b061a3d1b5f1373c380d Mon Sep 17 00:00:00 2001 From: Daniel Wiesmann Date: Thu, 21 Apr 2016 17:03:14 +0100 Subject: [PATCH] Fixed #25588 -- Added spatial lookups to RasterField. Thanks Tim Graham for the review. --- .../gis/db/backends/postgis/adapter.py | 44 +++- .../gis/db/backends/postgis/operations.py | 108 ++++++-- django/contrib/gis/db/models/fields.py | 121 +++++++-- django/contrib/gis/db/models/lookups.py | 69 ++++- docs/ref/contrib/gis/db-api.txt | 165 ++++++++---- docs/ref/contrib/gis/geoquerysets.txt | 132 +++++++--- docs/releases/1.10.txt | 3 + docs/spelling_wordlist | 1 + .../gis_migrations/test_operations.py | 2 +- tests/gis_tests/rasterapp/models.py | 12 + tests/gis_tests/rasterapp/test_rasterfield.py | 244 +++++++++++++++++- 11 files changed, 743 insertions(+), 158 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`. + +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(__=) >>> 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 ` -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. -A complete reference can be found in the :ref:`spatial lookup reference -`. +.. _spatial-lookup-raster: + +Raster Lookups +-------------- + +.. versionadded:: 1.10 + +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(__=) + >>> qs = Elevation.objects.filter(____=) + >>> qs = Elevation.objects.filter(__=() + +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 ` 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 `. 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. +for each spatial database backend. The PostGIS Raster (PGRaster) lookups are +divided into the three categories described in the :ref:`raster lookup details +`: native support ``N``, bilateral native support ``B``, +and geometry conversion support ``C``. -================================= ========= ======== ============ ========== -Lookup Type PostGIS Oracle MySQL [#]_ SpatiaLite -================================= ========= ======== ============ ========== -:lookup:`bbcontains` X X X -:lookup:`bboverlaps` X X X -:lookup:`contained` X X X -:lookup:`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 +================================= ========= ======== ============ ========== ======== +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 ` 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 `. For an overview of what lookups are compatible with a particular spatial backend, refer to the :ref:`spatial lookup compatibility table `. +.. 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 +`_ function. See also the +:ref:`introduction to raster lookups `. + +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 lookups take the following form:: - __=(, [, 'spheroid']) + __=(, [, 'spheroid']) + __=(, , [, 'spheroid']) + ____=(, , [, '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 -`). +`). 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 + `. + :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):