From 69535b7b1383c6d74fac4f5a1b1a89c2c9956e23 Mon Sep 17 00:00:00 2001 From: Justin Bronn Date: Sat, 17 Oct 2009 17:32:25 +0000 Subject: [PATCH] The `OGRGeometry.coord_dim` property may now be set; implemented a work-around for an OGR bug that changed geometries to 3D after transformation. Refs #11433. git-svn-id: http://code.djangoproject.com/svn/django/trunk@11628 bcc190cf-cafb-0310-a4f2-bffc1f526a37 --- django/contrib/gis/gdal/geometries.py | 53 +++++++++++++++------- django/contrib/gis/gdal/prototypes/geom.py | 3 +- django/contrib/gis/gdal/tests/test_geom.py | 12 +++++ 3 files changed, 51 insertions(+), 17 deletions(-) diff --git a/django/contrib/gis/gdal/geometries.py b/django/contrib/gis/gdal/geometries.py index 05b824b95d2..982ed5414ee 100644 --- a/django/contrib/gis/gdal/geometries.py +++ b/django/contrib/gis/gdal/geometries.py @@ -29,7 +29,7 @@ +proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs >>> print mpnt MULTIPOINT (-89.999930378602485 29.999797886557641,-89.999930378602485 29.999797886557641) - + The OGRGeomType class is to make it easy to specify an OGR geometry type: >>> from django.contrib.gis.gdal import OGRGeomType >>> gt1 = OGRGeomType(3) # Using an integer for the type @@ -78,7 +78,7 @@ class OGRGeometry(GDALBase): geom_input = buffer(a2b_hex(geom_input.upper())) str_instance = False - # Constructing the geometry, + # Constructing the geometry, if str_instance: # Checking if unicode if isinstance(geom_input, unicode): @@ -130,12 +130,12 @@ class OGRGeometry(GDALBase): self.__class__ = GEO_CLASSES[self.geom_type.num] @classmethod - def from_bbox(cls, bbox): + def from_bbox(cls, bbox): "Constructs a Polygon from a bounding box (4-tuple)." x0, y0, x1, y1 = bbox return OGRGeometry( 'POLYGON((%s %s, %s %s, %s %s, %s %s, %s %s))' % ( x0, y0, x0, y1, x1, y1, x1, y0, x0, y0) ) - + def __del__(self): "Deletes this Geometry." if self._ptr: capi.destroy_geom(self._ptr) @@ -179,10 +179,17 @@ class OGRGeometry(GDALBase): "Returns 0 for points, 1 for lines, and 2 for surfaces." return capi.get_dims(self.ptr) - @property - def coord_dim(self): + def _get_coord_dim(self): "Returns the coordinate dimension of the Geometry." - return capi.get_coord_dims(self.ptr) + return capi.get_coord_dim(self.ptr) + + def _set_coord_dim(self, dim): + "Sets the coordinate dimension of this Geometry." + if not dim in (2, 3): + raise ValueError('Geometry dimension must be either 2 or 3') + capi.set_coord_dim(self.ptr, dim) + + coord_dim = property(_get_coord_dim, _set_coord_dim) @property def geom_count(self): @@ -237,7 +244,7 @@ class OGRGeometry(GDALBase): return self.envelope.tuple #### SpatialReference-related Properties #### - + # The SRS property def _get_srs(self): "Returns the Spatial Reference for this Geometry." @@ -298,7 +305,7 @@ class OGRGeometry(GDALBase): Returns the GeoJSON representation of this Geometry (requires GDAL 1.5+). """ - if GEOJSON: + if GEOJSON: return capi.to_json(self.ptr) else: raise NotImplementedError('GeoJSON output only supported on GDAL 1.5+.') @@ -335,7 +342,7 @@ class OGRGeometry(GDALBase): def wkt(self): "Returns the WKT representation of the Geometry." return capi.to_wkt(self.ptr, byref(c_char_p())) - + #### Geometry Methods #### def clone(self): "Clones this OGR Geometry." @@ -363,6 +370,16 @@ class OGRGeometry(GDALBase): klone = self.clone() klone.transform(coord_trans) return klone + + # Have to get the coordinate dimension of the original geometry + # so it can be used to reset the transformed geometry's dimension + # afterwards. This is done because of GDAL bug (in versions prior + # to 1.7) that turns geometries 3D after transformation, see: + # http://trac.osgeo.org/gdal/changeset/17792 + orig_dim = self.coord_dim + + # Depending on the input type, use the appropriate OGR routine + # to perform the transformation. if isinstance(coord_trans, CoordTransform): capi.geom_transform(self.ptr, coord_trans.ptr) elif isinstance(coord_trans, SpatialReference): @@ -373,6 +390,10 @@ class OGRGeometry(GDALBase): else: raise TypeError('Transform only accepts CoordTransform, SpatialReference, string, and integer objects.') + # Setting with original dimension, see comment above. + if self.coord_dim != orig_dim: + self.coord_dim = orig_dim + def transform_to(self, srs): "For backwards-compatibility." self.transform(srs) @@ -391,7 +412,7 @@ class OGRGeometry(GDALBase): def intersects(self, other): "Returns True if this geometry intersects with the other." return self._topology(capi.ogr_intersects, other) - + def equals(self, other): "Returns True if this geometry is equivalent to the other." return self._topology(capi.ogr_equals, other) @@ -436,7 +457,7 @@ class OGRGeometry(GDALBase): @property def convex_hull(self): """ - Returns the smallest convex Polygon that contains all the points in + Returns the smallest convex Polygon that contains all the points in this Geometry. """ return self._geomgen(capi.geom_convex_hull) @@ -456,7 +477,7 @@ class OGRGeometry(GDALBase): return self._geomgen(capi.geom_intersection, other) def sym_difference(self, other): - """ + """ Returns a new geometry which is the symmetric difference of this geometry and the other. """ @@ -545,7 +566,7 @@ class LineString(OGRGeometry): def y(self): "Returns the Y coordinates in a list." return self._listarr(capi.gety) - + @property def z(self): "Returns the Z coordinates in a list." @@ -610,7 +631,7 @@ class GeometryCollection(OGRGeometry): raise OGRIndexError('index out of range: %s' % index) else: return OGRGeometry(capi.clone_geom(capi.get_geom_ref(self.ptr, index)), self.srs) - + def __iter__(self): "Iterates over each Geometry." for i in xrange(self.geom_count): @@ -658,5 +679,5 @@ GEO_CLASSES = {1 : Point, 5 : MultiLineString, 6 : MultiPolygon, 7 : GeometryCollection, - 101: LinearRing, + 101: LinearRing, } diff --git a/django/contrib/gis/gdal/prototypes/geom.py b/django/contrib/gis/gdal/prototypes/geom.py index 2898198d0d8..e40f33e7455 100644 --- a/django/contrib/gis/gdal/prototypes/geom.py +++ b/django/contrib/gis/gdal/prototypes/geom.py @@ -83,7 +83,8 @@ get_geom_srs = srs_output(lgdal.OGR_G_GetSpatialReference, [c_void_p]) get_area = double_output(lgdal.OGR_G_GetArea, [c_void_p]) get_centroid = void_output(lgdal.OGR_G_Centroid, [c_void_p, c_void_p]) get_dims = int_output(lgdal.OGR_G_GetDimension, [c_void_p]) -get_coord_dims = int_output(lgdal.OGR_G_GetCoordinateDimension, [c_void_p]) +get_coord_dim = int_output(lgdal.OGR_G_GetCoordinateDimension, [c_void_p]) +set_coord_dim = void_output(lgdal.OGR_G_SetCoordinateDimension, [c_void_p, c_int], errcheck=False) get_geom_count = int_output(lgdal.OGR_G_GetGeometryCount, [c_void_p]) get_geom_name = const_string_output(lgdal.OGR_G_GetGeometryName, [c_void_p]) diff --git a/django/contrib/gis/gdal/tests/test_geom.py b/django/contrib/gis/gdal/tests/test_geom.py index c920adc6c08..b5d8046c0e9 100644 --- a/django/contrib/gis/gdal/tests/test_geom.py +++ b/django/contrib/gis/gdal/tests/test_geom.py @@ -319,6 +319,18 @@ class OGRGeomTest(unittest.TestCase): self.assertAlmostEqual(trans.x, p.x, prec) self.assertAlmostEqual(trans.y, p.y, prec) + def test09c_transform_dim(self): + "Testing coordinate dimension is the same on transformed geometries." + ls_orig = OGRGeometry('LINESTRING(-104.609 38.255)', 4326) + ls_trans = OGRGeometry('LINESTRING(992385.4472045 481455.4944650)', 2774) + + prec = 3 + ls_orig.transform(ls_trans.srs) + # Making sure the coordinate dimension is still 2D. + self.assertEqual(2, ls_orig.coord_dim) + self.assertAlmostEqual(ls_trans.x[0], ls_orig.x[0], prec) + self.assertAlmostEqual(ls_trans.y[0], ls_orig.y[0], prec) + def test10_difference(self): "Testing difference()." for i in xrange(len(topology_geoms)):