From 986c7d522a4c2351108ba1dc5f58b5fd4ed75ba0 Mon Sep 17 00:00:00 2001
From: Sergey Fedoseev <fedoseev.sergey@gmail.com>
Date: Wed, 16 Nov 2016 17:07:36 +0500
Subject: [PATCH] Fixed #27497 -- Improved support of geodetic coordinates on
 SpatiaLite.

Area function, Distance function, and distance lookups now work with
geodetic coordinates on SpatiaLite.
---
 .../contrib/gis/db/backends/base/features.py  |  1 +
 .../contrib/gis/db/backends/mysql/features.py |  1 +
 .../gis/db/backends/spatialite/features.py    |  5 ++-
 .../gis/db/backends/spatialite/operations.py  | 36 ++++++++++++------
 django/contrib/gis/db/models/functions.py     | 37 +++++++++++++++----
 django/contrib/gis/db/models/query.py         |  2 +-
 docs/releases/1.11.txt                        |  4 ++
 tests/gis_tests/distapp/tests.py              | 26 ++++++++-----
 tests/gis_tests/geogapp/tests.py              | 30 +++++++++++----
 9 files changed, 105 insertions(+), 37 deletions(-)

diff --git a/django/contrib/gis/db/backends/base/features.py b/django/contrib/gis/db/backends/base/features.py
index 3d1dfba1dd..53de204cb9 100644
--- a/django/contrib/gis/db/backends/base/features.py
+++ b/django/contrib/gis/db/backends/base/features.py
@@ -32,6 +32,7 @@ class BaseSpatialFeatures(object):
     supports_distance_geodetic = True
     supports_length_geodetic = True
     supports_perimeter_geodetic = False
+    supports_area_geodetic = True
     # Is the database able to count vertices on polygons (with `num_points`)?
     supports_num_points_poly = True
 
diff --git a/django/contrib/gis/db/backends/mysql/features.py b/django/contrib/gis/db/backends/mysql/features.py
index ce5893f0a8..05affeecb7 100644
--- a/django/contrib/gis/db/backends/mysql/features.py
+++ b/django/contrib/gis/db/backends/mysql/features.py
@@ -8,6 +8,7 @@ class DatabaseFeatures(BaseSpatialFeatures, MySQLDatabaseFeatures):
     supports_add_srs_entry = False
     supports_distance_geodetic = False
     supports_length_geodetic = False
+    supports_area_geodetic = False
     supports_distances_lookups = False
     supports_transform = False
     supports_real_shape_operations = False
diff --git a/django/contrib/gis/db/backends/spatialite/features.py b/django/contrib/gis/db/backends/spatialite/features.py
index 19195f1a57..79517e8190 100644
--- a/django/contrib/gis/db/backends/spatialite/features.py
+++ b/django/contrib/gis/db/backends/spatialite/features.py
@@ -6,7 +6,6 @@ from django.utils.functional import cached_property
 
 class DatabaseFeatures(BaseSpatialFeatures, SQLiteDatabaseFeatures):
     supports_3d_storage = True
-    supports_distance_geodetic = False
     # SpatiaLite can only count vertices in LineStrings
     supports_num_points_poly = False
 
@@ -16,3 +15,7 @@ class DatabaseFeatures(BaseSpatialFeatures, SQLiteDatabaseFeatures):
         # which can result in a significant performance improvement when
         # creating the database.
         return self.connection.ops.spatial_version >= (4, 1, 0)
+
+    @cached_property
+    def supports_area_geodetic(self):
+        return bool(self.connection.ops.lwgeom_version())
diff --git a/django/contrib/gis/db/backends/spatialite/operations.py b/django/contrib/gis/db/backends/spatialite/operations.py
index 73274a2fba..cc46b4b91f 100644
--- a/django/contrib/gis/db/backends/spatialite/operations.py
+++ b/django/contrib/gis/db/backends/spatialite/operations.py
@@ -19,6 +19,20 @@ from django.utils import six
 from django.utils.functional import cached_property
 
 
+class SpatiaLiteDistanceOperator(SpatialOperator):
+    def as_sql(self, connection, lookup, template_params, sql_params):
+        if lookup.lhs.output_field.geodetic(connection):
+            # SpatiaLite returns NULL instead of zero on geodetic coordinates
+            sql_template = 'COALESCE(%(func)s(%(lhs)s, %(rhs)s, %%s), 0) %(op)s %(value)s'
+            template_params.update({
+                'op': self.op,
+                'func': connection.ops.spatial_function_name('Distance'),
+            })
+            sql_params.insert(1, len(lookup.rhs) == 3 and lookup.rhs[-1] == 'spheroid')
+            return sql_template % template_params, sql_params
+        return super(SpatiaLiteDistanceOperator, self).as_sql(connection, lookup, template_params, sql_params)
+
+
 class SpatiaLiteOperations(BaseSpatialOperations, DatabaseOperations):
     name = 'spatialite'
     spatialite = True
@@ -79,10 +93,10 @@ class SpatiaLiteOperations(BaseSpatialOperations, DatabaseOperations):
         'exact': SpatialOperator(func='Equals'),
         # Distance predicates
         'dwithin': SpatialOperator(func='PtDistWithin'),
-        'distance_gt': SpatialOperator(func='Distance', op='>'),
-        'distance_gte': SpatialOperator(func='Distance', op='>='),
-        'distance_lt': SpatialOperator(func='Distance', op='<'),
-        'distance_lte': SpatialOperator(func='Distance', op='<='),
+        'distance_gt': SpatiaLiteDistanceOperator(func='Distance', op='>'),
+        'distance_gte': SpatiaLiteDistanceOperator(func='Distance', op='>='),
+        'distance_lt': SpatiaLiteDistanceOperator(func='Distance', op='<'),
+        'distance_lte': SpatiaLiteDistanceOperator(func='Distance', op='<='),
     }
 
     disallowed_aggregates = (aggregates.Extent3D,)
@@ -140,19 +154,19 @@ class SpatiaLiteOperations(BaseSpatialOperations, DatabaseOperations):
     def get_distance(self, f, value, lookup_type, **kwargs):
         """
         Returns the distance parameters for the given geometry field,
-        lookup value, and lookup type.  SpatiaLite only supports regular
-        cartesian-based queries (no spheroid/sphere calculations for point
-        geometries like PostGIS).
+        lookup value, and lookup type.
         """
         if not value:
             return []
         value = value[0]
         if isinstance(value, Distance):
             if f.geodetic(self.connection):
-                raise ValueError('SpatiaLite does not support distance queries on '
-                                 'geometry fields with a geodetic coordinate system. '
-                                 'Distance objects; use a numeric value of your '
-                                 'distance in degrees instead.')
+                if lookup_type == 'dwithin':
+                    raise ValueError(
+                        'Only numeric values of degree units are allowed on '
+                        'geographic DWithin queries.'
+                    )
+                dist_param = value.m
             else:
                 dist_param = getattr(value, Distance.unit_attname(f.units_name(self.connection)))
         else:
diff --git a/django/contrib/gis/db/models/functions.py b/django/contrib/gis/db/models/functions.py
index dea042493a..06b44e7b9f 100644
--- a/django/contrib/gis/db/models/functions.py
+++ b/django/contrib/gis/db/models/functions.py
@@ -38,6 +38,10 @@ class GeoFunc(Func):
         except (AttributeError, FieldError):
             return None
 
+    @property
+    def geo_field(self):
+        return GeometryField(srid=self.srid) if self.srid else None
+
     def as_sql(self, compiler, connection, **extra_context):
         if self.function is None:
             self.function = connection.ops.spatial_function_name(self.name)
@@ -122,26 +126,34 @@ class Area(OracleToleranceMixin, GeoFunc):
     output_field_class = AreaField
     arity = 1
 
-    def as_sql(self, compiler, connection):
+    def as_sql(self, compiler, connection, **extra_context):
         if connection.ops.geography:
             self.output_field.area_att = 'sq_m'
         else:
             # Getting the area units of the geographic field.
-            source_fields = self.get_source_fields()
-            if len(source_fields):
-                source_field = source_fields[0]
-                if source_field.geodetic(connection):
+            geo_field = self.geo_field
+            if geo_field.geodetic(connection):
+                if connection.features.supports_area_geodetic:
+                    self.output_field.area_att = 'sq_m'
+                else:
                     # TODO: Do we want to support raw number areas for geodetic fields?
                     raise NotImplementedError('Area on geodetic coordinate systems not supported.')
-                units_name = source_field.units_name(connection)
+            else:
+                units_name = geo_field.units_name(connection)
                 if units_name:
                     self.output_field.area_att = AreaMeasure.unit_attname(units_name)
-        return super(Area, self).as_sql(compiler, connection)
+        return super(Area, self).as_sql(compiler, connection, **extra_context)
 
     def as_oracle(self, compiler, connection):
         self.output_field = AreaField('sq_m')  # Oracle returns area in units of meters.
         return super(Area, self).as_oracle(compiler, connection)
 
+    def as_sqlite(self, compiler, connection, **extra_context):
+        if self.geo_field.geodetic(connection):
+            extra_context['template'] = '%(function)s(%(expressions)s, %(spheroid)d)'
+            extra_context['spheroid'] = True
+        return self.as_sql(compiler, connection, **extra_context)
+
 
 class AsGeoJSON(GeoFunc):
     output_field_class = TextField
@@ -226,7 +238,7 @@ 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
+        geo_field = self.geo_field
         if geo_field.geodetic(connection):
             dist_att = 'm'
         else:
@@ -275,6 +287,15 @@ class Distance(DistanceResultMixin, OracleToleranceMixin, GeoFuncWithGeoParam):
             self.source_expressions.pop(2)
         return super(Distance, self).as_oracle(compiler, connection)
 
+    def as_sqlite(self, compiler, connection, **extra_context):
+        if self.spheroid:
+            self.source_expressions.pop(2)
+        if self.geo_field.geodetic(connection):
+            # SpatiaLite returns NULL instead of zero on geodetic coordinates
+            extra_context['template'] = 'COALESCE(%(function)s(%(expressions)s, %(spheroid)s), 0)'
+            extra_context['spheroid'] = int(bool(self.spheroid))
+        return super(Distance, self).as_sql(compiler, connection, **extra_context)
+
 
 class Envelope(GeoFunc):
     arity = 1
diff --git a/django/contrib/gis/db/models/query.py b/django/contrib/gis/db/models/query.py
index d0a80f9e56..d1933721b6 100644
--- a/django/contrib/gis/db/models/query.py
+++ b/django/contrib/gis/db/models/query.py
@@ -546,7 +546,7 @@ class GeoQuerySet(QuerySet):
                 u, unit_name, s = get_srid_info(srid, connection)
                 geodetic = unit_name.lower() in geo_field.geodetic_units
 
-            if geodetic and not connection.features.supports_distance_geodetic:
+            if geodetic and (not connection.features.supports_distance_geodetic or connection.ops.spatialite):
                 raise ValueError(
                     'This database does not support linear distance '
                     'calculations on geodetic coordinate systems.'
diff --git a/docs/releases/1.11.txt b/docs/releases/1.11.txt
index 4eea7f0df8..9cfd7ae2d8 100644
--- a/docs/releases/1.11.txt
+++ b/docs/releases/1.11.txt
@@ -149,6 +149,10 @@ Minor features
 
 * Added support for the :lookup:`dwithin` lookup on SpatiaLite.
 
+* The :class:`~django.contrib.gis.db.models.functions.Area` function,
+  :class:`~django.contrib.gis.db.models.functions.Distance` function, and
+  distance lookups now work with geodetic coordinates on SpatiaLite.
+
 * The OpenLayers-based form widgets now use ``OpenLayers.js`` from
   ``https://cdnjs.cloudflare.com`` which is more suitable for production use
   than the the old ``http://openlayers.org`` source.
diff --git a/tests/gis_tests/distapp/tests.py b/tests/gis_tests/distapp/tests.py
index c943351952..a62ea46280 100644
--- a/tests/gis_tests/distapp/tests.py
+++ b/tests/gis_tests/distapp/tests.py
@@ -1,5 +1,7 @@
 from __future__ import unicode_literals
 
+from unittest import skipIf
+
 from django.contrib.gis.db.models.functions import (
     Area, Distance, Length, Perimeter, Transform,
 )
@@ -143,6 +145,7 @@ class DistanceTest(TestCase):
                 self.assertAlmostEqual(m_distances[i], c.distance.m, tol)
                 self.assertAlmostEqual(ft_distances[i], c.distance.survey_ft, tol)
 
+    @skipIf(spatialite, "distance method doesn't support geodetic coordinates on SpatiaLite.")
     @skipUnlessDBFeature("has_distance_method", "supports_distance_geodetic")
     @ignore_warnings(category=RemovedInDjango20Warning)
     def test_distance_geodetic(self):
@@ -270,12 +273,15 @@ class DistanceTest(TestCase):
         # a 100km of that line (which should exclude only Hobart & Adelaide).
         line = GEOSGeometry('LINESTRING(144.9630 -37.8143,151.2607 -33.8870)', 4326)
         dist_qs = AustraliaCity.objects.filter(point__distance_lte=(line, D(km=100)))
-
-        self.assertEqual(9, dist_qs.count())
-        self.assertEqual(['Batemans Bay', 'Canberra', 'Hillsdale',
-                          'Melbourne', 'Mittagong', 'Shellharbour',
-                          'Sydney', 'Thirroul', 'Wollongong'],
-                         self.get_names(dist_qs))
+        expected_cities = [
+            'Batemans Bay', 'Canberra', 'Hillsdale',
+            'Melbourne', 'Mittagong', 'Shellharbour',
+            'Sydney', 'Thirroul', 'Wollongong',
+        ]
+        if spatialite:
+            # SpatiaLite is less accurate and returns 102.8km for Batemans Bay.
+            expected_cities.pop(0)
+        self.assertEqual(expected_cities, self.get_names(dist_qs))
 
         # Too many params (4 in this case) should raise a ValueError.
         queryset = AustraliaCity.objects.filter(point__distance_lte=('POINT(5 23)', D(km=100), 'spheroid', '4'))
@@ -355,6 +361,7 @@ class DistanceTest(TestCase):
         for i, z in enumerate(SouthTexasZipcode.objects.order_by('name').area()):
             self.assertAlmostEqual(area_sq_m[i], z.area.sq_m, tol)
 
+    @skipIf(spatialite, "length method doesn't support geodetic coordinates on SpatiaLite.")
     @skipUnlessDBFeature("has_length_method")
     @ignore_warnings(category=RemovedInDjango20Warning)
     def test_length(self):
@@ -544,8 +551,9 @@ class DistanceFunctionsTests(TestCase):
                      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)
+            # Testing equivalence to within a meter (kilometer on SpatiaLite).
+            tol = -3 if spatialite else 0
+            self.assertAlmostEqual(distance, city.distance.m, tol)
 
     @skipUnlessDBFeature("has_Distance_function", "supports_distance_geodetic")
     def test_distance_geodetic_spheroid(self):
@@ -575,7 +583,7 @@ class DistanceFunctionsTests(TestCase):
         ).order_by('id')
         for i, c in enumerate(qs):
             self.assertAlmostEqual(spheroid_distances[i], c.distance.m, tol)
-        if postgis:
+        if postgis or spatialite:
             # PostGIS uses sphere-only distances by default, testing these as well.
             qs = AustraliaCity.objects.exclude(id=hillsdale.id).annotate(
                 distance=Distance('point', hillsdale.point)
diff --git a/tests/gis_tests/geogapp/tests.py b/tests/gis_tests/geogapp/tests.py
index f1e6de2625..89a47405ce 100644
--- a/tests/gis_tests/geogapp/tests.py
+++ b/tests/gis_tests/geogapp/tests.py
@@ -4,18 +4,20 @@ Tests for geography support in PostGIS
 from __future__ import unicode_literals
 
 import os
-from unittest import skipUnless
+from unittest import skipIf, skipUnless
 
 from django.contrib.gis.db import models
 from django.contrib.gis.db.models.functions import Area, Distance
 from django.contrib.gis.measure import D
 from django.db import connection
 from django.db.models.functions import Cast
-from django.test import TestCase, ignore_warnings, skipUnlessDBFeature
+from django.test import (
+    TestCase, ignore_warnings, skipIfDBFeature, skipUnlessDBFeature,
+)
 from django.utils._os import upath
 from django.utils.deprecation import RemovedInDjango20Warning
 
-from ..utils import oracle, postgis
+from ..utils import oracle, postgis, spatialite
 from .models import City, County, Zipcode
 
 
@@ -27,6 +29,7 @@ class GeographyTest(TestCase):
         "Ensure geography features loaded properly."
         self.assertEqual(8, City.objects.count())
 
+    @skipIf(spatialite, "SpatiaLite doesn't support distance lookups with Distance objects.")
     @skipUnlessDBFeature("supports_distances_lookups", "supports_distance_geodetic")
     def test02_distance_lookup(self):
         "Testing GeoQuerySet distance lookup support on non-point geography fields."
@@ -42,6 +45,7 @@ class GeographyTest(TestCase):
         for cities in [cities1, cities2]:
             self.assertEqual(['Dallas', 'Houston', 'Oklahoma City'], cities)
 
+    @skipIf(spatialite, "distance() doesn't support geodetic coordinates on SpatiaLite.")
     @skipUnlessDBFeature("has_distance_method", "supports_distance_geodetic")
     @ignore_warnings(category=RemovedInDjango20Warning)
     def test03_distance_method(self):
@@ -97,6 +101,7 @@ class GeographyTest(TestCase):
             self.assertEqual(name, c.name)
             self.assertEqual(state, c.state)
 
+    @skipIf(spatialite, "area() doesn't support geodetic coordinates on SpatiaLite.")
     @skipUnlessDBFeature("has_area_method", "supports_distance_geodetic")
     @ignore_warnings(category=RemovedInDjango20Warning)
     def test06_geography_area(self):
@@ -136,17 +141,22 @@ class GeographyFunctionTests(TestCase):
         """
         if oracle:
             ref_dists = [0, 4899.68, 8081.30, 9115.15]
+        elif spatialite:
+            # SpatiaLite returns non-zero distance for polygons and points
+            # covered by that polygon.
+            ref_dists = [326.61, 4899.68, 8081.30, 9115.15]
         else:
             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)
-        # Distance function in combination with a lookup.
-        hzip = Zipcode.objects.get(code='77002')
-        self.assertEqual(qs.get(distance__lte=0), hzip)
+        if not spatialite:
+            # Distance function combined with a lookup.
+            hzip = Zipcode.objects.get(code='77002')
+            self.assertEqual(qs.get(distance__lte=0), hzip)
 
-    @skipUnlessDBFeature("has_Area_function", "supports_distance_geodetic")
+    @skipUnlessDBFeature("has_Area_function", "supports_area_geodetic")
     def test_geography_area(self):
         """
         Testing that Area calculations work on geography columns.
@@ -158,3 +168,9 @@ class GeographyFunctionTests(TestCase):
         rounded_value = z.area.sq_m
         rounded_value -= z.area.sq_m % 1000
         self.assertEqual(rounded_value, 5439000)
+
+    @skipUnlessDBFeature("has_Area_function")
+    @skipIfDBFeature("supports_area_geodetic")
+    def test_geodetic_area_raises_if_not_supported(self):
+        with self.assertRaisesMessage(NotImplementedError, 'Area on geodetic coordinate systems not supported.'):
+            Zipcode.objects.annotate(area=Area('poly')).get(code='77002')