PageRenderTime 88ms CodeModel.GetById 73ms app.highlight 12ms RepoModel.GetById 0ms app.codeStats 1ms

/django/contrib/gis/tests/geo3d/tests.py

https://code.google.com/p/mango-py/
Python | 231 lines | 151 code | 29 blank | 51 comment | 22 complexity | 19832037a00b73be7de890a623371492 MD5 | raw file
  1import os
  2import re
  3from django.utils.unittest import TestCase
  4from django.contrib.gis.db.models import Union, Extent3D
  5from django.contrib.gis.geos import GEOSGeometry, Point, Polygon
  6from django.contrib.gis.utils import LayerMapping, LayerMapError
  7
  8from models import City3D, Interstate2D, Interstate3D, \
  9    InterstateProj2D, InterstateProj3D, \
 10    Point2D, Point3D, MultiPoint3D, Polygon2D, Polygon3D
 11
 12data_path = os.path.realpath(os.path.join(os.path.dirname(__file__), '..', 'data'))
 13city_file = os.path.join(data_path, 'cities', 'cities.shp')
 14vrt_file = os.path.join(data_path, 'test_vrt', 'test_vrt.vrt')
 15
 16# The coordinates of each city, with Z values corresponding to their
 17# altitude in meters.
 18city_data = (
 19    ('Houston', (-95.363151, 29.763374, 18)),
 20    ('Dallas', (-96.801611, 32.782057, 147)),
 21    ('Oklahoma City', (-97.521157, 34.464642, 380)),
 22    ('Wellington', (174.783117, -41.315268, 14)),
 23    ('Pueblo', (-104.609252, 38.255001, 1433)),
 24    ('Lawrence', (-95.235060, 38.971823, 251)),
 25    ('Chicago', (-87.650175, 41.850385, 181)),
 26    ('Victoria', (-123.305196, 48.462611, 15)),
 27)
 28
 29# Reference mapping of city name to its altitude (Z value).
 30city_dict = dict((name, coords) for name, coords in city_data)
 31
 32# 3D freeway data derived from the National Elevation Dataset: 
 33#  http://seamless.usgs.gov/products/9arc.php
 34interstate_data = (
 35    ('I-45', 
 36     '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)',
 37     ( 11.339,   4.536,   9.762,  12.448,  10.457,   9.418,  14.858,
 38       15.386,  13.168,  15.104,  16.516,  13.923,  14.385,  15.16 ,
 39       15.544,  14.975,  15.688,  16.099,  15.197,  17.268,  19.857,
 40       15.435),
 41     ),
 42    )
 43
 44# Bounding box polygon for inner-loop of Houston (in projected coordinate
 45# system 32140), with elevation values from the National Elevation Dataset
 46# (see above).
 47bbox_wkt = 'POLYGON((941527.97 4225693.20,962596.48 4226349.75,963152.57 4209023.95,942051.75 4208366.38,941527.97 4225693.20))'
 48bbox_z = (21.71, 13.21, 9.12, 16.40, 21.71)
 49def gen_bbox():
 50    bbox_2d = GEOSGeometry(bbox_wkt, srid=32140)
 51    bbox_3d = Polygon(tuple((x, y, z) for (x, y), z in zip(bbox_2d[0].coords, bbox_z)), srid=32140)    
 52    return bbox_2d, bbox_3d
 53
 54class Geo3DTest(TestCase):
 55    """
 56    Only a subset of the PostGIS routines are 3D-enabled, and this TestCase
 57    tries to test the features that can handle 3D and that are also 
 58    available within GeoDjango.  For more information, see the PostGIS docs
 59    on the routines that support 3D:
 60
 61    http://postgis.refractions.net/documentation/manual-1.4/ch08.html#PostGIS_3D_Functions
 62    """
 63
 64    def test01_3d(self):
 65        "Test the creation of 3D models."
 66        # 3D models for the rest of the tests will be populated in here.
 67        # For each 3D data set create model (and 2D version if necessary), 
 68        # retrieve, and assert geometry is in 3D and contains the expected
 69        # 3D values.
 70        for name, pnt_data in city_data:
 71            x, y, z = pnt_data
 72            pnt = Point(x, y, z, srid=4326)
 73            City3D.objects.create(name=name, point=pnt)
 74            city = City3D.objects.get(name=name)
 75            self.assertTrue(city.point.hasz)
 76            self.assertEqual(z, city.point.z)
 77
 78        # Interstate (2D / 3D and Geographic/Projected variants)
 79        for name, line, exp_z in interstate_data:
 80            line_3d = GEOSGeometry(line, srid=4269)
 81            # Using `hex` attribute because it omits 3D.
 82            line_2d = GEOSGeometry(line_3d.hex, srid=4269)
 83
 84            # Creating a geographic and projected version of the
 85            # interstate in both 2D and 3D.
 86            Interstate3D.objects.create(name=name, line=line_3d)
 87            InterstateProj3D.objects.create(name=name, line=line_3d)
 88            Interstate2D.objects.create(name=name, line=line_2d)
 89            InterstateProj2D.objects.create(name=name, line=line_2d)
 90
 91            # Retrieving and making sure it's 3D and has expected
 92            # Z values -- shouldn't change because of coordinate system.
 93            interstate = Interstate3D.objects.get(name=name)
 94            interstate_proj = InterstateProj3D.objects.get(name=name)
 95            for i in [interstate, interstate_proj]:
 96                self.assertTrue(i.line.hasz)
 97                self.assertEqual(exp_z, tuple(i.line.z))
 98
 99        # Creating 3D Polygon.
100        bbox2d, bbox3d = gen_bbox()
101        Polygon2D.objects.create(name='2D BBox', poly=bbox2d)
102        Polygon3D.objects.create(name='3D BBox', poly=bbox3d)
103        p3d = Polygon3D.objects.get(name='3D BBox')
104        self.assertTrue(p3d.poly.hasz)
105        self.assertEqual(bbox3d, p3d.poly)
106
107    def test01a_3d_layermapping(self):
108        "Testing LayerMapping on 3D models."
109        from models import Point2D, Point3D
110
111        point_mapping = {'point' : 'POINT'}
112        mpoint_mapping = {'mpoint' : 'MULTIPOINT'}
113
114        # The VRT is 3D, but should still be able to map sans the Z.
115        lm = LayerMapping(Point2D, vrt_file, point_mapping, transform=False)
116        lm.save()
117        self.assertEqual(3, Point2D.objects.count())
118
119        # The city shapefile is 2D, and won't be able to fill the coordinates
120        # in the 3D model -- thus, a LayerMapError is raised.
121        self.assertRaises(LayerMapError, LayerMapping,
122                          Point3D, city_file, point_mapping, transform=False)
123        
124        # 3D model should take 3D data just fine.
125        lm = LayerMapping(Point3D, vrt_file, point_mapping, transform=False)
126        lm.save()
127        self.assertEqual(3, Point3D.objects.count())
128
129        # Making sure LayerMapping.make_multi works right, by converting
130        # a Point25D into a MultiPoint25D.
131        lm = LayerMapping(MultiPoint3D, vrt_file, mpoint_mapping, transform=False)
132        lm.save()
133        self.assertEqual(3, MultiPoint3D.objects.count())
134
135    def test02a_kml(self):
136        "Test GeoQuerySet.kml() with Z values."
137        h = City3D.objects.kml(precision=6).get(name='Houston')
138        # KML should be 3D.
139        # `SELECT ST_AsKML(point, 6) FROM geo3d_city3d WHERE name = 'Houston';`
140        ref_kml_regex = re.compile(r'^<Point><coordinates>-95.363\d+,29.763\d+,18</coordinates></Point>$')
141        self.assertTrue(ref_kml_regex.match(h.kml))
142
143    def test02b_geojson(self):
144        "Test GeoQuerySet.geojson() with Z values."
145        h = City3D.objects.geojson(precision=6).get(name='Houston')
146        # GeoJSON should be 3D
147        # `SELECT ST_AsGeoJSON(point, 6) FROM geo3d_city3d WHERE name='Houston';`
148        ref_json_regex = re.compile(r'^{"type":"Point","coordinates":\[-95.363151,29.763374,18(\.0+)?\]}$')
149        self.assertTrue(ref_json_regex.match(h.geojson))
150
151    def test03a_union(self):
152        "Testing the Union aggregate of 3D models."
153        # PostGIS query that returned the reference EWKT for this test:
154        #  `SELECT ST_AsText(ST_Union(point)) FROM geo3d_city3d;`
155        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)'
156        ref_union = GEOSGeometry(ref_ewkt)
157        union = City3D.objects.aggregate(Union('point'))['point__union']
158        self.assertTrue(union.hasz)
159        self.assertEqual(ref_union, union)
160
161    def test03b_extent(self):
162        "Testing the Extent3D aggregate for 3D models."
163        # `SELECT ST_Extent3D(point) FROM geo3d_city3d;`
164        ref_extent3d = (-123.305196, -41.315268, 14,174.783117, 48.462611, 1433)
165        extent1 = City3D.objects.aggregate(Extent3D('point'))['point__extent3d']
166        extent2 = City3D.objects.extent3d()
167
168        def check_extent3d(extent3d, tol=6):
169            for ref_val, ext_val in zip(ref_extent3d, extent3d):
170                self.assertAlmostEqual(ref_val, ext_val, tol)
171
172        for e3d in [extent1, extent2]:
173            check_extent3d(e3d)
174
175    def test04_perimeter(self):
176        "Testing GeoQuerySet.perimeter() on 3D fields."
177        # Reference query for values below:
178        #  `SELECT ST_Perimeter3D(poly), ST_Perimeter2D(poly) FROM geo3d_polygon3d;`
179        ref_perim_3d = 76859.2620451
180        ref_perim_2d = 76859.2577803
181        tol = 6
182        self.assertAlmostEqual(ref_perim_2d,
183                               Polygon2D.objects.perimeter().get(name='2D BBox').perimeter.m,
184                               tol)
185        self.assertAlmostEqual(ref_perim_3d,
186                               Polygon3D.objects.perimeter().get(name='3D BBox').perimeter.m,
187                               tol)
188
189    def test05_length(self):
190        "Testing GeoQuerySet.length() on 3D fields."
191        # ST_Length_Spheroid Z-aware, and thus does not need to use
192        # a separate function internally.
193        # `SELECT ST_Length_Spheroid(line, 'SPHEROID["GRS 1980",6378137,298.257222101]') 
194        #    FROM geo3d_interstate[2d|3d];`
195        tol = 3
196        ref_length_2d = 4368.1721949481
197        ref_length_3d = 4368.62547052088
198        self.assertAlmostEqual(ref_length_2d,
199                               Interstate2D.objects.length().get(name='I-45').length.m,
200                               tol)
201        self.assertAlmostEqual(ref_length_3d,
202                               Interstate3D.objects.length().get(name='I-45').length.m,
203                               tol)
204
205        # Making sure `ST_Length3D` is used on for a projected
206        # and 3D model rather than `ST_Length`.
207        # `SELECT ST_Length(line) FROM geo3d_interstateproj2d;`
208        ref_length_2d = 4367.71564892392
209        # `SELECT ST_Length3D(line) FROM geo3d_interstateproj3d;`
210        ref_length_3d = 4368.16897234101
211        self.assertAlmostEqual(ref_length_2d,
212                               InterstateProj2D.objects.length().get(name='I-45').length.m,
213                               tol)
214        self.assertAlmostEqual(ref_length_3d,
215                               InterstateProj3D.objects.length().get(name='I-45').length.m,
216                               tol)
217        
218    def test06_scale(self):
219        "Testing GeoQuerySet.scale() on Z values."
220        # Mapping of City name to reference Z values.
221        zscales = (-3, 4, 23)
222        for zscale in zscales:
223            for city in City3D.objects.scale(1.0, 1.0, zscale):
224                self.assertEqual(city_dict[city.name][2] * zscale, city.scale.z)
225
226    def test07_translate(self):
227        "Testing GeoQuerySet.translate() on Z values."
228        ztranslations = (5.23, 23, -17)
229        for ztrans in ztranslations:
230            for city in City3D.objects.translate(0, 0, ztrans):
231                self.assertEqual(city_dict[city.name][2] + ztrans, city.translate.z)