1 import os, unittest

     2 from decimal import Decimal

     3 

     4 from django.db.models import Q

     5 from django.contrib.gis.gdal import DataSource

     6 from django.contrib.gis.geos import GEOSGeometry, Point, LineString

     7 from django.contrib.gis.measure import D # alias for Distance

     8 from django.contrib.gis.db.models import GeoQ

     9 from django.contrib.gis.tests.utils import oracle, postgis, no_oracle

    10 

    11 from models import AustraliaCity, Interstate, SouthTexasCity, SouthTexasCityFt, CensusZipcode, SouthTexasZipcode

    12 from data import au_cities, interstates, stx_cities, stx_zips

    13 

    14 class DistanceTest(unittest.TestCase):

    15 

    16     # A point we are testing distances with -- using a WGS84

    17     # coordinate that'll be implicitly transormed to that to

    18     # the coordinate system of the field, EPSG:32140 (Texas South Central

    19     # w/units in meters)

    20     stx_pnt = GEOSGeometry('POINT (-95.370401017314293 29.704867409475465)', 4326)

    21     # Another one for Australia

    22     au_pnt = GEOSGeometry('POINT (150.791 -34.4919)', 4326)

    23 

    24     def get_names(self, qs):

    25         cities = [c.name for c in qs]

    26         cities.sort()

    27         return cities

    28 

    29     def test01_init(self):

    30         "Initialization of distance models."

    31 

    32         # Loading up the cities.

    33         def load_cities(city_model, data_tup):

    34             for name, x, y in data_tup:

    35                 c = city_model(name=name, point=Point(x, y, srid=4326))

    36                 c.save()

    37         

    38         load_cities(SouthTexasCity, stx_cities)

    39         load_cities(SouthTexasCityFt, stx_cities)

    40         load_cities(AustraliaCity, au_cities)

    41 

    42         self.assertEqual(9, SouthTexasCity.objects.count())

    43         self.assertEqual(9, SouthTexasCityFt.objects.count())

    44         self.assertEqual(11, AustraliaCity.objects.count())

    45         

    46         # Loading up the South Texas Zip Codes.

    47         for name, wkt in stx_zips:

    48             poly = GEOSGeometry(wkt, srid=4269)

    49             SouthTexasZipcode(name=name, poly=poly).save()

    50             CensusZipcode(name=name, poly=poly).save()

    51         self.assertEqual(4, SouthTexasZipcode.objects.count())

    52         self.assertEqual(4, CensusZipcode.objects.count())

    53 

    54         # Loading up the Interstates.

    55         for name, wkt in interstates:

    56             Interstate(name=name, line=GEOSGeometry(wkt, srid=4326)).save()

    57         self.assertEqual(1, Interstate.objects.count())

    58 

    59     def test02_dwithin(self):

    60         "Testing the `dwithin` lookup type."

    61         # Distances -- all should be equal (except for the

    62         # degree/meter pair in au_cities, that's somewhat

    63         # approximate).

    64         tx_dists = [(7000, 22965.83), D(km=7), D(mi=4.349)]

    65         au_dists = [(0.5, 32000), D(km=32), D(mi=19.884)]

    66         

    67         # Expected cities for Australia and Texas.

    68         tx_cities = ['Downtown Houston', 'Southside Place']

    69         au_cities = ['Mittagong', 'Shellharbour', 'Thirroul', 'Wollongong']

    70 

    71         # Performing distance queries on two projected coordinate systems one

    72         # with units in meters and the other in units of U.S. survey feet.

    73         for dist in tx_dists:

    74             if isinstance(dist, tuple): dist1, dist2 = dist

    75             else: dist1 = dist2 = dist

    76             qs1 = SouthTexasCity.objects.filter(point__dwithin=(self.stx_pnt, dist1))

    77             qs2 = SouthTexasCityFt.objects.filter(point__dwithin=(self.stx_pnt, dist2))

    78             for qs in qs1, qs2:

    79                 self.assertEqual(tx_cities, self.get_names(qs))

    80 

    81         # Now performing the `dwithin` queries on a geodetic coordinate system.

    82         for dist in au_dists:

    83             if isinstance(dist, D) and not oracle: type_error = True

    84             else: type_error = False

    85 

    86             if isinstance(dist, tuple):

    87                 if oracle: dist = dist[1]

    88                 else: dist = dist[0]

    89                 

    90             # Creating the query set.

    91             qs = AustraliaCity.objects.order_by('name')

    92             if type_error:

    93                 # A TypeError should be raised on PostGIS when trying to pass

    94                 # Distance objects into a DWithin query using a geodetic field.  

    95                 self.assertRaises(TypeError, AustraliaCity.objects.filter, point__dwithin=(self.au_pnt, dist))

    96             else:

    97                 self.assertEqual(au_cities, self.get_names(qs.filter(point__dwithin=(self.au_pnt, dist))))

    98                                  

    99     def test03a_distance_method(self):

   100         "Testing the `distance` GeoQuerySet method on projected coordinate systems."

   101         # The point for La Grange, TX

   102         lagrange = GEOSGeometry('POINT(-96.876369 29.905320)', 4326)

   103         # Reference distances in feet and in meters. Got these values from 

   104         # using the provided raw SQL statements.

   105         #  SELECT ST_Distance(point, ST_Transform(ST_GeomFromText('POINT(-96.876369 29.905320)', 4326), 32140)) FROM distapp_southtexascity;

   106         m_distances = [147075.069813, 139630.198056, 140888.552826,

   107                        138809.684197, 158309.246259, 212183.594374,

   108                        70870.188967, 165337.758878, 139196.085105]

   109         #  SELECT ST_Distance(point, ST_Transform(ST_GeomFromText('POINT(-96.876369 29.905320)', 4326), 2278)) FROM distapp_southtexascityft;

   110         ft_distances = [482528.79154625, 458103.408123001, 462231.860397575,

   111                         455411.438904354, 519386.252102563, 696139.009211594,

   112                         232513.278304279, 542445.630586414, 456679.155883207]

   113 

   114         # Testing using different variations of parameters and using models

   115         # with different projected coordinate systems.

   116         dist1 = SouthTexasCity.objects.distance(lagrange, field_name='point')

   117         dist2 = SouthTexasCity.objects.distance(lagrange)  # Using GEOSGeometry parameter

   118         dist3 = SouthTexasCityFt.objects.distance(lagrange.ewkt) # Using EWKT string parameter.

   119         dist4 = SouthTexasCityFt.objects.distance(lagrange)

   120 

   121         # Original query done on PostGIS, have to adjust AlmostEqual tolerance

   122         # for Oracle.

   123         if oracle: tol = 2

   124         else: tol = 5

   125 

   126         # Ensuring expected distances are returned for each distance queryset.

   127         for qs in [dist1, dist2, dist3, dist4]:

   128             for i, c in enumerate(qs):

   129                 self.assertAlmostEqual(m_distances[i], c.distance.m, tol)

   130                 self.assertAlmostEqual(ft_distances[i], c.distance.survey_ft, tol)

   131 

   132     def test03b_distance_method(self):

   133         "Testing the `distance` GeoQuerySet method on geodetic coordnate systems."

   134         if oracle: tol = 2

   135         else: tol = 5

   136 

   137         # Now testing geodetic distance aggregation.

   138         hillsdale = AustraliaCity.objects.get(name='Hillsdale')

   139         if not oracle:

   140             # PostGIS is limited to disance queries only to/from point geometries,

   141             # ensuring a TypeError is raised if something else is put in.

   142             self.assertRaises(TypeError, AustraliaCity.objects.distance, 'LINESTRING(0 0, 1 1)')

   143             self.assertRaises(TypeError, AustraliaCity.objects.distance, LineString((0, 0), (1, 1)))

   144 

   145         # Got the reference distances using the raw SQL statements:

   146         #  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));

   147         spheroid_distances = [60504.0628825298, 77023.948962654, 49154.8867507115, 90847.435881812, 217402.811862568, 709599.234619957, 640011.483583758, 7772.00667666425, 1047861.7859506, 1165126.55237647]

   148         #  SELECT ST_distance_sphere(point, ST_GeomFromText('POINT(151.231341 -33.952685)', 4326)) FROM distapp_australiacity WHERE (NOT (id = 11));  st_distance_sphere

   149         sphere_distances = [60580.7612632291, 77143.7785056615, 49199.2725132184, 90804.4414289463, 217712.63666124, 709131.691061906, 639825.959074112, 7786.80274606706, 1049200.46122281, 1162619.7297006]

   150 

   151         # Testing with spheroid distances first.

   152         qs = AustraliaCity.objects.exclude(id=hillsdale.id).distance(hillsdale.point, spheroid=True)

   153         for i, c in enumerate(qs):

   154             self.assertAlmostEqual(spheroid_distances[i], c.distance.m, tol)

   155         if postgis:

   156             # PostGIS uses sphere-only distances by default, testing these as well.

   157             qs =  AustraliaCity.objects.exclude(id=hillsdale.id).distance(hillsdale.point)

   158             for i, c in enumerate(qs):

   159                 self.assertAlmostEqual(sphere_distances[i], c.distance.m, tol)

   160 

   161     @no_oracle # Oracle already handles geographic distance calculation.

   162     def test03c_distance_method(self):

   163         "Testing the `distance` GeoQuerySet method used with `transform` on a geographic field."

   164         # Normally you can't compute distances from a geometry field

   165         # that is not a PointField (on PostGIS).

   166         self.assertRaises(TypeError, CensusZipcode.objects.distance, self.stx_pnt)

   167         

   168         # We'll be using a Polygon (created by buffering the centroid

   169         # of 77005 to 100m) -- which aren't allowed in geographic distance

   170         # queries normally, however our field has been transformed to

   171         # a non-geographic system.

   172         z = SouthTexasZipcode.objects.get(name='77005')

   173 

   174         # Reference query:

   175         # SELECT ST_Distance(ST_Transform("distapp_censuszipcode"."poly", 32140), ST_GeomFromText('<buffer_wkt>', 32140)) FROM "distapp_censuszipcode";

   176         dists_m = [3553.30384972258, 1243.18391525602, 2186.15439472242]

   177 

   178         # Having our buffer in the SRID of the transformation and of the field

   179         # -- should get the same results. The first buffer has no need for

   180         # transformation SQL because it is the same SRID as what was given

   181         # to `transform()`.  The second buffer will need to be transformed,

   182         # however.

   183         buf1 = z.poly.centroid.buffer(100)

   184         buf2 = buf1.transform(4269, clone=True)

   185         for buf in [buf1, buf2]:

   186             qs = CensusZipcode.objects.exclude(name='77005').transform(32140).distance(buf)

   187             self.assertEqual(['77002', '77025', '77401'], self.get_names(qs))

   188             for i, z in enumerate(qs):

   189                 self.assertAlmostEqual(z.distance.m, dists_m[i], 5)

   190 

   191     def test04_distance_lookups(self):

   192         "Testing the `distance_lt`, `distance_gt`, `distance_lte`, and `distance_gte` lookup types."

   193         # Retrieving the cities within a 20km 'donut' w/a 7km radius 'hole'

   194         # (thus, Houston and Southside place will be excluded as tested in

   195         # the `test02_dwithin` above).

   196         qs1 = SouthTexasCity.objects.filter(point__distance_gte=(self.stx_pnt, D(km=7))).filter(point__distance_lte=(self.stx_pnt, D(km=20)))

   197         qs2 = SouthTexasCityFt.objects.filter(point__distance_gte=(self.stx_pnt, D(km=7))).filter(point__distance_lte=(self.stx_pnt, D(km=20)))

   198         for qs in qs1, qs2:

   199             cities = self.get_names(qs)

   200             self.assertEqual(cities, ['Bellaire', 'Pearland', 'West University Place'])

   201 

   202         # Doing a distance query using Polygons instead of a Point.

   203         z = SouthTexasZipcode.objects.get(name='77005')

   204         qs = SouthTexasZipcode.objects.exclude(name='77005').filter(poly__distance_lte=(z.poly, D(m=275)))

   205         self.assertEqual(['77025', '77401'], self.get_names(qs))

   206         # If we add a little more distance 77002 should be included.

   207         qs = SouthTexasZipcode.objects.exclude(name='77005').filter(poly__distance_lte=(z.poly, D(m=300)))

   208         self.assertEqual(['77002', '77025', '77401'], self.get_names(qs))

   209         

   210     def test05_geodetic_distance_lookups(self):

   211         "Testing distance lookups on geodetic coordinate systems."

   212         if not oracle:

   213             # Oracle doesn't have this limitation -- PostGIS only allows geodetic

   214             # distance queries from Points to PointFields.

   215             mp = GEOSGeometry('MULTIPOINT(0 0, 5 23)')

   216             self.assertRaises(TypeError,

   217                               AustraliaCity.objects.filter(point__distance_lte=(mp, D(km=100))))

   218             # Too many params (4 in this case) should raise a ValueError.

   219             self.assertRaises(ValueError, 

   220                               AustraliaCity.objects.filter, point__distance_lte=('POINT(5 23)', D(km=100), 'spheroid', '4'))

   221 

   222         # Not enough params should raise a ValueError.

   223         self.assertRaises(ValueError,

   224                           AustraliaCity.objects.filter, point__distance_lte=('POINT(5 23)',))

   225 

   226         # Getting all cities w/in 550 miles of Hobart.

   227         hobart = AustraliaCity.objects.get(name='Hobart')

   228         qs = AustraliaCity.objects.exclude(name='Hobart').filter(point__distance_lte=(hobart.point, D(mi=550)))

   229         cities = self.get_names(qs)

   230         self.assertEqual(cities, ['Batemans Bay', 'Canberra', 'Melbourne'])

   231 

   232         # Cities that are either really close or really far from Wollongong --

   233         # and using different units of distance.

   234         wollongong = AustraliaCity.objects.get(name='Wollongong')

   235         d1, d2 = D(yd=19500), D(nm=400) # Yards (~17km) & Nautical miles.

   236 

   237         # Normal geodetic distance lookup (uses `distance_sphere` on PostGIS.

   238         gq1 = GeoQ(point__distance_lte=(wollongong.point, d1))

   239         gq2 = GeoQ(point__distance_gte=(wollongong.point, d2))

   240         qs1 = AustraliaCity.objects.exclude(name='Wollongong').filter(gq1 | gq2)

   241 

   242         # Geodetic distance lookup but telling GeoDjango to use `distance_spheroid`

   243         # instead (we should get the same results b/c accuracy variance won't matter

   244         # in this test case). Using `Q` instead of `GeoQ` to be different (post-qsrf

   245         # it doesn't matter).

   246         if postgis:

   247             gq3 = Q(point__distance_lte=(wollongong.point, d1, 'spheroid'))

   248             gq4 = Q(point__distance_gte=(wollongong.point, d2, 'spheroid'))

   249             qs2 = AustraliaCity.objects.exclude(name='Wollongong').filter(gq3 | gq4)

   250             querysets = [qs1, qs2]

   251         else:

   252             querysets = [qs1]

   253 

   254         for qs in querysets:

   255             cities = self.get_names(qs)

   256             self.assertEqual(cities, ['Adelaide', 'Hobart', 'Shellharbour', 'Thirroul'])

   257 

   258     def test06_area(self):

   259         "Testing the `area` GeoQuerySet method."

   260         # Reference queries:

   261         # SELECT ST_Area(poly) FROM distapp_southtexaszipcode;

   262         area_sq_m = [5437908.90234375, 10183031.4389648, 11254471.0073242, 9881708.91772461]

   263         # Tolerance has to be lower for Oracle and differences

   264         # with GEOS 3.0.0RC4

   265         tol = 2

   266         for i, z in enumerate(SouthTexasZipcode.objects.area()):

   267             self.assertAlmostEqual(area_sq_m[i], z.area.sq_m, tol)

   268 

   269     def test07_length(self):

   270         "Testing the `length` GeoQuerySet method."

   271         # Reference query (should use `length_spheroid`).

   272         # SELECT ST_length_spheroid(ST_GeomFromText('<wkt>', 4326) 'SPHEROID["WGS 84",6378137,298.257223563, AUTHORITY["EPSG","7030"]]');

   273         len_m = 473504.769553813

   274         qs = Interstate.objects.length()

   275         if oracle: tol = 2

   276         else: tol = 5

   277         self.assertAlmostEqual(len_m, qs[0].length.m, tol)

   278 

   279     def test08_perimeter(self):

   280         "Testing the `perimeter` GeoQuerySet method."

   281         # Reference query:

   282         # SELECT ST_Perimeter(distapp_southtexaszipcode.poly) FROM distapp_southtexaszipcode;

   283         perim_m = [18404.3550889361, 15627.2108551001, 20632.5588368978, 17094.5996143697]

   284         if oracle: tol = 2

   285         else: tol = 7

   286         for i, z in enumerate(SouthTexasZipcode.objects.perimeter()):

   287             self.assertAlmostEqual(perim_m[i], z.perimeter.m, tol)

   288 

   289         # Running on points; should return 0.

   290         for i, c in enumerate(SouthTexasCity.objects.perimeter(model_att='perim')):

   291             self.assertEqual(0, c.perim.m)

   292 

   293 def suite():

   294     s = unittest.TestSuite()

   295     s.addTest(unittest.makeSuite(DistanceTest))

   296     return s