1 import os, re, unittest

     2 from django.contrib.gis.db.models import Union, Extent3D

     3 from django.contrib.gis.geos import GEOSGeometry, Point, Polygon

     4 from django.contrib.gis.utils import LayerMapping, LayerMapError

     5 

     6 from models import City3D, Interstate2D, Interstate3D, \

     7     InterstateProj2D, InterstateProj3D, \

     8     Point2D, Point3D, MultiPoint3D, Polygon2D, Polygon3D

     9 

    10 data_path = os.path.realpath(os.path.join(os.path.dirname(__file__), '..', 'data'))

    11 city_file = os.path.join(data_path, 'cities', 'cities.shp')

    12 vrt_file = os.path.join(data_path, 'test_vrt', 'test_vrt.vrt')

    13 

    14 # The coordinates of each city, with Z values corresponding to their

    15 # altitude in meters.

    16 city_data = (

    17     ('Houston', (-95.363151, 29.763374, 18)),

    18     ('Dallas', (-96.801611, 32.782057, 147)),

    19     ('Oklahoma City', (-97.521157, 34.464642, 380)),

    20     ('Wellington', (174.783117, -41.315268, 14)),

    21     ('Pueblo', (-104.609252, 38.255001, 1433)),

    22     ('Lawrence', (-95.235060, 38.971823, 251)),

    23     ('Chicago', (-87.650175, 41.850385, 181)),

    24     ('Victoria', (-123.305196, 48.462611, 15)),

    25 )

    26 

    27 # Reference mapping of city name to its altitude (Z value).

    28 city_dict = dict((name, coords) for name, coords in city_data)

    29 

    30 # 3D freeway data derived from the National Elevation Dataset: 

    31 #  http://seamless.usgs.gov/products/9arc.php

    32 interstate_data = (

    33     ('I-45', 

    34      'LINESTRING(-95.3708481 29.7765870 11.339,-95.3694580 29.7787980 4.536,-95.3690305 29.7797359 9.762,-95.3691886 29.7812450 12.448,-95.3696447 29.7850144 10.457,-95.3702511 29.7868518 9.418,-95.3706724 29.7881286 14.858,-95.3711632 29.7896157 15.386,-95.3714525 29.7936267 13.168,-95.3717848 29.7955007 15.104,-95.3717719 29.7969804 16.516,-95.3717305 29.7982117 13.923,-95.3717254 29.8000778 14.385,-95.3719875 29.8013539 15.160,-95.3720575 29.8026785 15.544,-95.3721321 29.8040912 14.975,-95.3722074 29.8050998 15.688,-95.3722779 29.8060430 16.099,-95.3733818 29.8076750 15.197,-95.3741563 29.8103686 17.268,-95.3749458 29.8129927 19.857,-95.3763564 29.8144557 15.435)',

    35      ( 11.339,   4.536,   9.762,  12.448,  10.457,   9.418,  14.858,

    36        15.386,  13.168,  15.104,  16.516,  13.923,  14.385,  15.16 ,

    37        15.544,  14.975,  15.688,  16.099,  15.197,  17.268,  19.857,

    38        15.435),

    39      ),

    40     )

    41 

    42 # Bounding box polygon for inner-loop of Houston (in projected coordinate

    43 # system 32140), with elevation values from the National Elevation Dataset

    44 # (see above).

    45 bbox_wkt = 'POLYGON((941527.97 4225693.20,962596.48 4226349.75,963152.57 4209023.95,942051.75 4208366.38,941527.97 4225693.20))'

    46 bbox_z = (21.71, 13.21, 9.12, 16.40, 21.71)

    47 def gen_bbox():

    48     bbox_2d = GEOSGeometry(bbox_wkt, srid=32140)

    49     bbox_3d = Polygon(tuple((x, y, z) for (x, y), z in zip(bbox_2d[0].coords, bbox_z)), srid=32140)    

    50     return bbox_2d, bbox_3d

    51 

    52 class Geo3DTest(unittest.TestCase):

    53     """

    54     Only a subset of the PostGIS routines are 3D-enabled, and this TestCase

    55     tries to test the features that can handle 3D and that are also 

    56     available within GeoDjango.  For more information, see the PostGIS docs

    57     on the routines that support 3D:

    58 

    59     http://postgis.refractions.net/documentation/manual-1.4/ch08.html#PostGIS_3D_Functions

    60     """

    61 

    62     def test01_3d(self):

    63         "Test the creation of 3D models."

    64         # 3D models for the rest of the tests will be populated in here.

    65         # For each 3D data set create model (and 2D version if necessary), 

    66         # retrieve, and assert geometry is in 3D and contains the expected

    67         # 3D values.

    68         for name, pnt_data in city_data:

    69             x, y, z = pnt_data

    70             pnt = Point(x, y, z, srid=4326)

    71             City3D.objects.create(name=name, point=pnt)

    72             city = City3D.objects.get(name=name)

    73             self.failUnless(city.point.hasz)

    74             self.assertEqual(z, city.point.z)

    75 

    76         # Interstate (2D / 3D and Geographic/Projected variants)

    77         for name, line, exp_z in interstate_data:

    78             line_3d = GEOSGeometry(line, srid=4269)

    79             # Using `hex` attribute because it omits 3D.

    80             line_2d = GEOSGeometry(line_3d.hex, srid=4269)

    81 

    82             # Creating a geographic and projected version of the

    83             # interstate in both 2D and 3D.

    84             Interstate3D.objects.create(name=name, line=line_3d)

    85             InterstateProj3D.objects.create(name=name, line=line_3d)

    86             Interstate2D.objects.create(name=name, line=line_2d)

    87             InterstateProj2D.objects.create(name=name, line=line_2d)

    88 

    89             # Retrieving and making sure it's 3D and has expected

    90             # Z values -- shouldn't change because of coordinate system.

    91             interstate = Interstate3D.objects.get(name=name)

    92             interstate_proj = InterstateProj3D.objects.get(name=name)

    93             for i in [interstate, interstate_proj]:

    94                 self.failUnless(i.line.hasz)

    95                 self.assertEqual(exp_z, tuple(i.line.z))

    96 

    97         # Creating 3D Polygon.

    98         bbox2d, bbox3d = gen_bbox()

    99         Polygon2D.objects.create(name='2D BBox', poly=bbox2d)

   100         Polygon3D.objects.create(name='3D BBox', poly=bbox3d)

   101         p3d = Polygon3D.objects.get(name='3D BBox')

   102         self.failUnless(p3d.poly.hasz)

   103         self.assertEqual(bbox3d, p3d.poly)

   104 

   105     def test01a_3d_layermapping(self):

   106         "Testing LayerMapping on 3D models."

   107         from models import Point2D, Point3D

   108 

   109         point_mapping = {'point' : 'POINT'}

   110         mpoint_mapping = {'mpoint' : 'MULTIPOINT'}

   111 

   112         # The VRT is 3D, but should still be able to map sans the Z.

   113         lm = LayerMapping(Point2D, vrt_file, point_mapping, transform=False)

   114         lm.save()

   115         self.assertEqual(3, Point2D.objects.count())

   116 

   117         # The city shapefile is 2D, and won't be able to fill the coordinates

   118         # in the 3D model -- thus, a LayerMapError is raised.

   119         self.assertRaises(LayerMapError, LayerMapping,

   120                           Point3D, city_file, point_mapping, transform=False)

   121         

   122         # 3D model should take 3D data just fine.

   123         lm = LayerMapping(Point3D, vrt_file, point_mapping, transform=False)

   124         lm.save()

   125         self.assertEqual(3, Point3D.objects.count())

   126 

   127         # Making sure LayerMapping.make_multi works right, by converting

   128         # a Point25D into a MultiPoint25D.

   129         lm = LayerMapping(MultiPoint3D, vrt_file, mpoint_mapping, transform=False)

   130         lm.save()

   131         self.assertEqual(3, MultiPoint3D.objects.count())

   132 

   133     def test02a_kml(self):

   134         "Test GeoQuerySet.kml() with Z values."

   135         h = City3D.objects.kml(precision=6).get(name='Houston')

   136         # KML should be 3D.

   137         # `SELECT ST_AsKML(point, 6) FROM geo3d_city3d WHERE name = 'Houston';`

   138         ref_kml_regex = re.compile(r'^<Point><coordinates>-95.363\d+,29.763\d+,18</coordinates></Point>$')

   139         self.failUnless(ref_kml_regex.match(h.kml))

   140 

   141     def test02b_geojson(self):

   142         "Test GeoQuerySet.geojson() with Z values."

   143         h = City3D.objects.geojson(precision=6).get(name='Houston')

   144         # GeoJSON should be 3D

   145         # `SELECT ST_AsGeoJSON(point, 6) FROM geo3d_city3d WHERE name='Houston';`

   146         ref_json_regex = re.compile(r'^{"type":"Point","coordinates":\[-95.363151,29.763374,18(\.0+)?\]}$')

   147         self.failUnless(ref_json_regex.match(h.geojson))

   148 

   149     def test03a_union(self):

   150         "Testing the Union aggregate of 3D models."

   151         # PostGIS query that returned the reference EWKT for this test:

   152         #  `SELECT ST_AsText(ST_Union(point)) FROM geo3d_city3d;`

   153         ref_ewkt = 'SRID=4326;MULTIPOINT(-123.305196 48.462611 15,-104.609252 38.255001 1433,-97.521157 34.464642 380,-96.801611 32.782057 147,-95.363151 29.763374 18,-95.23506 38.971823 251,-87.650175 41.850385 181,174.783117 -41.315268 14)'

   154         ref_union = GEOSGeometry(ref_ewkt)

   155         union = City3D.objects.aggregate(Union('point'))['point__union']

   156         self.failUnless(union.hasz)

   157         self.assertEqual(ref_union, union)

   158 

   159     def test03b_extent(self):

   160         "Testing the Extent3D aggregate for 3D models."

   161         # `SELECT ST_Extent3D(point) FROM geo3d_city3d;`

   162         ref_extent3d = (-123.305196, -41.315268, 14,174.783117, 48.462611, 1433)

   163         extent1 = City3D.objects.aggregate(Extent3D('point'))['point__extent3d']

   164         extent2 = City3D.objects.extent3d()

   165 

   166         def check_extent3d(extent3d, tol=6):

   167             for ref_val, ext_val in zip(ref_extent3d, extent3d):

   168                 self.assertAlmostEqual(ref_val, ext_val, tol)

   169 

   170         for e3d in [extent1, extent2]:

   171             check_extent3d(e3d)

   172 

   173     def test04_perimeter(self):

   174         "Testing GeoQuerySet.perimeter() on 3D fields."

   175         # Reference query for values below:

   176         #  `SELECT ST_Perimeter3D(poly), ST_Perimeter2D(poly) FROM geo3d_polygon3d;`

   177         ref_perim_3d = 76859.2620451

   178         ref_perim_2d = 76859.2577803

   179         tol = 6

   180         self.assertAlmostEqual(ref_perim_2d,

   181                                Polygon2D.objects.perimeter().get(name='2D BBox').perimeter.m,

   182                                tol)

   183         self.assertAlmostEqual(ref_perim_3d,

   184                                Polygon3D.objects.perimeter().get(name='3D BBox').perimeter.m,

   185                                tol)

   186 

   187     def test05_length(self):

   188         "Testing GeoQuerySet.length() on 3D fields."

   189         # ST_Length_Spheroid Z-aware, and thus does not need to use

   190         # a separate function internally.

   191         # `SELECT ST_Length_Spheroid(line, 'SPHEROID["GRS 1980",6378137,298.257222101]') 

   192         #    FROM geo3d_interstate[2d|3d];`

   193         tol = 3

   194         ref_length_2d = 4368.1721949481

   195         ref_length_3d = 4368.62547052088

   196         self.assertAlmostEqual(ref_length_2d,

   197                                Interstate2D.objects.length().get(name='I-45').length.m,

   198                                tol)

   199         self.assertAlmostEqual(ref_length_3d,

   200                                Interstate3D.objects.length().get(name='I-45').length.m,

   201                                tol)

   202 

   203         # Making sure `ST_Length3D` is used on for a projected

   204         # and 3D model rather than `ST_Length`.

   205         # `SELECT ST_Length(line) FROM geo3d_interstateproj2d;`

   206         ref_length_2d = 4367.71564892392

   207         # `SELECT ST_Length3D(line) FROM geo3d_interstateproj3d;`

   208         ref_length_3d = 4368.16897234101

   209         self.assertAlmostEqual(ref_length_2d,

   210                                InterstateProj2D.objects.length().get(name='I-45').length.m,

   211                                tol)

   212         self.assertAlmostEqual(ref_length_3d,

   213                                InterstateProj3D.objects.length().get(name='I-45').length.m,

   214                                tol)

   215         

   216     def test06_scale(self):

   217         "Testing GeoQuerySet.scale() on Z values."

   218         # Mapping of City name to reference Z values.

   219         zscales = (-3, 4, 23)

   220         for zscale in zscales:

   221             for city in City3D.objects.scale(1.0, 1.0, zscale):

   222                 self.assertEqual(city_dict[city.name][2] * zscale, city.scale.z)

   223 

   224     def test07_translate(self):

   225         "Testing GeoQuerySet.translate() on Z values."

   226         ztranslations = (5.23, 23, -17)

   227         for ztrans in ztranslations:

   228             for city in City3D.objects.translate(0, 0, ztrans):

   229                 self.assertEqual(city_dict[city.name][2] + ztrans, city.translate.z)

   230 

   231 def suite():

   232     s = unittest.TestSuite()

   233     s.addTest(unittest.makeSuite(Geo3DTest))

   234     return s