295 lines
12 KiB
Python
295 lines
12 KiB
Python
import os
|
|
import re
|
|
|
|
from django.contrib.gis.db.models import Extent3D, Union
|
|
from django.contrib.gis.db.models.functions import (
|
|
AsGeoJSON, AsKML, Length, Perimeter, Scale, Translate,
|
|
)
|
|
from django.contrib.gis.geos import GEOSGeometry, LineString, Point, Polygon
|
|
from django.test import TestCase, skipUnlessDBFeature
|
|
from django.utils._os import upath
|
|
|
|
from .models import (
|
|
City3D, Interstate2D, Interstate3D, InterstateProj2D, InterstateProj3D,
|
|
MultiPoint3D, Point2D, Point3D, Polygon2D, Polygon3D,
|
|
)
|
|
|
|
data_path = os.path.realpath(os.path.join(os.path.dirname(upath(__file__)), '..', 'data'))
|
|
city_file = os.path.join(data_path, 'cities', 'cities.shp')
|
|
vrt_file = os.path.join(data_path, 'test_vrt', 'test_vrt.vrt')
|
|
|
|
# The coordinates of each city, with Z values corresponding to their
|
|
# altitude in meters.
|
|
city_data = (
|
|
('Houston', (-95.363151, 29.763374, 18)),
|
|
('Dallas', (-96.801611, 32.782057, 147)),
|
|
('Oklahoma City', (-97.521157, 34.464642, 380)),
|
|
('Wellington', (174.783117, -41.315268, 14)),
|
|
('Pueblo', (-104.609252, 38.255001, 1433)),
|
|
('Lawrence', (-95.235060, 38.971823, 251)),
|
|
('Chicago', (-87.650175, 41.850385, 181)),
|
|
('Victoria', (-123.305196, 48.462611, 15)),
|
|
)
|
|
|
|
# Reference mapping of city name to its altitude (Z value).
|
|
city_dict = {name: coords for name, coords in city_data}
|
|
|
|
# 3D freeway data derived from the National Elevation Dataset:
|
|
# http://seamless.usgs.gov/products/9arc.php
|
|
interstate_data = (
|
|
('I-45',
|
|
'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)',
|
|
(11.339, 4.536, 9.762, 12.448, 10.457, 9.418, 14.858,
|
|
15.386, 13.168, 15.104, 16.516, 13.923, 14.385, 15.16,
|
|
15.544, 14.975, 15.688, 16.099, 15.197, 17.268, 19.857,
|
|
15.435),
|
|
),
|
|
)
|
|
|
|
# Bounding box polygon for inner-loop of Houston (in projected coordinate
|
|
# system 32140), with elevation values from the National Elevation Dataset
|
|
# (see above).
|
|
bbox_data = (
|
|
'POLYGON((941527.97 4225693.20,962596.48 4226349.75,963152.57 4209023.95,'
|
|
'942051.75 4208366.38,941527.97 4225693.20))',
|
|
(21.71, 13.21, 9.12, 16.40, 21.71)
|
|
)
|
|
|
|
|
|
class Geo3DLoadingHelper:
|
|
def _load_interstate_data(self):
|
|
# Interstate (2D / 3D and Geographic/Projected variants)
|
|
for name, line, exp_z in interstate_data:
|
|
line_3d = GEOSGeometry(line, srid=4269)
|
|
line_2d = LineString([l[:2] for l in line_3d.coords], srid=4269)
|
|
|
|
# Creating a geographic and projected version of the
|
|
# interstate in both 2D and 3D.
|
|
Interstate3D.objects.create(name=name, line=line_3d)
|
|
InterstateProj3D.objects.create(name=name, line=line_3d)
|
|
Interstate2D.objects.create(name=name, line=line_2d)
|
|
InterstateProj2D.objects.create(name=name, line=line_2d)
|
|
|
|
def _load_city_data(self):
|
|
for name, pnt_data in city_data:
|
|
City3D.objects.create(
|
|
name=name, point=Point(*pnt_data, srid=4326), pointg=Point(*pnt_data, srid=4326),
|
|
)
|
|
|
|
def _load_polygon_data(self):
|
|
bbox_wkt, bbox_z = bbox_data
|
|
bbox_2d = GEOSGeometry(bbox_wkt, srid=32140)
|
|
bbox_3d = Polygon(tuple((x, y, z) for (x, y), z in zip(bbox_2d[0].coords, bbox_z)), srid=32140)
|
|
Polygon2D.objects.create(name='2D BBox', poly=bbox_2d)
|
|
Polygon3D.objects.create(name='3D BBox', poly=bbox_3d)
|
|
|
|
|
|
@skipUnlessDBFeature("gis_enabled", "supports_3d_storage")
|
|
class Geo3DTest(Geo3DLoadingHelper, TestCase):
|
|
"""
|
|
Only a subset of the PostGIS routines are 3D-enabled, and this TestCase
|
|
tries to test the features that can handle 3D and that are also
|
|
available within GeoDjango. For more information, see the PostGIS docs
|
|
on the routines that support 3D:
|
|
|
|
http://postgis.net/docs/PostGIS_Special_Functions_Index.html#PostGIS_3D_Functions
|
|
"""
|
|
|
|
def test_3d_hasz(self):
|
|
"""
|
|
Make sure data is 3D and has expected Z values -- shouldn't change
|
|
because of coordinate system.
|
|
"""
|
|
self._load_interstate_data()
|
|
for name, line, exp_z in interstate_data:
|
|
interstate = Interstate3D.objects.get(name=name)
|
|
interstate_proj = InterstateProj3D.objects.get(name=name)
|
|
for i in [interstate, interstate_proj]:
|
|
self.assertTrue(i.line.hasz)
|
|
self.assertEqual(exp_z, tuple(i.line.z))
|
|
|
|
self._load_city_data()
|
|
for name, pnt_data in city_data:
|
|
city = City3D.objects.get(name=name)
|
|
# Testing both geometry and geography fields
|
|
self.assertTrue(city.point.hasz)
|
|
self.assertTrue(city.pointg.hasz)
|
|
self.assertEqual(city.point.z, pnt_data[2])
|
|
self.assertEqual(city.pointg.z, pnt_data[2])
|
|
|
|
def test_3d_polygons(self):
|
|
"""
|
|
Test the creation of polygon 3D models.
|
|
"""
|
|
self._load_polygon_data()
|
|
p3d = Polygon3D.objects.get(name='3D BBox')
|
|
self.assertTrue(p3d.poly.hasz)
|
|
self.assertIsInstance(p3d.poly, Polygon)
|
|
self.assertEqual(p3d.poly.srid, 32140)
|
|
|
|
def test_3d_layermapping(self):
|
|
"""
|
|
Testing LayerMapping on 3D models.
|
|
"""
|
|
# Import here as GDAL is required for those imports
|
|
from django.contrib.gis.utils import LayerMapping, LayerMapError
|
|
|
|
point_mapping = {'point': 'POINT'}
|
|
mpoint_mapping = {'mpoint': 'MULTIPOINT'}
|
|
|
|
# The VRT is 3D, but should still be able to map sans the Z.
|
|
lm = LayerMapping(Point2D, vrt_file, point_mapping, transform=False)
|
|
lm.save()
|
|
self.assertEqual(3, Point2D.objects.count())
|
|
|
|
# The city shapefile is 2D, and won't be able to fill the coordinates
|
|
# in the 3D model -- thus, a LayerMapError is raised.
|
|
with self.assertRaises(LayerMapError):
|
|
LayerMapping(Point3D, city_file, point_mapping, transform=False)
|
|
|
|
# 3D model should take 3D data just fine.
|
|
lm = LayerMapping(Point3D, vrt_file, point_mapping, transform=False)
|
|
lm.save()
|
|
self.assertEqual(3, Point3D.objects.count())
|
|
|
|
# Making sure LayerMapping.make_multi works right, by converting
|
|
# a Point25D into a MultiPoint25D.
|
|
lm = LayerMapping(MultiPoint3D, vrt_file, mpoint_mapping, transform=False)
|
|
lm.save()
|
|
self.assertEqual(3, MultiPoint3D.objects.count())
|
|
|
|
@skipUnlessDBFeature("supports_3d_functions")
|
|
def test_union(self):
|
|
"""
|
|
Testing the Union aggregate of 3D models.
|
|
"""
|
|
# PostGIS query that returned the reference EWKT for this test:
|
|
# `SELECT ST_AsText(ST_Union(point)) FROM geo3d_city3d;`
|
|
self._load_city_data()
|
|
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)'
|
|
)
|
|
ref_union = GEOSGeometry(ref_ewkt)
|
|
union = City3D.objects.aggregate(Union('point'))['point__union']
|
|
self.assertTrue(union.hasz)
|
|
# Ordering of points in the resulting geometry may vary between implementations
|
|
self.assertSetEqual({p.ewkt for p in ref_union}, {p.ewkt for p in union})
|
|
|
|
@skipUnlessDBFeature("supports_3d_functions")
|
|
def test_extent(self):
|
|
"""
|
|
Testing the Extent3D aggregate for 3D models.
|
|
"""
|
|
self._load_city_data()
|
|
# `SELECT ST_Extent3D(point) FROM geo3d_city3d;`
|
|
ref_extent3d = (-123.305196, -41.315268, 14, 174.783117, 48.462611, 1433)
|
|
extent = City3D.objects.aggregate(Extent3D('point'))['point__extent3d']
|
|
|
|
def check_extent3d(extent3d, tol=6):
|
|
for ref_val, ext_val in zip(ref_extent3d, extent3d):
|
|
self.assertAlmostEqual(ref_val, ext_val, tol)
|
|
|
|
check_extent3d(extent)
|
|
self.assertIsNone(City3D.objects.none().aggregate(Extent3D('point'))['point__extent3d'])
|
|
|
|
|
|
@skipUnlessDBFeature("gis_enabled", "supports_3d_functions")
|
|
class Geo3DFunctionsTests(Geo3DLoadingHelper, TestCase):
|
|
def test_kml(self):
|
|
"""
|
|
Test KML() function with Z values.
|
|
"""
|
|
self._load_city_data()
|
|
h = City3D.objects.annotate(kml=AsKML('point', precision=6)).get(name='Houston')
|
|
# KML should be 3D.
|
|
# `SELECT ST_AsKML(point, 6) FROM geo3d_city3d WHERE name = 'Houston';`
|
|
ref_kml_regex = re.compile(r'^<Point><coordinates>-95.363\d+,29.763\d+,18</coordinates></Point>$')
|
|
self.assertTrue(ref_kml_regex.match(h.kml))
|
|
|
|
def test_geojson(self):
|
|
"""
|
|
Test GeoJSON() function with Z values.
|
|
"""
|
|
self._load_city_data()
|
|
h = City3D.objects.annotate(geojson=AsGeoJSON('point', precision=6)).get(name='Houston')
|
|
# GeoJSON should be 3D
|
|
# `SELECT ST_AsGeoJSON(point, 6) FROM geo3d_city3d WHERE name='Houston';`
|
|
ref_json_regex = re.compile(r'^{"type":"Point","coordinates":\[-95.363151,29.763374,18(\.0+)?\]}$')
|
|
self.assertTrue(ref_json_regex.match(h.geojson))
|
|
|
|
def test_perimeter(self):
|
|
"""
|
|
Testing Perimeter() function on 3D fields.
|
|
"""
|
|
self._load_polygon_data()
|
|
# Reference query for values below:
|
|
# `SELECT ST_Perimeter3D(poly), ST_Perimeter2D(poly) FROM geo3d_polygon3d;`
|
|
ref_perim_3d = 76859.2620451
|
|
ref_perim_2d = 76859.2577803
|
|
tol = 6
|
|
poly2d = Polygon2D.objects.annotate(perimeter=Perimeter('poly')).get(name='2D BBox')
|
|
self.assertAlmostEqual(ref_perim_2d, poly2d.perimeter.m, tol)
|
|
poly3d = Polygon3D.objects.annotate(perimeter=Perimeter('poly')).get(name='3D BBox')
|
|
self.assertAlmostEqual(ref_perim_3d, poly3d.perimeter.m, tol)
|
|
|
|
def test_length(self):
|
|
"""
|
|
Testing Length() function on 3D fields.
|
|
"""
|
|
# ST_Length_Spheroid Z-aware, and thus does not need to use
|
|
# a separate function internally.
|
|
# `SELECT ST_Length_Spheroid(line, 'SPHEROID["GRS 1980",6378137,298.257222101]')
|
|
# FROM geo3d_interstate[2d|3d];`
|
|
self._load_interstate_data()
|
|
tol = 3
|
|
ref_length_2d = 4368.1721949481
|
|
ref_length_3d = 4368.62547052088
|
|
inter2d = Interstate2D.objects.annotate(length=Length('line')).get(name='I-45')
|
|
self.assertAlmostEqual(ref_length_2d, inter2d.length.m, tol)
|
|
inter3d = Interstate3D.objects.annotate(length=Length('line')).get(name='I-45')
|
|
self.assertAlmostEqual(ref_length_3d, inter3d.length.m, tol)
|
|
|
|
# Making sure `ST_Length3D` is used on for a projected
|
|
# and 3D model rather than `ST_Length`.
|
|
# `SELECT ST_Length(line) FROM geo3d_interstateproj2d;`
|
|
ref_length_2d = 4367.71564892392
|
|
# `SELECT ST_Length3D(line) FROM geo3d_interstateproj3d;`
|
|
ref_length_3d = 4368.16897234101
|
|
inter2d = InterstateProj2D.objects.annotate(length=Length('line')).get(name='I-45')
|
|
self.assertAlmostEqual(ref_length_2d, inter2d.length.m, tol)
|
|
inter3d = InterstateProj3D.objects.annotate(length=Length('line')).get(name='I-45')
|
|
self.assertAlmostEqual(ref_length_3d, inter3d.length.m, tol)
|
|
|
|
def test_scale(self):
|
|
"""
|
|
Testing Scale() function on Z values.
|
|
"""
|
|
self._load_city_data()
|
|
# Mapping of City name to reference Z values.
|
|
zscales = (-3, 4, 23)
|
|
for zscale in zscales:
|
|
for city in City3D.objects.annotate(scale=Scale('point', 1.0, 1.0, zscale)):
|
|
self.assertEqual(city_dict[city.name][2] * zscale, city.scale.z)
|
|
|
|
def test_translate(self):
|
|
"""
|
|
Testing Translate() function on Z values.
|
|
"""
|
|
self._load_city_data()
|
|
ztranslations = (5.23, 23, -17)
|
|
for ztrans in ztranslations:
|
|
for city in City3D.objects.annotate(translate=Translate('point', 0, 0, ztrans)):
|
|
self.assertEqual(city_dict[city.name][2] + ztrans, city.translate.z)
|