Fixed #11433 -- 3D geometry fields are now supported with PostGIS; EWKB is now…

Fixed #11433 -- 3D geometry fields are now supported with PostGIS; EWKB is now used by `PostGISAdaptor`.


git-svn-id: http://code.djangoproject.com/svn/django/trunk@11742 bcc190cf-cafb-0310-a4f2-bffc1f526a37

django/contrib/gis/db/backend/postgis/__init__.py

@@ -18,18 +18,21 @@ SpatialBackend = BaseSpatialBackend(name='postgis', postgis=True,
                                     distance_spheroid=DISTANCE_SPHEROID,
                                     envelope=ENVELOPE,
                                     extent=EXTENT,
+                                    extent3d=EXTENT3D,
                                     gis_terms=POSTGIS_TERMS,
                                     geojson=ASGEOJSON,
                                     gml=ASGML,
                                     intersection=INTERSECTION,
                                     kml=ASKML,
                                     length=LENGTH,
+                                    length3d=LENGTH3D,
                                     length_spheroid=LENGTH_SPHEROID,
                                     make_line=MAKE_LINE,
                                     mem_size=MEM_SIZE,
                                     num_geom=NUM_GEOM,
                                     num_points=NUM_POINTS,
                                     perimeter=PERIMETER,
+                                    perimeter3d=PERIMETER3D,
                                     point_on_surface=POINT_ON_SURFACE,
                                     scale=SCALE,
                                     select=GEOM_SELECT,

django/contrib/gis/db/backend/postgis/adaptor.py

@@ -2,7 +2,7 @@
  This object provides quoting for GEOS geometries into PostgreSQL/PostGIS.
 """
 
-from django.contrib.gis.db.backend.postgis.query import GEOM_FROM_WKB
+from django.contrib.gis.db.backend.postgis.query import GEOM_FROM_EWKB
 from psycopg2 import Binary
 from psycopg2.extensions import ISQLQuote
 
@@ -11,7 +11,7 @@ class PostGISAdaptor(object):
         "Initializes on the geometry."
         # Getting the WKB (in string form, to allow easy pickling of
         # the adaptor) and the SRID from the geometry.
-        self.wkb = str(geom.wkb)
+        self.ewkb = str(geom.ewkb)
         self.srid = geom.srid
 
     def __conform__(self, proto):
@@ -30,7 +30,7 @@ class PostGISAdaptor(object):
     def getquoted(self):
         "Returns a properly quoted string for use in PostgreSQL/PostGIS."
         # Want to use WKB, so wrap with psycopg2 Binary() to quote properly.
-        return "%s(E%s, %s)" % (GEOM_FROM_WKB, Binary(self.wkb), self.srid or -1)
+        return "%s(E%s)" % (GEOM_FROM_EWKB, Binary(self.ewkb))
 
     def prepare_database_save(self, unused):
         return self

django/contrib/gis/db/backend/postgis/query.py

@@ -63,17 +63,21 @@ if MAJOR_VERSION >= 1:
     DISTANCE_SPHERE = get_func('distance_sphere')
     DISTANCE_SPHEROID = get_func('distance_spheroid')
     ENVELOPE = get_func('Envelope')
-    EXTENT = get_func('extent')
+    EXTENT = get_func('Extent')
+    EXTENT3D = get_func('Extent3D')
     GEOM_FROM_TEXT = get_func('GeomFromText')
+    GEOM_FROM_EWKB = get_func('GeomFromEWKB')
     GEOM_FROM_WKB = get_func('GeomFromWKB')
     INTERSECTION = get_func('Intersection')
     LENGTH = get_func('Length')
+    LENGTH3D = get_func('Length3D')
     LENGTH_SPHEROID = get_func('length_spheroid')
     MAKE_LINE = get_func('MakeLine')
     MEM_SIZE = get_func('mem_size')
     NUM_GEOM = get_func('NumGeometries')
     NUM_POINTS = get_func('npoints')
     PERIMETER = get_func('Perimeter')
+    PERIMETER3D = get_func('Perimeter3D')
     POINT_ON_SURFACE = get_func('PointOnSurface')
     SCALE = get_func('Scale')
     SNAP_TO_GRID = get_func('SnapToGrid')

django/contrib/gis/db/models/aggregates.py

@@ -24,6 +24,9 @@ class Collect(GeoAggregate):
 class Extent(GeoAggregate):
     name = 'Extent'
 
+class Extent3D(GeoAggregate):
+    name = 'Extent3D'
+
 class MakeLine(GeoAggregate):
     name = 'MakeLine'
 

django/contrib/gis/db/models/manager.py

@@ -34,6 +34,9 @@ class GeoManager(Manager):
     def extent(self, *args, **kwargs):
         return self.get_query_set().extent(*args, **kwargs)
 
+    def extent3d(self, *args, **kwargs):
+        return self.get_query_set().extent3d(*args, **kwargs)
+
     def geojson(self, *args, **kwargs):
         return self.get_query_set().geojson(*args, **kwargs)
 

django/contrib/gis/db/models/query.py

@@ -110,6 +110,14 @@ class GeoQuerySet(QuerySet):
         """
         return self._spatial_aggregate(aggregates.Extent, **kwargs)
 
+    def extent3d(self, **kwargs):
+        """
+        Returns the aggregate extent, in 3D, of the features in the
+        GeoQuerySet. It is returned as a 6-tuple, comprising:
+          (xmin, ymin, zmin, xmax, ymax, zmax).
+        """
+        return self._spatial_aggregate(aggregates.Extent3D, **kwargs)
+
     def geojson(self, precision=8, crs=False, bbox=False, **kwargs):
         """
         Returns a GeoJSON representation of the geomtry field in a `geojson`
@@ -524,12 +532,14 @@ class GeoQuerySet(QuerySet):
         else:
             dist_att = Distance.unit_attname(geo_field.units_name)
 
-        # Shortcut booleans for what distance function we're using.
+        # Shortcut booleans for what distance function we're using and
+        # whether the geometry field is 3D.
         distance = func == 'distance'
         length = func == 'length'
         perimeter = func == 'perimeter'
         if not (distance or length or perimeter):
             raise ValueError('Unknown distance function: %s' % func)
+        geom_3d = geo_field.dim == 3
 
         # The field's get_db_prep_lookup() is used to get any
         # extra distance parameters.  Here we set up the
@@ -604,7 +614,7 @@ class GeoQuerySet(QuerySet):
                     # some error checking is required.
                     if not isinstance(geo_field, PointField):
                         raise ValueError('Spherical distance calculation only supported on PointFields.')
-                    if not str(SpatialBackend.Geometry(buffer(params[0].wkb)).geom_type) == 'Point':
+                    if not str(SpatialBackend.Geometry(buffer(params[0].ewkb)).geom_type) == 'Point':
                         raise ValueError('Spherical distance calculation only supported with Point Geometry parameters')
                     # The `function` procedure argument needs to be set differently for
                     # geodetic distance calculations.
@@ -617,9 +627,16 @@ class GeoQuerySet(QuerySet):
             elif length or perimeter:
                 procedure_fmt = '%(geo_col)s'
                 if geodetic and length:
-                    # There's no `length_sphere`
+                    # There's no `length_sphere`, and `length_spheroid` also
+                    # works on 3D geometries.
                     procedure_fmt += ',%(spheroid)s'
                     procedure_args.update({'function' : SpatialBackend.length_spheroid, 'spheroid' : where[1]})
+                elif geom_3d and SpatialBackend.postgis:
+                    # Use 3D variants of perimeter and length routines on PostGIS.
+                    if perimeter:
+                        procedure_args.update({'function' : SpatialBackend.perimeter3d})
+                    elif length:
+                        procedure_args.update({'function' : SpatialBackend.length3d})
 
         # Setting up the settings for `_spatial_attribute`.
         s = {'select_field' : DistanceField(dist_att),

django/contrib/gis/db/models/sql/aggregates.py

@@ -11,6 +11,9 @@ geo_template = '%(function)s(%(field)s)'
 def convert_extent(box):
     raise NotImplementedError('Aggregate extent not implemented for this spatial backend.')
 
+def convert_extent3d(box):
+    raise NotImplementedError('Aggregate 3D extent not implemented for this spatial backend.')
+
 def convert_geom(wkt, geo_field):
     raise NotImplementedError('Aggregate method not implemented for this spatial backend.')
 
@@ -23,6 +26,14 @@ if SpatialBackend.postgis:
         xmax, ymax = map(float, ur.split())
         return (xmin, ymin, xmax, ymax)
 
+    def convert_extent3d(box3d):
+        # Box text will be something like "BOX3D(-90.0 30.0 1, -85.0 40.0 2)";
+        # parsing out and returning as a 4-tuple.
+        ll, ur = box3d[6:-1].split(',')
+        xmin, ymin, zmin = map(float, ll.split())
+        xmax, ymax, zmax = map(float, ur.split())
+        return (xmin, ymin, zmin, xmax, ymax, zmax)
+
     def convert_geom(hex, geo_field):
         if hex: return SpatialBackend.Geometry(hex)
         else: return None
@@ -94,7 +105,7 @@ class Collect(GeoAggregate):
     sql_function = SpatialBackend.collect
 
 class Extent(GeoAggregate):
-    is_extent = True
+    is_extent = '2D'
     sql_function = SpatialBackend.extent
 
 if SpatialBackend.oracle:
@@ -102,6 +113,10 @@ if SpatialBackend.oracle:
     Extent.conversion_class = GeomField
     Extent.sql_template = '%(function)s(%(field)s)'
 
+class Extent3D(GeoAggregate):
+    is_extent = '3D'
+    sql_function = SpatialBackend.extent3d
+
 class MakeLine(GeoAggregate):
     conversion_class = GeomField
     sql_function = SpatialBackend.make_line

django/contrib/gis/db/models/sql/query.py

@@ -262,7 +262,10 @@ class GeoQuery(sql.Query):
         """
         if isinstance(aggregate, self.aggregates_module.GeoAggregate):
             if aggregate.is_extent:
-                return self.aggregates_module.convert_extent(value)
+                if aggregate.is_extent == '3D':
+                    return self.aggregates_module.convert_extent3d(value)
+                else:
+                    return self.aggregates_module.convert_extent(value)
             else:
                 return self.aggregates_module.convert_geom(value, aggregate.source)
         else:

django/contrib/gis/geos/geometry.py

@@ -373,9 +373,10 @@ class GEOSGeometry(GEOSBase, ListMixin):
     @property
     def hex(self):
         """
-        Returns the HEX of the Geometry -- please note that the SRID is not
-        included in this representation, because it is not a part of the
-        OGC specification (use the `hexewkb` property instead).
+        Returns the WKB of this Geometry in hexadecimal form.  Please note
+        that the SRID and Z values are not included in this representation
+        because it is not a part of the OGC specification (use the `hexewkb` 
+        property instead).
         """
         # A possible faster, all-python, implementation:
         #  str(self.wkb).encode('hex')
@@ -384,9 +385,9 @@ class GEOSGeometry(GEOSBase, ListMixin):
     @property
     def hexewkb(self):
         """
-        Returns the HEXEWKB of this Geometry.  This is an extension of the WKB
-        specification that includes SRID and Z values taht are a part of this
-        geometry.
+        Returns the EWKB of this Geometry in hexadecimal form.  This is an 
+        extension of the WKB specification that includes SRID and Z values 
+        that are a part of this geometry.
         """
         if self.hasz:
             if not GEOS_PREPARE:

django/contrib/gis/tests/__init__.py

@@ -9,9 +9,10 @@ def geo_suite():
     some backends).
     """
     from django.conf import settings
+    from django.contrib.gis.geos import GEOS_PREPARE
     from django.contrib.gis.gdal import HAS_GDAL
     from django.contrib.gis.utils import HAS_GEOIP
-    from django.contrib.gis.tests.utils import mysql
+    from django.contrib.gis.tests.utils import postgis, mysql
 
     # The test suite.
     s = unittest.TestSuite()
@@ -32,6 +33,10 @@ def geo_suite():
     if not mysql:
         test_apps.append('distapp')
 
+    # Only PostGIS using GEOS 3.1+ can support 3D so far.
+    if postgis and GEOS_PREPARE:
+        test_apps.append('geo3d')
+
     if HAS_GDAL:
         # These tests require GDAL.
         test_suite_names.extend(['test_spatialrefsys', 'test_geoforms'])

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

+from django.contrib.gis.db import models
+
+class City3D(models.Model):
+    name = models.CharField(max_length=30)
+    point = models.PointField(dim=3)
+    objects = models.GeoManager()
+
+    def __unicode__(self):
+        return self.name
+
+class Interstate2D(models.Model):
+    name = models.CharField(max_length=30)
+    line = models.LineStringField(srid=4269)
+    objects = models.GeoManager()
+
+    def __unicode__(self):
+        return self.name
+
+class Interstate3D(models.Model):
+    name = models.CharField(max_length=30)
+    line = models.LineStringField(dim=3, srid=4269)
+    objects = models.GeoManager()
+
+    def __unicode__(self):
+        return self.name
+
+class InterstateProj2D(models.Model):
+    name = models.CharField(max_length=30)
+    line = models.LineStringField(srid=32140)
+    objects = models.GeoManager()
+
+    def __unicode__(self):
+        return self.name
+
+class InterstateProj3D(models.Model):
+    name = models.CharField(max_length=30)
+    line = models.LineStringField(dim=3, srid=32140)
+    objects = models.GeoManager()
+
+    def __unicode__(self):
+        return self.name
+
+class Polygon2D(models.Model):
+    name = models.CharField(max_length=30)
+    poly = models.PolygonField(srid=32140)
+    objects = models.GeoManager()
+    
+    def __unicode__(self):
+        return self.name
+
+class Polygon3D(models.Model):
+    name = models.CharField(max_length=30)
+    poly = models.PolygonField(dim=3, srid=32140)
+    objects = models.GeoManager()
+    
+    def __unicode__(self):
+        return self.name
+
+class Point2D(models.Model):
+    point = models.PointField()
+    objects = models.GeoManager()
+
+class Point3D(models.Model):
+    point = models.PointField(dim=3)
+    objects = models.GeoManager()
+
+class MultiPoint3D(models.Model):
+    mpoint = models.MultiPointField(dim=3)
+    objects = models.GeoManager()

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

+# Create your views here.

django/contrib/gis/utils/__init__.py

@@ -10,7 +10,7 @@ if HAS_GDAL:
     try:
         # LayerMapping requires DJANGO_SETTINGS_MODULE to be set, 
         # so this needs to be in try/except.
-        from django.contrib.gis.utils.layermapping import LayerMapping
+        from django.contrib.gis.utils.layermapping import LayerMapping, LayerMapError
     except:
         pass
     

django/contrib/gis/utils/layermapping.py

@@ -133,6 +133,9 @@ class LayerMapping(object):
     MULTI_TYPES = {1 : OGRGeomType('MultiPoint'),
                    2 : OGRGeomType('MultiLineString'),
                    3 : OGRGeomType('MultiPolygon'),
+                   OGRGeomType('Point25D').num : OGRGeomType('MultiPoint25D'),
+                   OGRGeomType('LineString25D').num : OGRGeomType('MultiLineString25D'),
+                   OGRGeomType('Polygon25D').num : OGRGeomType('MultiPolygon25D'),
                    }
 
     # Acceptable Django field types and corresponding acceptable OGR
@@ -282,19 +285,28 @@ class LayerMapping(object):
                 if self.geom_field:
                     raise LayerMapError('LayerMapping does not support more than one GeometryField per model.')
 
+                # Getting the coordinate dimension of the geometry field.
+                coord_dim = model_field.dim
+
                 try:
-                    gtype = OGRGeomType(ogr_name)
+                    if coord_dim == 3:
+                        gtype = OGRGeomType(ogr_name + '25D')
+                    else:
+                        gtype = OGRGeomType(ogr_name)
                 except OGRException:
                     raise LayerMapError('Invalid mapping for GeometryField "%s".' % field_name)
 
                 # Making sure that the OGR Layer's Geometry is compatible.
                 ltype = self.layer.geom_type
-                if not (gtype == ltype or self.make_multi(ltype, model_field)):
-                    raise LayerMapError('Invalid mapping geometry; model has %s, feature has %s.' % (fld_name, gtype))
+                if not (ltype.name.startswith(gtype.name) or self.make_multi(ltype, model_field)):
+                    raise LayerMapError('Invalid mapping geometry; model has %s%s, layer is %s.' % 
+                                        (fld_name, (coord_dim == 3 and '(dim=3)') or '', ltype))
 
                 # Setting the `geom_field` attribute w/the name of the model field
-                # that is a Geometry.
+                # that is a Geometry.  Also setting the coordinate dimension
+                # attribute.
                 self.geom_field = field_name
+                self.coord_dim = coord_dim
                 fields_val = model_field
             elif isinstance(model_field, models.ForeignKey):
                 if isinstance(ogr_name, dict):
@@ -482,6 +494,10 @@ class LayerMapping(object):
         if necessary (for example if the model field is MultiPolygonField while
         the mapped shapefile only contains Polygons).
         """
+        # Downgrade a 3D geom to a 2D one, if necessary.
+        if self.coord_dim != geom.coord_dim:
+            geom.coord_dim = self.coord_dim
+
         if self.make_multi(geom.geom_type, model_field):
             # Constructing a multi-geometry type to contain the single geometry
             multi_type = self.MULTI_TYPES[geom.geom_type.num]