1 import os, unittest

     2 from decimal import Decimal

     3 

     4 from django.db import connection

     5 from django.db.models import Q

     6 from django.contrib.gis.gdal import DataSource

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

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

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

    10 

    11 from models import AustraliaCity, Interstate, SouthTexasInterstate, \

    12     SouthTexasCity, SouthTexasCityFt, CensusZipcode, SouthTexasZipcode

    13 from data import au_cities, interstates, stx_interstates, stx_cities, stx_zips

    14 

    15 class DistanceTest(unittest.TestCase):

    16 

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

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

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

    20     # w/units in meters)

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

    22     # Another one for Australia

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

    24 

    25     def get_names(self, qs):

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

    27         cities.sort()

    28         return cities

    29 

    30     def test01_init(self):

    31         "Initialization of distance models."

    32 

    33         # Loading up the cities.

    34         def load_cities(city_model, data_tup):

    35             for name, x, y in data_tup:

    36                 city_model(name=name, point=Point(x, y, srid=4326)).save()

    37 

    38         def load_interstates(imodel, data_tup):

    39             for name, wkt in data_tup:

    40                 imodel(name=name, path=wkt).save()

    41 

    42         load_cities(SouthTexasCity, stx_cities)

    43         load_cities(SouthTexasCityFt, stx_cities)

    44         load_cities(AustraliaCity, au_cities)

    45 

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

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

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

    49 

    50         # Loading up the South Texas Zip Codes.

    51         for name, wkt in stx_zips:

    52             poly = GEOSGeometry(wkt, srid=4269)

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

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

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

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

    57 

    58         # Loading up the Interstates.

    59         load_interstates(Interstate, interstates)

    60         load_interstates(SouthTexasInterstate, stx_interstates)

    61 

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

    63         self.assertEqual(1, SouthTexasInterstate.objects.count())

    64 

    65     @no_spatialite

    66     def test02_dwithin(self):

    67         "Testing the `dwithin` lookup type."

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

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

    70         # approximate).

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

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

    73 

    74         # Expected cities for Australia and Texas.

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

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

    77 

    78         # Performing distance queries on two projected coordinate systems one

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

    80         for dist in tx_dists:

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

    82             else: dist1 = dist2 = dist

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

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

    85             for qs in qs1, qs2:

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

    87 

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

    89         for dist in au_dists:

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

    91             else: type_error = False

    92 

    93             if isinstance(dist, tuple):

    94                 if oracle: dist = dist[1]

    95                 else: dist = dist[0]

    96 

    97             # Creating the query set.

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

    99             if type_error:

   100                 # A ValueError should be raised on PostGIS when trying to pass

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

   102                 self.assertRaises(ValueError, AustraliaCity.objects.filter(point__dwithin=(self.au_pnt, dist)).count)

   103             else:

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

   105 

   106     def test03a_distance_method(self):

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

   108         # The point for La Grange, TX

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

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

   111         # using the provided raw SQL statements.

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

   113         m_distances = [147075.069813, 139630.198056, 140888.552826,

   114                        138809.684197, 158309.246259, 212183.594374,

   115                        70870.188967, 165337.758878, 139196.085105]

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

   117         # Oracle 11 thinks this is not a projected coordinate system, so it's s

   118         # not tested.

   119         ft_distances = [482528.79154625, 458103.408123001, 462231.860397575,

   120                         455411.438904354, 519386.252102563, 696139.009211594,

   121                         232513.278304279, 542445.630586414, 456679.155883207]

   122 

   123         # Testing using different variations of parameters and using models

   124         # with different projected coordinate systems.

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

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

   127         if spatialite or oracle:

   128             dist_qs = [dist1, dist2]

   129         else:

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

   131             dist4 = SouthTexasCityFt.objects.distance(lagrange)

   132             dist_qs = [dist1, dist2, dist3, dist4]

   133 

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

   135         # for Oracle.

   136         if oracle: tol = 2

   137         else: tol = 5

   138 

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

   140         for qs in dist_qs:

   141             for i, c in enumerate(qs):

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

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

   144 

   145     @no_spatialite

   146     def test03b_distance_method(self):

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

   148         if oracle: tol = 2

   149         else: tol = 5

   150 

   151         # Testing geodetic distance calculation with a non-point geometry

   152         # (a LineString of Wollongong and Shellharbour coords).

   153         ls = LineString( ( (150.902, -34.4245), (150.87, -34.5789) ) )

   154         if oracle or connection.ops.geography:

   155             # Reference query:

   156             #  SELECT ST_distance_sphere(point, ST_GeomFromText('LINESTRING(150.9020 -34.4245,150.8700 -34.5789)', 4326)) FROM distapp_australiacity ORDER BY name;

   157             distances = [1120954.92533513, 140575.720018241, 640396.662906304,

   158                          60580.9693849269, 972807.955955075, 568451.8357838,

   159                          40435.4335201384, 0, 68272.3896586844, 12375.0643697706, 0]

   160             qs = AustraliaCity.objects.distance(ls).order_by('name')

   161             for city, distance in zip(qs, distances):

   162                 # Testing equivalence to within a meter.

   163                 self.assertAlmostEqual(distance, city.distance.m, 0)

   164         else:

   165             # PostGIS 1.4 and below is limited to disance queries only

   166             # to/from point geometries, check for raising of ValueError.

   167             self.assertRaises(ValueError, AustraliaCity.objects.distance, ls)

   168             self.assertRaises(ValueError, AustraliaCity.objects.distance, ls.wkt)

   169 

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

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

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

   173         if connection.ops.postgis and connection.ops.proj_version_tuple() >= (4, 7, 0):

   174             # PROJ.4 versions 4.7+ have updated datums, and thus different

   175             # distance values.

   176             spheroid_distances = [60504.0628957201, 77023.9489850262, 49154.8867574404,

   177                                   90847.4358768573, 217402.811919332, 709599.234564757,

   178                                   640011.483550888, 7772.00667991925, 1047861.78619339,

   179                                   1165126.55236034]

   180             sphere_distances = [60580.9693849267, 77144.0435286473, 49199.4415344719,

   181                                 90804.7533823494, 217713.384600405, 709134.127242793,

   182                                 639828.157159169, 7786.82949717788, 1049204.06569028,

   183                                 1162623.7238134]

   184 

   185         else:

   186             spheroid_distances = [60504.0628825298, 77023.948962654, 49154.8867507115,

   187                                   90847.435881812, 217402.811862568, 709599.234619957,

   188                                   640011.483583758, 7772.00667666425, 1047861.7859506,

   189                                   1165126.55237647]

   190             sphere_distances = [60580.7612632291, 77143.7785056615, 49199.2725132184,

   191                                 90804.4414289463, 217712.63666124, 709131.691061906,

   192                                 639825.959074112, 7786.80274606706, 1049200.46122281,

   193                                 1162619.7297006]

   194 

   195         # Testing with spheroid distances first.

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

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

   198         for i, c in enumerate(qs):

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

   200         if postgis:

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

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

   203             for i, c in enumerate(qs):

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

   205 

   206     @no_oracle # Oracle already handles geographic distance calculation.

   207     def test03c_distance_method(self):

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

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

   210         # that is not a PointField (on PostGIS 1.4 and below).

   211         if not connection.ops.geography:

   212             self.assertRaises(ValueError, CensusZipcode.objects.distance, self.stx_pnt)

   213 

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

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

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

   217         # a non-geographic system.

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

   219 

   220         # Reference query:

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

   222         dists_m = [3553.30384972258, 1243.18391525602, 2186.15439472242]

   223 

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

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

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

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

   228         # however.

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

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

   231         ref_zips = ['77002', '77025', '77401']

   232 

   233         for buf in [buf1, buf2]:

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

   235             self.assertEqual(ref_zips, self.get_names(qs))

   236             for i, z in enumerate(qs):

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

   238 

   239     def test04_distance_lookups(self):

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

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

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

   243         # the `test02_dwithin` above).

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

   245 

   246         # Can't determine the units on SpatiaLite from PROJ.4 string, and

   247         # Oracle 11 incorrectly thinks it is not projected.

   248         if spatialite or oracle:

   249             dist_qs = (qs1,)

   250         else:

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

   252             dist_qs = (qs1, qs2)

   253 

   254         for qs in dist_qs:

   255             cities = self.get_names(qs)

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

   257 

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

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

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

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

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

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

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

   265 

   266     def test05_geodetic_distance_lookups(self):

   267         "Testing distance lookups on geodetic coordinate systems."

   268         # Line is from Canberra to Sydney.  Query is for all other cities within

   269         # a 100km of that line (which should exclude only Hobart & Adelaide).

   270         line = GEOSGeometry('LINESTRING(144.9630 -37.8143,151.2607 -33.8870)', 4326)

   271         dist_qs = AustraliaCity.objects.filter(point__distance_lte=(line, D(km=100)))

   272 

   273         if oracle or connection.ops.geography:

   274             # Oracle and PostGIS 1.5 can do distance lookups on arbitrary geometries.

   275             self.assertEqual(9, dist_qs.count())

   276             self.assertEqual(['Batemans Bay', 'Canberra', 'Hillsdale',

   277                               'Melbourne', 'Mittagong', 'Shellharbour',

   278                               'Sydney', 'Thirroul', 'Wollongong'],

   279                              self.get_names(dist_qs))

   280         else:

   281             # PostGIS 1.4 and below only allows geodetic distance queries (utilizing

   282             # ST_Distance_Sphere/ST_Distance_Spheroid) from Points to PointFields

   283             # on geometry columns.

   284             self.assertRaises(ValueError, dist_qs.count)

   285 

   286             # Ensured that a ValueError was raised, none of the rest of the test is

   287             # support on this backend, so bail now.

   288             if spatialite: return

   289 

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

   291         self.assertRaises(ValueError, len,

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

   293 

   294         # Not enough params should raise a ValueError.

   295         self.assertRaises(ValueError, len,

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

   297 

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

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

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

   301         cities = self.get_names(qs)

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

   303 

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

   305         # and using different units of distance.

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

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

   308 

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

   310         gq1 = Q(point__distance_lte=(wollongong.point, d1))

   311         gq2 = Q(point__distance_gte=(wollongong.point, d2))

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

   313 

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

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

   316         # in this test case).

   317         if postgis:

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

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

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

   321             querysets = [qs1, qs2]

   322         else:

   323             querysets = [qs1]

   324 

   325         for qs in querysets:

   326             cities = self.get_names(qs)

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

   328 

   329     def test06_area(self):

   330         "Testing the `area` GeoQuerySet method."

   331         # Reference queries:

   332         # SELECT ST_Area(poly) FROM distapp_southtexaszipcode;

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

   334         # Tolerance has to be lower for Oracle and differences

   335         # with GEOS 3.0.0RC4

   336         tol = 2

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

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

   339 

   340     def test07_length(self):

   341         "Testing the `length` GeoQuerySet method."

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

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

   344         len_m1 = 473504.769553813

   345         len_m2 = 4617.668

   346 

   347         if spatialite:

   348             # Does not support geodetic coordinate systems.

   349             self.assertRaises(ValueError, Interstate.objects.length)

   350         else:

   351             qs = Interstate.objects.length()

   352             if oracle: tol = 2

   353             else: tol = 5

   354             self.assertAlmostEqual(len_m1, qs[0].length.m, tol)

   355 

   356         # Now doing length on a projected coordinate system.

   357         i10 = SouthTexasInterstate.objects.length().get(name='I-10')

   358         self.assertAlmostEqual(len_m2, i10.length.m, 2)

   359 

   360     @no_spatialite

   361     def test08_perimeter(self):

   362         "Testing the `perimeter` GeoQuerySet method."

   363         # Reference query:

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

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

   366         if oracle: tol = 2

   367         else: tol = 7

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

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

   370 

   371         # Running on points; should return 0.

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

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

   374 

   375     def test09_measurement_null_fields(self):

   376         "Testing the measurement GeoQuerySet methods on fields with NULL values."

   377         # Creating SouthTexasZipcode w/NULL value.

   378         SouthTexasZipcode.objects.create(name='78212')

   379         # Performing distance/area queries against the NULL PolygonField,

   380         # and ensuring the result of the operations is None.

   381         htown = SouthTexasCity.objects.get(name='Downtown Houston')

   382         z = SouthTexasZipcode.objects.distance(htown.point).area().get(name='78212')

   383         self.assertEqual(None, z.distance)

   384         self.assertEqual(None, z.area)

   385 

   386 def suite():

   387     s = unittest.TestSuite()

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

   389     return s