Fixed #24214 -- Added GIS functions to replace geoqueryset's methods

Thanks Simon Charette and Tim Graham for the reviews.

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

@@ -1,3 +1,4 @@
+import re
 from functools import partial
 
 from django.contrib.gis.db.models import aggregates
@@ -59,11 +60,11 @@ class BaseSpatialFeatures(object):
     # `has_<name>_method` (defined in __init__) which accesses connection.ops
     # to determine GIS method availability.
     geoqueryset_methods = (
-        'area', 'centroid', 'difference', 'distance', 'distance_spheroid',
+        'area', 'bounding_circle', 'centroid', 'difference', 'distance',
-        'envelope', 'force_rhr', 'geohash', 'gml', 'intersection', 'kml',
+        'distance_spheroid', 'envelope', 'force_rhr', 'geohash', 'gml',
-        'length', 'num_geom', 'perimeter', 'point_on_surface', 'reverse',
+        'intersection', 'kml', 'length', 'mem_size', 'num_geom', 'num_points',
-        'scale', 'snap_to_grid', 'svg', 'sym_difference', 'transform',
+        'perimeter', 'point_on_surface', 'reverse', 'scale', 'snap_to_grid',
-        'translate', 'union', 'unionagg',
+        'svg', 'sym_difference', 'transform', 'translate', 'union', 'unionagg',
     )
 
     # Specifies whether the Collect and Extent aggregates are supported by the database
@@ -86,5 +87,13 @@ class BaseSpatialFeatures(object):
             setattr(self.__class__, 'has_%s_method' % method,
                     property(partial(BaseSpatialFeatures.has_ops_method, method=method)))
 
+    def __getattr__(self, name):
+        m = re.match(r'has_(\w*)_function$', name)
+        if m:
+            func_name = m.group(1)
+            if func_name not in self.connection.ops.unsupported_functions:
+                return True
+        return False
+
     def has_ops_method(self, method):
         return getattr(self.connection.ops, method, False)

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

@@ -22,6 +22,7 @@ class BaseSpatialOperations(object):
     geometry = False
 
     area = False
+    bounding_circle = False
     centroid = False
     difference = False
     distance = False
@@ -30,7 +31,6 @@ class BaseSpatialOperations(object):
     envelope = False
     force_rhr = False
     mem_size = False
-    bounding_circle = False
     num_geom = False
     num_points = False
     perimeter = False
@@ -48,6 +48,22 @@ class BaseSpatialOperations(object):
     # Aggregates
     disallowed_aggregates = ()
 
+    geom_func_prefix = ''
+
+    # Mapping between Django function names and backend names, when names do not
+    # match; used in spatial_function_name().
+    function_names = {}
+
+    # Blacklist/set of known unsupported functions of the backend
+    unsupported_functions = {
+        'Area', 'AsGeoHash', 'AsGeoJSON', 'AsGML', 'AsKML', 'AsSVG',
+        'BoundingCircle', 'Centroid', 'Difference', 'Distance', 'Envelope',
+        'ForceRHR', 'Intersection', 'Length', 'MemSize', 'NumGeometries',
+        'NumPoints', 'Perimeter', 'PointOnSurface', 'Reverse', 'Scale',
+        'SnapToGrid', 'SymDifference', 'Transform', 'Translate',
+        'Union',
+    }
+
     # Serialization
     geohash = False
     geojson = False
@@ -108,9 +124,14 @@ class BaseSpatialOperations(object):
     def spatial_aggregate_name(self, agg_name):
         raise NotImplementedError('Aggregate support not implemented for this spatial backend.')
 
+    def spatial_function_name(self, func_name):
+        if func_name in self.unsupported_functions:
+            raise NotImplementedError("This backend doesn't support the %s function." % func_name)
+        return self.function_names.get(func_name, self.geom_func_prefix + func_name)
+
     # Routines for getting the OGC-compliant models.
     def geometry_columns(self):
-        raise NotImplementedError('subclasses of BaseSpatialOperations must a provide geometry_columns() method')
+        raise NotImplementedError('Subclasses of BaseSpatialOperations must provide a geometry_columns() method.')
 
     def spatial_ref_sys(self):
         raise NotImplementedError('subclasses of BaseSpatialOperations must a provide spatial_ref_sys() method')

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

@@ -8,12 +8,13 @@ from psycopg2.extensions import ISQLQuote
 
 
 class PostGISAdapter(object):
-    def __init__(self, geom):
+    def __init__(self, geom, geography=False):
         "Initializes on the 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
+        self.geography = geography
         self._adapter = Binary(self.ewkb)
 
     def __conform__(self, proto):
@@ -44,4 +45,7 @@ class PostGISAdapter(object):
     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('ST_GeomFromEWKB(%s)' % self._adapter.getquoted().decode())
+        return str('%s(%s)' % (
+            'ST_GeogFromWKB' if self.geography else 'ST_GeomFromEWKB',
+            self._adapter.getquoted().decode())
+        )

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

@@ -88,6 +88,13 @@ class PostGISOperations(BaseSpatialOperations, DatabaseOperations):
         'distance_lte': PostGISDistanceOperator(func='ST_Distance', op='<=', geography=True),
     }
 
+    unsupported_functions = set()
+    function_names = {
+        'BoundingCircle': 'ST_MinimumBoundingCircle',
+        'MemSize': 'ST_Mem_Size',
+        'NumPoints': 'ST_NPoints',
+    }
+
     def __init__(self, connection):
         super(PostGISOperations, self).__init__(connection)
 

django/contrib/gis/db/models/functions.py

@@ -0,0 +1,351 @@
+from decimal import Decimal
+
+from django.contrib.gis.db.models.fields import GeometryField
+from django.contrib.gis.db.models.sql import AreaField
+from django.contrib.gis.geos.geometry import GEOSGeometry
+from django.contrib.gis.measure import (
+    Area as AreaMeasure, Distance as DistanceMeasure,
+)
+from django.core.exceptions import FieldError
+from django.db.models import FloatField, IntegerField, TextField
+from django.db.models.expressions import Func, Value
+from django.utils import six
+
+NUMERIC_TYPES = six.integer_types + (float, Decimal)
+
+
+class GeoFunc(Func):
+    function = None
+    output_field_class = None
+    geom_param_pos = 0
+
+    def __init__(self, *expressions, **extra):
+        if 'output_field' not in extra and self.output_field_class:
+            extra['output_field'] = self.output_field_class()
+        super(GeoFunc, self).__init__(*expressions, **extra)
+
+    @property
+    def name(self):
+        return self.__class__.__name__
+
+    @property
+    def srid(self):
+        expr = self.source_expressions[self.geom_param_pos]
+        if hasattr(expr, 'srid'):
+            return expr.srid
+        try:
+            return expr.field.srid
+        except (AttributeError, FieldError):
+            return None
+
+    def as_sql(self, compiler, connection):
+        if self.function is None:
+            self.function = connection.ops.spatial_function_name(self.name)
+        return super(GeoFunc, self).as_sql(compiler, connection)
+
+    def resolve_expression(self, *args, **kwargs):
+        res = super(GeoFunc, self).resolve_expression(*args, **kwargs)
+        base_srid = res.srid
+        if not base_srid:
+            raise TypeError("Geometry functions can only operate on geometric content.")
+
+        for pos, expr in enumerate(res.source_expressions[1:], start=1):
+            if isinstance(expr, GeomValue) and expr.srid != base_srid:
+                # Automatic SRID conversion so objects are comparable
+                res.source_expressions[pos] = Transform(expr, base_srid).resolve_expression(*args, **kwargs)
+        return res
+
+    def _handle_param(self, value, param_name='', check_types=None):
+        if not hasattr(value, 'resolve_expression'):
+            if check_types and not isinstance(value, check_types):
+                raise TypeError(
+                    "The %s parameter has the wrong type: should be %s." % (
+                        param_name, str(check_types))
+                )
+        return value
+
+
+class GeomValue(Value):
+    geography = False
+
+    @property
+    def srid(self):
+        return self.value.srid
+
+    def as_sql(self, compiler, connection):
+        if self.geography:
+            self.value = connection.ops.Adapter(self.value, geography=self.geography)
+        else:
+            self.value = connection.ops.Adapter(self.value)
+        return super(GeomValue, self).as_sql(compiler, connection)
+
+
+class GeoFuncWithGeoParam(GeoFunc):
+    def __init__(self, expression, geom, *expressions, **extra):
+        if not hasattr(geom, 'srid'):
+            # Try to interpret it as a geometry input
+            try:
+                geom = GEOSGeometry(geom)
+            except Exception:
+                raise ValueError("This function requires a geometric parameter.")
+        if not geom.srid:
+            raise ValueError("Please provide a geometry attribute with a defined SRID.")
+        geom = GeomValue(geom)
+        super(GeoFuncWithGeoParam, self).__init__(expression, geom, *expressions, **extra)
+
+
+class Area(GeoFunc):
+    def as_sql(self, compiler, connection):
+        if connection.ops.oracle:
+            self.output_field = AreaField('sq_m')  # Oracle returns area in units of meters.
+        else:
+            if connection.ops.geography:
+                # Geography fields support area calculation, returns square meters.
+                self.output_field = AreaField('sq_m')
+            elif not self.output_field.geodetic(connection):
+                # Getting the area units of the geographic field.
+                self.output_field = AreaField(
+                    AreaMeasure.unit_attname(self.output_field.units_name(connection)))
+            else:
+                # TODO: Do we want to support raw number areas for geodetic fields?
+                raise NotImplementedError('Area on geodetic coordinate systems not supported.')
+        return super(Area, self).as_sql(compiler, connection)
+
+
+class AsGeoJSON(GeoFunc):
+    output_field_class = TextField
+
+    def __init__(self, expression, bbox=False, crs=False, precision=8, **extra):
+        expressions = [expression]
+        if precision is not None:
+            expressions.append(self._handle_param(precision, 'precision', six.integer_types))
+        options = 0
+        if crs and bbox:
+            options = 3
+        elif bbox:
+            options = 1
+        elif crs:
+            options = 2
+        if options:
+            expressions.append(options)
+        super(AsGeoJSON, self).__init__(*expressions, **extra)
+
+
+class AsGML(GeoFunc):
+    geom_param_pos = 1
+    output_field_class = TextField
+
+    def __init__(self, expression, version=2, precision=8, **extra):
+        expressions = [version, expression]
+        if precision is not None:
+            expressions.append(self._handle_param(precision, 'precision', six.integer_types))
+        super(AsGML, self).__init__(*expressions, **extra)
+
+
+class AsKML(AsGML):
+    pass
+
+
+class AsSVG(GeoFunc):
+    output_field_class = TextField
+
+    def __init__(self, expression, relative=False, precision=8, **extra):
+        relative = relative if hasattr(relative, 'resolve_expression') else int(relative)
+        expressions = [
+            expression,
+            relative,
+            self._handle_param(precision, 'precision', six.integer_types),
+        ]
+        super(AsSVG, self).__init__(*expressions, **extra)
+
+
+class BoundingCircle(GeoFunc):
+    def __init__(self, expression, num_seg=48, **extra):
+        super(BoundingCircle, self).__init__(*[expression, num_seg], **extra)
+
+
+class Centroid(GeoFunc):
+    pass
+
+
+class Difference(GeoFuncWithGeoParam):
+    pass
+
+
+class DistanceResultMixin(object):
+    def convert_value(self, value, expression, connection, context):
+        if value is None:
+            return None
+        geo_field = GeometryField(srid=self.srid)  # Fake field to get SRID info
+        if geo_field.geodetic(connection):
+            dist_att = 'm'
+        else:
+            dist_att = DistanceMeasure.unit_attname(geo_field.units_name(connection))
+        return DistanceMeasure(**{dist_att: value})
+
+
+class Distance(DistanceResultMixin, GeoFuncWithGeoParam):
+    output_field_class = FloatField
+    spheroid = None
+
+    def __init__(self, expr1, expr2, spheroid=None, **extra):
+        expressions = [expr1, expr2]
+        if spheroid is not None:
+            self.spheroid = spheroid
+            expressions += (self._handle_param(spheroid, 'spheroid', bool),)
+        super(Distance, self).__init__(*expressions, **extra)
+
+    def as_postgresql(self, compiler, connection):
+        geo_field = GeometryField(srid=self.srid)  # Fake field to get SRID info
+        src_field = self.get_source_fields()[0]
+        geography = src_field.geography and self.srid == 4326
+        if geography:
+            # Set parameters as geography if base field is geography
+            for pos, expr in enumerate(
+                    self.source_expressions[self.geom_param_pos + 1:], start=self.geom_param_pos + 1):
+                if isinstance(expr, GeomValue):
+                    expr.geography = True
+        elif geo_field.geodetic(connection):
+            # Geometry fields with geodetic (lon/lat) coordinates need special distance functions
+            if self.spheroid:
+                self.function = 'ST_Distance_Spheroid'  # More accurate, resource intensive
+                # Replace boolean param by the real spheroid of the base field
+                self.source_expressions[2] = Value(geo_field._spheroid)
+            else:
+                self.function = 'ST_Distance_Sphere'
+        return super(Distance, self).as_sql(compiler, connection)
+
+
+class Envelope(GeoFunc):
+    pass
+
+
+class ForceRHR(GeoFunc):
+    pass
+
+
+class GeoHash(GeoFunc):
+    output_field_class = TextField
+
+    def __init__(self, expression, precision=None, **extra):
+        expressions = [expression]
+        if precision is not None:
+            expressions.append(self._handle_param(precision, 'precision', six.integer_types))
+        super(GeoHash, self).__init__(*expressions, **extra)
+
+
+class Intersection(GeoFuncWithGeoParam):
+    pass
+
+
+class Length(DistanceResultMixin, GeoFunc):
+    output_field_class = FloatField
+
+    def __init__(self, expr1, spheroid=True, **extra):
+        self.spheroid = spheroid
+        super(Length, self).__init__(expr1, **extra)
+
+    def as_postgresql(self, compiler, connection):
+        geo_field = GeometryField(srid=self.srid)  # Fake field to get SRID info
+        src_field = self.get_source_fields()[0]
+        geography = src_field.geography and self.srid == 4326
+        if geography:
+            self.source_expressions.append(Value(self.spheroid))
+        elif geo_field.geodetic(connection):
+            # Geometry fields with geodetic (lon/lat) coordinates need length_spheroid
+            self.function = 'ST_Length_Spheroid'
+            self.source_expressions.append(Value(geo_field._spheroid))
+        else:
+            dim = min(f.dim for f in self.get_source_fields() if f)
+            if dim > 2:
+                self.function = connection.ops.length3d
+        return super(Length, self).as_sql(compiler, connection)
+
+
+class MemSize(GeoFunc):
+    output_field_class = IntegerField
+
+
+class NumGeometries(GeoFunc):
+    output_field_class = IntegerField
+
+
+class NumPoints(GeoFunc):
+    output_field_class = IntegerField
+
+
+class Perimeter(DistanceResultMixin, GeoFunc):
+    output_field_class = FloatField
+
+    def as_postgresql(self, compiler, connection):
+        dim = min(f.dim for f in self.get_source_fields())
+        if dim > 2:
+            self.function = connection.ops.perimeter3d
+        return super(Perimeter, self).as_sql(compiler, connection)
+
+
+class PointOnSurface(GeoFunc):
+    pass
+
+
+class Reverse(GeoFunc):
+    pass
+
+
+class Scale(GeoFunc):
+    def __init__(self, expression, x, y, z=0.0, **extra):
+        expressions = [
+            expression,
+            self._handle_param(x, 'x', NUMERIC_TYPES),
+            self._handle_param(y, 'y', NUMERIC_TYPES),
+        ]
+        if z != 0.0:
+            expressions.append(self._handle_param(z, 'z', NUMERIC_TYPES))
+        super(Scale, self).__init__(*expressions, **extra)
+
+
+class SnapToGrid(GeoFunc):
+    def __init__(self, expression, *args, **extra):
+        nargs = len(args)
+        expressions = [expression]
+        if nargs in (1, 2):
+            expressions.extend(
+                [self._handle_param(arg, '', NUMERIC_TYPES) for arg in args]
+            )
+        elif nargs == 4:
+            # Reverse origin and size param ordering
+            expressions.extend(
+                [self._handle_param(arg, '', NUMERIC_TYPES) for arg in args[2:]]
+            )
+            expressions.extend(
+                [self._handle_param(arg, '', NUMERIC_TYPES) for arg in args[0:2]]
+            )
+        else:
+            raise ValueError('Must provide 1, 2, or 4 arguments to `SnapToGrid`.')
+        super(SnapToGrid, self).__init__(*expressions, **extra)
+
+
+class SymDifference(GeoFuncWithGeoParam):
+    pass
+
+
+class Transform(GeoFunc):
+    def __init__(self, expression, srid, **extra):
+        expressions = [
+            expression,
+            self._handle_param(srid, 'srid', six.integer_types),
+        ]
+        super(Transform, self).__init__(*expressions, **extra)
+
+    @property
+    def srid(self):
+        # Make srid the resulting srid of the transformation
+        return self.source_expressions[self.geom_param_pos + 1].value
+
+
+class Translate(Scale):
+    pass
+
+
+class Union(GeoFuncWithGeoParam):
+    pass

tests/gis_tests/distapp/tests.py

@@ -1,5 +1,8 @@
 from __future__ import unicode_literals
 
+from django.contrib.gis.db.models.functions import (
+    Area, Distance, Length, Perimeter, Transform,
+)
 from django.contrib.gis.geos import HAS_GEOS
 from django.contrib.gis.measure import D  # alias for Distance
 from django.db import connection
@@ -390,3 +393,275 @@ class DistanceTest(TestCase):
             'distance'
         ).values_list('name', flat=True).filter(name__in=('San Antonio', 'Pearland'))
         self.assertQuerysetEqual(qs, ['San Antonio', 'Pearland'], lambda x: x)
+
+
+'''
+=============================
+Distance functions on PostGIS
+=============================
+
+                                              | Projected Geometry | Lon/lat Geometry | Geography (4326)
+
+ST_Distance(geom1, geom2)                     |    OK (meters)     |   :-( (degrees)  |    OK (meters)
+
+ST_Distance(geom1, geom2, use_spheroid=False) |    N/A             |   N/A            |    OK (meters), less accurate, quick
+
+Distance_Sphere(geom1, geom2)                 |    N/A             |   OK (meters)    |    N/A
+
+Distance_Spheroid(geom1, geom2, spheroid)     |    N/A             |   OK (meters)    |    N/A
+
+
+================================
+Distance functions on Spatialite
+================================
+
+                                                | Projected Geometry | Lon/lat Geometry
+
+ST_Distance(geom1, geom2)                       |    OK (meters)     |      N/A
+
+ST_Distance(geom1, geom2, use_ellipsoid=True)   |    N/A             |      OK (meters)
+
+ST_Distance(geom1, geom2, use_ellipsoid=False)  |    N/A             |      OK (meters), less accurate, quick
+
+'''
+
+
+@skipUnlessDBFeature("gis_enabled")
+class DistanceFunctionsTests(TestCase):
+    fixtures = ['initial']
+
+    @skipUnlessDBFeature("has_Area_function")
+    def test_area(self):
+        # Reference queries:
+        # SELECT ST_Area(poly) FROM distapp_southtexaszipcode;
+        area_sq_m = [5437908.90234375, 10183031.4389648, 11254471.0073242, 9881708.91772461]
+        # Tolerance has to be lower for Oracle
+        tol = 2
+        for i, z in enumerate(SouthTexasZipcode.objects.annotate(area=Area('poly')).order_by('name')):
+            self.assertAlmostEqual(area_sq_m[i], z.area.sq_m, tol)
+
+    @skipUnlessDBFeature("has_Distance_function")
+    def test_distance_simple(self):
+        """
+        Test a simple distance query, with projected coordinates and without
+        transformation.
+        """
+        lagrange = GEOSGeometry('POINT(805066.295722839 4231496.29461335)', 32140)
+        houston = SouthTexasCity.objects.annotate(dist=Distance('point', lagrange)).order_by('id').first()
+        tol = 2 if oracle else 5
+        self.assertAlmostEqual(
+            houston.dist.m if hasattr(houston.dist, 'm') else houston.dist,
+            147075.069813,
+            tol
+        )
+
+    @skipUnlessDBFeature("has_Distance_function", "has_Transform_function")
+    def test_distance_projected(self):
+        """
+        Test the `Distance` function on projected coordinate systems.
+        """
+        # The point for La Grange, TX
+        lagrange = GEOSGeometry('POINT(-96.876369 29.905320)', 4326)
+        # Reference distances in feet and in meters. Got these values from
+        # using the provided raw SQL statements.
+        #  SELECT ST_Distance(point, ST_Transform(ST_GeomFromText('POINT(-96.876369 29.905320)', 4326), 32140))
+        #  FROM distapp_southtexascity;
+        m_distances = [147075.069813, 139630.198056, 140888.552826,
+                       138809.684197, 158309.246259, 212183.594374,
+                       70870.188967, 165337.758878, 139196.085105]
+        #  SELECT ST_Distance(point, ST_Transform(ST_GeomFromText('POINT(-96.876369 29.905320)', 4326), 2278))
+        #  FROM distapp_southtexascityft;
+        # Oracle 11 thinks this is not a projected coordinate system, so it's
+        # not tested.
+        ft_distances = [482528.79154625, 458103.408123001, 462231.860397575,
+                        455411.438904354, 519386.252102563, 696139.009211594,
+                        232513.278304279, 542445.630586414, 456679.155883207]
+
+        # Testing using different variations of parameters and using models
+        # with different projected coordinate systems.
+        dist1 = SouthTexasCity.objects.annotate(distance=Distance('point', lagrange)).order_by('id')
+        if spatialite or oracle:
+            dist_qs = [dist1]
+        else:
+            dist2 = SouthTexasCityFt.objects.annotate(distance=Distance('point', lagrange)).order_by('id')
+            # Using EWKT string parameter.
+            dist3 = SouthTexasCityFt.objects.annotate(distance=Distance('point', lagrange.ewkt)).order_by('id')
+            dist_qs = [dist1, dist2, dist3]
+
+        # Original query done on PostGIS, have to adjust AlmostEqual tolerance
+        # for Oracle.
+        tol = 2 if oracle else 5
+
+        # Ensuring expected distances are returned for each distance queryset.
+        for qs in dist_qs:
+            for i, c in enumerate(qs):
+                self.assertAlmostEqual(m_distances[i], c.distance.m, tol)
+                self.assertAlmostEqual(ft_distances[i], c.distance.survey_ft, tol)
+
+    @skipUnlessDBFeature("has_Distance_function", "supports_distance_geodetic")
+    def test_distance_geodetic(self):
+        """
+        Test the `Distance` function on geodetic coordinate systems.
+        """
+        # Testing geodetic distance calculation with a non-point geometry
+        # (a LineString of Wollongong and Shellharbour coords).
+        ls = LineString(((150.902, -34.4245), (150.87, -34.5789)), srid=4326)
+
+        # Reference query:
+        #  SELECT ST_distance_sphere(point, ST_GeomFromText('LINESTRING(150.9020 -34.4245,150.8700 -34.5789)', 4326))
+        #  FROM distapp_australiacity ORDER BY name;
+        distances = [1120954.92533513, 140575.720018241, 640396.662906304,
+                     60580.9693849269, 972807.955955075, 568451.8357838,
+                     40435.4335201384, 0, 68272.3896586844, 12375.0643697706, 0]
+        qs = AustraliaCity.objects.annotate(distance=Distance('point', ls)).order_by('name')
+        for city, distance in zip(qs, distances):
+            # Testing equivalence to within a meter.
+            self.assertAlmostEqual(distance, city.distance.m, 0)
+
+    @skipUnlessDBFeature("has_Distance_function", "supports_distance_geodetic")
+    def test_distance_geodetic_spheroid(self):
+        tol = 2 if oracle else 5
+
+        # Got the reference distances using the raw SQL statements:
+        #  SELECT ST_distance_spheroid(point, ST_GeomFromText('POINT(151.231341 -33.952685)', 4326),
+        #    'SPHEROID["WGS 84",6378137.0,298.257223563]') FROM distapp_australiacity WHERE (NOT (id = 11));
+        #  SELECT ST_distance_sphere(point, ST_GeomFromText('POINT(151.231341 -33.952685)', 4326))
+        #  FROM distapp_australiacity WHERE (NOT (id = 11));  st_distance_sphere
+        if connection.ops.postgis and connection.ops.proj_version_tuple() >= (4, 7, 0):
+            # PROJ.4 versions 4.7+ have updated datums, and thus different
+            # distance values.
+            spheroid_distances = [60504.0628957201, 77023.9489850262, 49154.8867574404,
+                                  90847.4358768573, 217402.811919332, 709599.234564757,
+                                  640011.483550888, 7772.00667991925, 1047861.78619339,
+                                  1165126.55236034]
+            sphere_distances = [60580.9693849267, 77144.0435286473, 49199.4415344719,
+                                90804.7533823494, 217713.384600405, 709134.127242793,
+                                639828.157159169, 7786.82949717788, 1049204.06569028,
+                                1162623.7238134]
+
+        else:
+            spheroid_distances = [60504.0628825298, 77023.948962654, 49154.8867507115,
+                                  90847.435881812, 217402.811862568, 709599.234619957,
+                                  640011.483583758, 7772.00667666425, 1047861.7859506,
+                                  1165126.55237647]
+            sphere_distances = [60580.7612632291, 77143.7785056615, 49199.2725132184,
+                                90804.4414289463, 217712.63666124, 709131.691061906,
+                                639825.959074112, 7786.80274606706, 1049200.46122281,
+                                1162619.7297006]
+
+        # Testing with spheroid distances first.
+        hillsdale = AustraliaCity.objects.get(name='Hillsdale')
+        qs = AustraliaCity.objects.exclude(id=hillsdale.id).annotate(
+            distance=Distance('point', hillsdale.point, spheroid=True)
+        ).order_by('id')
+        for i, c in enumerate(qs):
+            self.assertAlmostEqual(spheroid_distances[i], c.distance.m, tol)
+        if postgis:
+            # PostGIS uses sphere-only distances by default, testing these as well.
+            qs = AustraliaCity.objects.exclude(id=hillsdale.id).annotate(
+                distance=Distance('point', hillsdale.point)
+            ).order_by('id')
+            for i, c in enumerate(qs):
+                self.assertAlmostEqual(sphere_distances[i], c.distance.m, tol)
+
+    @no_oracle  # Oracle already handles geographic distance calculation.
+    @skipUnlessDBFeature("has_Distance_function", 'has_Transform_function')
+    def test_distance_transform(self):
+        """
+        Test the `Distance` function used with `Transform` on a geographic field.
+        """
+        # We'll be using a Polygon (created by buffering the centroid
+        # of 77005 to 100m) -- which aren't allowed in geographic distance
+        # queries normally, however our field has been transformed to
+        # a non-geographic system.
+        z = SouthTexasZipcode.objects.get(name='77005')
+
+        # Reference query:
+        # SELECT ST_Distance(ST_Transform("distapp_censuszipcode"."poly", 32140),
+        #   ST_GeomFromText('<buffer_wkt>', 32140))
+        # FROM "distapp_censuszipcode";
+        dists_m = [3553.30384972258, 1243.18391525602, 2186.15439472242]
+
+        # Having our buffer in the SRID of the transformation and of the field
+        # -- should get the same results. The first buffer has no need for
+        # transformation SQL because it is the same SRID as what was given
+        # to `transform()`.  The second buffer will need to be transformed,
+        # however.
+        buf1 = z.poly.centroid.buffer(100)
+        buf2 = buf1.transform(4269, clone=True)
+        ref_zips = ['77002', '77025', '77401']
+
+        for buf in [buf1, buf2]:
+            qs = CensusZipcode.objects.exclude(name='77005').annotate(
+                distance=Distance(Transform('poly', 32140), buf)
+            ).order_by('name')
+            self.assertEqual(ref_zips, sorted([c.name for c in qs]))
+            for i, z in enumerate(qs):
+                self.assertAlmostEqual(z.distance.m, dists_m[i], 5)
+
+    @skipUnlessDBFeature("has_Distance_function")
+    def test_distance_order_by(self):
+        qs = SouthTexasCity.objects.annotate(distance=Distance('point', Point(3, 3, srid=32140))).order_by(
+            'distance'
+        ).values_list('name', flat=True).filter(name__in=('San Antonio', 'Pearland'))
+        self.assertQuerysetEqual(qs, ['San Antonio', 'Pearland'], lambda x: x)
+
+    @skipUnlessDBFeature("has_Length_function")
+    def test_length(self):
+        """
+        Test the `Length` function.
+        """
+        # Reference query (should use `length_spheroid`).
+        # SELECT ST_length_spheroid(ST_GeomFromText('<wkt>', 4326) 'SPHEROID["WGS 84",6378137,298.257223563,
+        #   AUTHORITY["EPSG","7030"]]');
+        len_m1 = 473504.769553813
+        len_m2 = 4617.668
+
+        if connection.features.supports_distance_geodetic:
+            qs = Interstate.objects.annotate(length=Length('path'))
+            tol = 2 if oracle else 3
+            self.assertAlmostEqual(len_m1, qs[0].length.m, tol)
+        else:
+            # Does not support geodetic coordinate systems.
+            self.assertRaises(ValueError, Interstate.objects.annotate(length=Length('path')))
+
+        # Now doing length on a projected coordinate system.
+        i10 = SouthTexasInterstate.objects.annotate(length=Length('path')).get(name='I-10')
+        self.assertAlmostEqual(len_m2, i10.length.m, 2)
+        self.assertTrue(
+            SouthTexasInterstate.objects.annotate(length=Length('path')).filter(length__gt=4000).exists()
+        )
+
+    @skipUnlessDBFeature("has_Perimeter_function")
+    def test_perimeter(self):
+        """
+        Test the `Perimeter` function.
+        """
+        # Reference query:
+        # SELECT ST_Perimeter(distapp_southtexaszipcode.poly) FROM distapp_southtexaszipcode;
+        perim_m = [18404.3550889361, 15627.2108551001, 20632.5588368978, 17094.5996143697]
+        tol = 2 if oracle else 7
+        qs = SouthTexasZipcode.objects.annotate(perimeter=Perimeter('poly')).order_by('name')
+        for i, z in enumerate(qs):
+            self.assertAlmostEqual(perim_m[i], z.perimeter.m, tol)
+
+        # Running on points; should return 0.
+        qs = SouthTexasCity.objects.annotate(perim=Perimeter('point'))
+        for city in qs:
+            self.assertEqual(0, city.perim.m)
+
+    @skipUnlessDBFeature("has_Area_function", "has_Distance_function")
+    def test_measurement_null_fields(self):
+        """
+        Test the measurement functions on fields with NULL values.
+        """
+        # Creating SouthTexasZipcode w/NULL value.
+        SouthTexasZipcode.objects.create(name='78212')
+        # Performing distance/area queries against the NULL PolygonField,
+        # and ensuring the result of the operations is None.
+        htown = SouthTexasCity.objects.get(name='Downtown Houston')
+        z = SouthTexasZipcode.objects.annotate(
+            distance=Distance('poly', htown.point), area=Area('poly')
+        ).get(name='78212')
+        self.assertIsNone(z.distance)
+        self.assertIsNone(z.area)

tests/gis_tests/geo3d/tests.py

@@ -4,6 +4,9 @@ import os
 import re
 from unittest import skipUnless
 
+from django.contrib.gis.db.models.functions import (
+    AsGeoJSON, AsKML, Length, Perimeter, Scale, Translate,
+)
 from django.contrib.gis.gdal import HAS_GDAL
 from django.contrib.gis.geos import HAS_GEOS
 from django.test import TestCase, ignore_warnings, skipUnlessDBFeature
@@ -73,18 +76,7 @@ bbox_data = (
 )
 
 
-@skipUnless(HAS_GDAL, "GDAL is required for Geo3DTest.")
+class Geo3DLoadingHelper(object):
-@skipUnlessDBFeature("gis_enabled", "supports_3d_storage")
-class Geo3DTest(TestCase):
-    """
-    Only a subset of the PostGIS routines are 3D-enabled, and this TestCase
-    tries to test the features that can handle 3D and that are also
-    available within GeoDjango.  For more information, see the PostGIS docs
-    on the routines that support 3D:
-
-    http://postgis.net/docs/PostGIS_Special_Functions_Index.html#PostGIS_3D_Functions
-    """
-
     def _load_interstate_data(self):
         # Interstate (2D / 3D and Geographic/Projected variants)
         for name, line, exp_z in interstate_data:
@@ -109,6 +101,19 @@ class Geo3DTest(TestCase):
         Polygon2D.objects.create(name='2D BBox', poly=bbox_2d)
         Polygon3D.objects.create(name='3D BBox', poly=bbox_3d)
 
+
+@skipUnless(HAS_GDAL, "GDAL is required for Geo3DTest.")
+@skipUnlessDBFeature("gis_enabled", "supports_3d_storage")
+class Geo3DTest(Geo3DLoadingHelper, TestCase):
+    """
+    Only a subset of the PostGIS routines are 3D-enabled, and this TestCase
+    tries to test the features that can handle 3D and that are also
+    available within GeoDjango.  For more information, see the PostGIS docs
+    on the routines that support 3D:
+
+    http://postgis.net/docs/PostGIS_Special_Functions_Index.html#PostGIS_3D_Functions
+    """
+
     def test_3d_hasz(self):
         """
         Make sure data is 3D and has expected Z values -- shouldn't change
@@ -302,3 +307,93 @@ class Geo3DTest(TestCase):
         for ztrans in ztranslations:
             for city in City3D.objects.translate(0, 0, ztrans):
                 self.assertEqual(city_dict[city.name][2] + ztrans, city.translate.z)
+
+
+@skipUnless(HAS_GDAL, "GDAL is required for Geo3DTest.")
+@skipUnlessDBFeature("gis_enabled", "supports_3d_functions")
+class Geo3DFunctionsTests(Geo3DLoadingHelper, TestCase):
+    def test_kml(self):
+        """
+        Test KML() function with Z values.
+        """
+        self._load_city_data()
+        h = City3D.objects.annotate(kml=AsKML('point', precision=6)).get(name='Houston')
+        # KML should be 3D.
+        # `SELECT ST_AsKML(point, 6) FROM geo3d_city3d WHERE name = 'Houston';`
+        ref_kml_regex = re.compile(r'^<Point><coordinates>-95.363\d+,29.763\d+,18</coordinates></Point>$')
+        self.assertTrue(ref_kml_regex.match(h.kml))
+
+    def test_geojson(self):
+        """
+        Test GeoJSON() function with Z values.
+        """
+        self._load_city_data()
+        h = City3D.objects.annotate(geojson=AsGeoJSON('point', precision=6)).get(name='Houston')
+        # GeoJSON should be 3D
+        # `SELECT ST_AsGeoJSON(point, 6) FROM geo3d_city3d WHERE name='Houston';`
+        ref_json_regex = re.compile(r'^{"type":"Point","coordinates":\[-95.363151,29.763374,18(\.0+)?\]}$')
+        self.assertTrue(ref_json_regex.match(h.geojson))
+
+    def test_perimeter(self):
+        """
+        Testing Perimeter() function on 3D fields.
+        """
+        self._load_polygon_data()
+        # Reference query for values below:
+        #  `SELECT ST_Perimeter3D(poly), ST_Perimeter2D(poly) FROM geo3d_polygon3d;`
+        ref_perim_3d = 76859.2620451
+        ref_perim_2d = 76859.2577803
+        tol = 6
+        poly2d = Polygon2D.objects.annotate(perimeter=Perimeter('poly')).get(name='2D BBox')
+        self.assertAlmostEqual(ref_perim_2d, poly2d.perimeter.m, tol)
+        poly3d = Polygon3D.objects.annotate(perimeter=Perimeter('poly')).get(name='3D BBox')
+        self.assertAlmostEqual(ref_perim_3d, poly3d.perimeter.m, tol)
+
+    def test_length(self):
+        """
+        Testing Length() function on 3D fields.
+        """
+        # ST_Length_Spheroid Z-aware, and thus does not need to use
+        # a separate function internally.
+        # `SELECT ST_Length_Spheroid(line, 'SPHEROID["GRS 1980",6378137,298.257222101]')
+        #    FROM geo3d_interstate[2d|3d];`
+        self._load_interstate_data()
+        tol = 3
+        ref_length_2d = 4368.1721949481
+        ref_length_3d = 4368.62547052088
+        inter2d = Interstate2D.objects.annotate(length=Length('line')).get(name='I-45')
+        self.assertAlmostEqual(ref_length_2d, inter2d.length.m, tol)
+        inter3d = Interstate3D.objects.annotate(length=Length('line')).get(name='I-45')
+        self.assertAlmostEqual(ref_length_3d, inter3d.length.m, tol)
+
+        # Making sure `ST_Length3D` is used on for a projected
+        # and 3D model rather than `ST_Length`.
+        # `SELECT ST_Length(line) FROM geo3d_interstateproj2d;`
+        ref_length_2d = 4367.71564892392
+        # `SELECT ST_Length3D(line) FROM geo3d_interstateproj3d;`
+        ref_length_3d = 4368.16897234101
+        inter2d = InterstateProj2D.objects.annotate(length=Length('line')).get(name='I-45')
+        self.assertAlmostEqual(ref_length_2d, inter2d.length.m, tol)
+        inter3d = InterstateProj3D.objects.annotate(length=Length('line')).get(name='I-45')
+        self.assertAlmostEqual(ref_length_3d, inter3d.length.m, tol)
+
+    def test_scale(self):
+        """
+        Testing Scale() function on Z values.
+        """
+        self._load_city_data()
+        # Mapping of City name to reference Z values.
+        zscales = (-3, 4, 23)
+        for zscale in zscales:
+            for city in City3D.objects.annotate(scale=Scale('point', 1.0, 1.0, zscale)):
+                self.assertEqual(city_dict[city.name][2] * zscale, city.scale.z)
+
+    def test_translate(self):
+        """
+        Testing Translate() function on Z values.
+        """
+        self._load_city_data()
+        ztranslations = (5.23, 23, -17)
+        for ztrans in ztranslations:
+            for city in City3D.objects.annotate(translate=Translate('point', 0, 0, ztrans)):
+                self.assertEqual(city_dict[city.name][2] + ztrans, city.translate.z)

tests/gis_tests/geoapp/test_functions.py

@@ -0,0 +1,447 @@
+from __future__ import unicode_literals
+
+import re
+from decimal import Decimal
+
+from django.contrib.gis.db.models import functions
+from django.contrib.gis.geos import HAS_GEOS
+from django.db import connection
+from django.test import TestCase, skipUnlessDBFeature
+from django.utils import six
+
+from ..utils import oracle, postgis, spatialite
+
+if HAS_GEOS:
+    from django.contrib.gis.geos import LineString, Point, Polygon, fromstr
+    from .models import Country, City, State, Track
+
+
+@skipUnlessDBFeature("gis_enabled")
+class GISFunctionsTests(TestCase):
+    """
+    Testing functions from django/contrib/gis/db/models/functions.py.
+    Several tests are taken and adapted from GeoQuerySetTest.
+    Area/Distance/Length/Perimeter are tested in distapp/tests.
+
+    Please keep the tests in function's alphabetic order.
+    """
+    fixtures = ['initial']
+
+    def test_asgeojson(self):
+        # Only PostGIS and SpatiaLite 3.0+ support GeoJSON.
+        if not connection.ops.geojson:
+            with self.assertRaises(NotImplementedError):
+                list(Country.objects.annotate(json=functions.AsGeoJSON('mpoly')))
+            return
+
+        pueblo_json = '{"type":"Point","coordinates":[-104.609252,38.255001]}'
+        houston_json = (
+            '{"type":"Point","crs":{"type":"name","properties":'
+            '{"name":"EPSG:4326"}},"coordinates":[-95.363151,29.763374]}'
+        )
+        victoria_json = (
+            '{"type":"Point","bbox":[-123.30519600,48.46261100,-123.30519600,48.46261100],'
+            '"coordinates":[-123.305196,48.462611]}'
+        )
+        chicago_json = (
+            '{"type":"Point","crs":{"type":"name","properties":{"name":"EPSG:4326"}},'
+            '"bbox":[-87.65018,41.85039,-87.65018,41.85039],"coordinates":[-87.65018,41.85039]}'
+        )
+        if spatialite:
+            victoria_json = (
+                '{"type":"Point","bbox":[-123.305196,48.462611,-123.305196,48.462611],'
+                '"coordinates":[-123.305196,48.462611]}'
+            )
+
+        # Precision argument should only be an integer
+        with self.assertRaises(TypeError):
+            City.objects.annotate(geojson=functions.AsGeoJSON('point', precision='foo'))
+
+        # Reference queries and values.
+        # SELECT ST_AsGeoJson("geoapp_city"."point", 8, 0)
+        # FROM "geoapp_city" WHERE "geoapp_city"."name" = 'Pueblo';
+        self.assertEqual(
+            pueblo_json,
+            City.objects.annotate(geojson=functions.AsGeoJSON('point')).get(name='Pueblo').geojson
+        )
+
+        # SELECT ST_AsGeoJson("geoapp_city"."point", 8, 2) FROM "geoapp_city"
+        # WHERE "geoapp_city"."name" = 'Houston';
+        # This time we want to include the CRS by using the `crs` keyword.
+        self.assertEqual(
+            houston_json,
+            City.objects.annotate(json=functions.AsGeoJSON('point', crs=True)).get(name='Houston').json
+        )
+
+        # SELECT ST_AsGeoJson("geoapp_city"."point", 8, 1) FROM "geoapp_city"
+        # WHERE "geoapp_city"."name" = 'Houston';
+        # This time we include the bounding box by using the `bbox` keyword.
+        self.assertEqual(
+            victoria_json,
+            City.objects.annotate(
+                geojson=functions.AsGeoJSON('point', bbox=True)
+            ).get(name='Victoria').geojson
+        )
+
+        # SELECT ST_AsGeoJson("geoapp_city"."point", 5, 3) FROM "geoapp_city"
+        # WHERE "geoapp_city"."name" = 'Chicago';
+        # Finally, we set every available keyword.
+        self.assertEqual(
+            chicago_json,
+            City.objects.annotate(
+                geojson=functions.AsGeoJSON('point', bbox=True, crs=True, precision=5)
+                ).get(name='Chicago').geojson
+        )
+
+    @skipUnlessDBFeature("has_AsGML_function")
+    def test_asgml(self):
+        # Should throw a TypeError when tyring to obtain GML from a
+        # non-geometry field.
+        qs = City.objects.all()
+        with self.assertRaises(TypeError):
+            qs.annotate(gml=functions.AsGML('name'))
+        ptown = City.objects.annotate(gml=functions.AsGML('point', precision=9)).get(name='Pueblo')
+
+        if oracle:
+            # No precision parameter for Oracle :-/
+            gml_regex = re.compile(
+                r'^<gml:Point srsName="SDO:4326" xmlns:gml="http://www.opengis.net/gml">'
+                r'<gml:coordinates decimal="\." cs="," ts=" ">-104.60925\d+,38.25500\d+ '
+                r'</gml:coordinates></gml:Point>'
+            )
+        elif spatialite and connection.ops.spatial_version < (3, 0, 0):
+            # Spatialite before 3.0 has extra colon in SrsName
+            gml_regex = re.compile(
+                r'^<gml:Point SrsName="EPSG::4326"><gml:coordinates decimal="\." '
+                r'cs="," ts=" ">-104.609251\d+,38.255001</gml:coordinates></gml:Point>'
+            )
+        else:
+            gml_regex = re.compile(
+                r'^<gml:Point srsName="EPSG:4326"><gml:coordinates>'
+                r'-104\.60925\d+,38\.255001</gml:coordinates></gml:Point>'
+            )
+
+        self.assertTrue(gml_regex.match(ptown.gml))
+
+        if postgis:
+            self.assertIn(
+                '<gml:pos srsDimension="2">',
+                City.objects.annotate(gml=functions.AsGML('point', version=3)).get(name='Pueblo').gml
+            )
+
+    @skipUnlessDBFeature("has_AsKML_function")
+    def test_askml(self):
+        # Should throw a TypeError when trying to obtain KML from a
+        # non-geometry field.
+        with self.assertRaises(TypeError):
+            City.objects.annotate(kml=functions.AsKML('name'))
+
+        # Ensuring the KML is as expected.
+        ptown = City.objects.annotate(kml=functions.AsKML('point', precision=9)).get(name='Pueblo')
+        self.assertEqual('<Point><coordinates>-104.609252,38.255001</coordinates></Point>', ptown.kml)
+
+    @skipUnlessDBFeature("has_AsSVG_function")
+    def test_assvg(self):
+        with self.assertRaises(TypeError):
+            City.objects.annotate(svg=functions.AsSVG('point', precision='foo'))
+        # SELECT AsSVG(geoapp_city.point, 0, 8) FROM geoapp_city WHERE name = 'Pueblo';
+        svg1 = 'cx="-104.609252" cy="-38.255001"'
+        # Even though relative, only one point so it's practically the same except for
+        # the 'c' letter prefix on the x,y values.
+        svg2 = svg1.replace('c', '')
+        self.assertEqual(svg1, City.objects.annotate(svg=functions.AsSVG('point')).get(name='Pueblo').svg)
+        self.assertEqual(svg2, City.objects.annotate(svg=functions.AsSVG('point', relative=5)).get(name='Pueblo').svg)
+
+    @skipUnlessDBFeature("has_BoundingCircle_function")
+    def test_bounding_circle(self):
+        qs = Country.objects.annotate(circle=functions.BoundingCircle('mpoly')).order_by('name')
+        self.assertAlmostEqual(qs[0].circle.area, 168.89, 2)
+        self.assertAlmostEqual(qs[1].circle.area, 135.95, 2)
+
+        qs = Country.objects.annotate(circle=functions.BoundingCircle('mpoly', num_seg=12)).order_by('name')
+        self.assertAlmostEqual(qs[0].circle.area, 168.44, 2)
+        self.assertAlmostEqual(qs[1].circle.area, 135.59, 2)
+
+    @skipUnlessDBFeature("has_Centroid_function")
+    def test_centroid(self):
+        qs = State.objects.exclude(poly__isnull=True).annotate(centroid=functions.Centroid('poly'))
+        for state in qs:
+            tol = 0.1  # High tolerance due to oracle
+            self.assertTrue(state.poly.centroid.equals_exact(state.centroid, tol))
+
+    @skipUnlessDBFeature("has_Difference_function")
+    def test_difference(self):
+        geom = Point(5, 23, srid=4326)
+        qs = Country.objects.annotate(difference=functions.Difference('mpoly', geom))
+        for c in qs:
+            self.assertEqual(c.mpoly.difference(geom), c.difference)
+
+    @skipUnlessDBFeature("has_Difference_function")
+    def test_difference_mixed_srid(self):
+        """Testing with mixed SRID (Country has default 4326)."""
+        geom = Point(556597.4, 2632018.6, srid=3857)  # Spherical mercator
+        qs = Country.objects.annotate(difference=functions.Difference('mpoly', geom))
+        for c in qs:
+            self.assertEqual(c.mpoly.difference(geom), c.difference)
+
+    @skipUnlessDBFeature("has_Envelope_function")
+    def test_envelope(self):
+        countries = Country.objects.annotate(envelope=functions.Envelope('mpoly'))
+        for country in countries:
+            self.assertIsInstance(country.envelope, Polygon)
+
+    @skipUnlessDBFeature("has_ForceRHR_function")
+    def test_force_rhr(self):
+        rings = (
+            ((0, 0), (5, 0), (0, 5), (0, 0)),
+            ((1, 1), (1, 3), (3, 1), (1, 1)),
+        )
+        rhr_rings = (
+            ((0, 0), (0, 5), (5, 0), (0, 0)),
+            ((1, 1), (3, 1), (1, 3), (1, 1)),
+        )
+        State.objects.create(name='Foo', poly=Polygon(*rings))
+        st = State.objects.annotate(force_rhr=functions.ForceRHR('poly')).get(name='Foo')
+        self.assertEqual(rhr_rings, st.force_rhr.coords)
+
+    @skipUnlessDBFeature("has_GeoHash_function")
+    def test_geohash(self):
+        # Reference query:
+        # SELECT ST_GeoHash(point) FROM geoapp_city WHERE name='Houston';
+        # SELECT ST_GeoHash(point, 5) FROM geoapp_city WHERE name='Houston';
+        ref_hash = '9vk1mfq8jx0c8e0386z6'
+        h1 = City.objects.annotate(geohash=functions.GeoHash('point')).get(name='Houston')
+        h2 = City.objects.annotate(geohash=functions.GeoHash('point', precision=5)).get(name='Houston')
+        self.assertEqual(ref_hash, h1.geohash)
+        self.assertEqual(ref_hash[:5], h2.geohash)
+
+    @skipUnlessDBFeature("has_Intersection_function")
+    def test_intersection(self):
+        geom = Point(5, 23, srid=4326)
+        qs = Country.objects.annotate(inter=functions.Intersection('mpoly', geom))
+        for c in qs:
+            self.assertEqual(c.mpoly.intersection(geom), c.inter)
+
+    @skipUnlessDBFeature("has_MemSize_function")
+    def test_memsize(self):
+        ptown = City.objects.annotate(size=functions.MemSize('point')).get(name='Pueblo')
+        self.assertTrue(20 <= ptown.size <= 40)  # Exact value may depend on PostGIS version
+
+    @skipUnlessDBFeature("has_NumGeom_function")
+    def test_num_geom(self):
+        # Both 'countries' only have two geometries.
+        for c in Country.objects.annotate(num_geom=functions.NumGeometries('mpoly')):
+            self.assertEqual(2, c.num_geom)
+
+        qs = City.objects.filter(point__isnull=False).annotate(num_geom=functions.NumGeometries('point'))
+        for city in qs:
+            # Oracle and PostGIS 2.0+ will return 1 for the number of
+            # geometries on non-collections, whereas PostGIS < 2.0.0
+            # will return None.
+            if postgis and connection.ops.spatial_version < (2, 0, 0):
+                self.assertIsNone(city.num_geom)
+            else:
+                self.assertEqual(1, city.num_geom)
+
+    @skipUnlessDBFeature("has_NumPoint_function")
+    def test_num_points(self):
+        coords = [(-95.363151, 29.763374), (-95.448601, 29.713803)]
+        Track.objects.create(name='Foo', line=LineString(coords))
+        qs = Track.objects.annotate(num_points=functions.NumPoints('line'))
+        self.assertEqual(qs.first().num_points, 2)
+        if spatialite:
+            # Spatialite can only count points on LineStrings
+            return
+
+        for c in Country.objects.annotate(num_points=functions.NumPoints('mpoly')):
+            self.assertEqual(c.mpoly.num_points, c.num_points)
+
+        if not oracle:
+            # Oracle cannot count vertices in Point geometries.
+            for c in City.objects.annotate(num_points=functions.NumPoints('point')):
+                self.assertEqual(1, c.num_points)
+
+    @skipUnlessDBFeature("has_PointOnSurface_function")
+    def test_point_on_surface(self):
+        # Reference values.
+        if oracle:
+            # SELECT SDO_UTIL.TO_WKTGEOMETRY(SDO_GEOM.SDO_POINTONSURFACE(GEOAPP_COUNTRY.MPOLY, 0.05))
+            # FROM GEOAPP_COUNTRY;
+            ref = {'New Zealand': fromstr('POINT (174.616364 -36.100861)', srid=4326),
+                   'Texas': fromstr('POINT (-103.002434 36.500397)', srid=4326),
+                   }
+        else:
+            # Using GEOSGeometry to compute the reference point on surface values
+            # -- since PostGIS also uses GEOS these should be the same.
+            ref = {'New Zealand': Country.objects.get(name='New Zealand').mpoly.point_on_surface,
+                   'Texas': Country.objects.get(name='Texas').mpoly.point_on_surface
+                   }
+
+        qs = Country.objects.annotate(point_on_surface=functions.PointOnSurface('mpoly'))
+        for country in qs:
+            tol = 0.00001  # Spatialite might have WKT-translation-related precision issues
+            self.assertTrue(ref[country.name].equals_exact(country.point_on_surface, tol))
+
+    @skipUnlessDBFeature("has_Reverse_function")
+    def test_reverse_geom(self):
+        coords = [(-95.363151, 29.763374), (-95.448601, 29.713803)]
+        Track.objects.create(name='Foo', line=LineString(coords))
+        track = Track.objects.annotate(reverse_geom=functions.Reverse('line')).get(name='Foo')
+        coords.reverse()
+        self.assertEqual(tuple(coords), track.reverse_geom.coords)
+
+    @skipUnlessDBFeature("has_Scale_function")
+    def test_scale(self):
+        xfac, yfac = 2, 3
+        tol = 5  # The low precision tolerance is for SpatiaLite
+        qs = Country.objects.annotate(scaled=functions.Scale('mpoly', xfac, yfac))
+        for country in qs:
+            for p1, p2 in zip(country.mpoly, country.scaled):
+                for r1, r2 in zip(p1, p2):
+                    for c1, c2 in zip(r1.coords, r2.coords):
+                        self.assertAlmostEqual(c1[0] * xfac, c2[0], tol)
+                        self.assertAlmostEqual(c1[1] * yfac, c2[1], tol)
+        # Test float/Decimal values
+        qs = Country.objects.annotate(scaled=functions.Scale('mpoly', 1.5, Decimal('2.5')))
+        self.assertGreater(qs[0].scaled.area, qs[0].mpoly.area)
+
+    @skipUnlessDBFeature("has_SnapToGrid_function")
+    def test_snap_to_grid(self):
+        # Let's try and break snap_to_grid() with bad combinations of arguments.
+        for bad_args in ((), range(3), range(5)):
+            with self.assertRaises(ValueError):
+                Country.objects.annotate(snap=functions.SnapToGrid('mpoly', *bad_args))
+        for bad_args in (('1.0',), (1.0, None), tuple(map(six.text_type, range(4)))):
+            with self.assertRaises(TypeError):
+                Country.objects.annotate(snap=functions.SnapToGrid('mpoly', *bad_args))
+
+        # Boundary for San Marino, courtesy of Bjorn Sandvik of thematicmapping.org
+        # from the world borders dataset he provides.
+        wkt = ('MULTIPOLYGON(((12.41580 43.95795,12.45055 43.97972,12.45389 43.98167,'
+               '12.46250 43.98472,12.47167 43.98694,12.49278 43.98917,'
+               '12.50555 43.98861,12.51000 43.98694,12.51028 43.98277,'
+               '12.51167 43.94333,12.51056 43.93916,12.49639 43.92333,'
+               '12.49500 43.91472,12.48778 43.90583,12.47444 43.89722,'
+               '12.46472 43.89555,12.45917 43.89611,12.41639 43.90472,'
+               '12.41222 43.90610,12.40782 43.91366,12.40389 43.92667,'
+               '12.40500 43.94833,12.40889 43.95499,12.41580 43.95795)))')
+        Country.objects.create(name='San Marino', mpoly=fromstr(wkt))
+
+        # Because floating-point arithmetic isn't exact, we set a tolerance
+        # to pass into GEOS `equals_exact`.
+        tol = 0.000000001
+
+        # SELECT AsText(ST_SnapToGrid("geoapp_country"."mpoly", 0.1)) FROM "geoapp_country"
+        # WHERE "geoapp_country"."name" = 'San Marino';
+        ref = fromstr('MULTIPOLYGON(((12.4 44,12.5 44,12.5 43.9,12.4 43.9,12.4 44)))')
+        self.assertTrue(
+            ref.equals_exact(
+                Country.objects.annotate(
+                    snap=functions.SnapToGrid('mpoly', 0.1)
+                ).get(name='San Marino').snap,
+                tol
+            )
+        )
+
+        # SELECT AsText(ST_SnapToGrid("geoapp_country"."mpoly", 0.05, 0.23)) FROM "geoapp_country"
+        # WHERE "geoapp_country"."name" = 'San Marino';
+        ref = fromstr('MULTIPOLYGON(((12.4 43.93,12.45 43.93,12.5 43.93,12.45 43.93,12.4 43.93)))')
+        self.assertTrue(
+            ref.equals_exact(
+                Country.objects.annotate(
+                    snap=functions.SnapToGrid('mpoly', 0.05, 0.23)
+                ).get(name='San Marino').snap,
+                tol
+            )
+        )
+
+        # SELECT AsText(ST_SnapToGrid("geoapp_country"."mpoly", 0.5, 0.17, 0.05, 0.23)) FROM "geoapp_country"
+        # WHERE "geoapp_country"."name" = 'San Marino';
+        ref = fromstr(
+            'MULTIPOLYGON(((12.4 43.87,12.45 43.87,12.45 44.1,12.5 44.1,12.5 43.87,12.45 43.87,12.4 43.87)))'
+        )
+        self.assertTrue(
+            ref.equals_exact(
+                Country.objects.annotate(
+                    snap=functions.SnapToGrid('mpoly', 0.05, 0.23, 0.5, 0.17)
+                ).get(name='San Marino').snap,
+                tol
+            )
+        )
+
+    @skipUnlessDBFeature("has_SymDifference_function")
+    def test_sym_difference(self):
+        geom = Point(5, 23, srid=4326)
+        qs = Country.objects.annotate(sym_difference=functions.SymDifference('mpoly', geom))
+        for country in qs:
+            # Ordering might differ in collections
+            self.assertSetEqual(set(g.wkt for g in country.mpoly.sym_difference(geom)),
+                                set(g.wkt for g in country.sym_difference))
+
+    @skipUnlessDBFeature("has_Transform_function")
+    def test_transform(self):
+        # Pre-transformed points for Houston and Pueblo.
+        ptown = fromstr('POINT(992363.390841912 481455.395105533)', srid=2774)
+        prec = 3  # Precision is low due to version variations in PROJ and GDAL.
+
+        # Asserting the result of the transform operation with the values in
+        #  the pre-transformed points.
+        h = City.objects.annotate(pt=functions.Transform('point', ptown.srid)).get(name='Pueblo')
+        self.assertEqual(2774, h.pt.srid)
+        self.assertAlmostEqual(ptown.x, h.pt.x, prec)
+        self.assertAlmostEqual(ptown.y, h.pt.y, prec)
+
+    @skipUnlessDBFeature("has_Translate_function")
+    def test_translate(self):
+        xfac, yfac = 5, -23
+        qs = Country.objects.annotate(translated=functions.Translate('mpoly', xfac, yfac))
+        for c in qs:
+            for p1, p2 in zip(c.mpoly, c.translated):
+                for r1, r2 in zip(p1, p2):
+                    for c1, c2 in zip(r1.coords, r2.coords):
+                        # The low precision is for SpatiaLite
+                        self.assertAlmostEqual(c1[0] + xfac, c2[0], 5)
+                        self.assertAlmostEqual(c1[1] + yfac, c2[1], 5)
+
+    # Some combined function tests
+    @skipUnlessDBFeature(
+        "has_Difference_function", "has_Intersection_function",
+        "has_SymDifference_function", "has_Union_function")
+    def test_diff_intersection_union(self):
+        "Testing the `difference`, `intersection`, `sym_difference`, and `union` GeoQuerySet methods."
+        geom = Point(5, 23, srid=4326)
+        qs = Country.objects.all().annotate(
+            difference=functions.Difference('mpoly', geom),
+            sym_difference=functions.SymDifference('mpoly', geom),
+            union=functions.Union('mpoly', geom),
+        )
+
+        # XXX For some reason SpatiaLite does something screwey with the Texas geometry here.  Also,
+        # XXX it doesn't like the null intersection.
+        if spatialite:
+            qs = qs.exclude(name='Texas')
+        else:
+            qs = qs.annotate(intersection=functions.Intersection('mpoly', geom))
+
+        if oracle:
+            # Should be able to execute the queries; however, they won't be the same
+            # as GEOS (because Oracle doesn't use GEOS internally like PostGIS or
+            # SpatiaLite).
+            return
+        for c in qs:
+            self.assertEqual(c.mpoly.difference(geom), c.difference)
+            if not spatialite:
+                self.assertEqual(c.mpoly.intersection(geom), c.intersection)
+            # Ordering might differ in collections
+            self.assertSetEqual(set(g.wkt for g in c.mpoly.sym_difference(geom)),
+                                set(g.wkt for g in c.sym_difference))
+            self.assertSetEqual(set(g.wkt for g in c.mpoly.union(geom)),
+                                set(g.wkt for g in c.union))
+
+    @skipUnlessDBFeature("has_Union_function")
+    def test_union(self):
+        geom = Point(-95.363151, 29.763374, srid=4326)
+        ptown = City.objects.annotate(union=functions.Union('point', geom)).get(name='Dallas')
+        tol = 0.00001
+        expected = fromstr('MULTIPOINT(-96.801611 32.782057,-95.363151 29.763374)', srid=4326)
+        self.assertTrue(expected.equals_exact(ptown.union, tol))

tests/gis_tests/geogapp/tests.py

@@ -6,6 +6,7 @@ from __future__ import unicode_literals
 import os
 from unittest import skipUnless
 
+from django.contrib.gis.db.models.functions import Area, Distance
 from django.contrib.gis.gdal import HAS_GDAL
 from django.contrib.gis.geos import HAS_GEOS
 from django.contrib.gis.measure import D
@@ -101,3 +102,30 @@ class GeographyTest(TestCase):
         tol = 5
         z = Zipcode.objects.area().get(code='77002')
         self.assertAlmostEqual(z.area.sq_m, ref_area, tol)
+
+
+@skipUnlessDBFeature("gis_enabled")
+class GeographyFunctionTests(TestCase):
+    fixtures = ['initial']
+
+    @skipUnlessDBFeature("has_Distance_function", "supports_distance_geodetic")
+    def test_distance_function(self):
+        """
+        Testing Distance() support on non-point geography fields.
+        """
+        ref_dists = [0, 4891.20, 8071.64, 9123.95]
+        htown = City.objects.get(name='Houston')
+        qs = Zipcode.objects.annotate(distance=Distance('poly', htown.point))
+        for z, ref in zip(qs, ref_dists):
+            self.assertAlmostEqual(z.distance.m, ref, 2)
+
+    @skipUnlessDBFeature("has_Area_function", "supports_distance_geodetic")
+    def test_geography_area(self):
+        """
+        Testing that Area calculations work on geography columns.
+        """
+        # SELECT ST_Area(poly) FROM geogapp_zipcode WHERE code='77002';
+        ref_area = 5439100.95415646 if oracle else 5439084.70637573
+        tol = 5
+        z = Zipcode.objects.annotate(area=Area('poly')).get(code='77002')
+        self.assertAlmostEqual(z.area.sq_m, ref_area, tol)