[gdal-dev] Difficulty finding all parts of shapefile polygon using Python and osgeo/ogr
Bruce Raup
bruceraup at gmail.com
Thu May 15 09:13:14 PDT 2014
Hello,
I have a Python script that reads a shapefile, converts holes in polygons
to top-level polygons, and writes a new polygon. The new top-level
polygons inherit attributes from their parents, with optional modifications
of specific attributes. For a simple test shapefile, the output of shpdump
shows what I mean:
$ shpdump simple_shapes_with_holes.shp
Shapefile Type: Polygon # of Shapes: 3
File Bounds: ( -1.113, 0.146,0,0)
to ( 0.751, 1.101,0,0)
Shape:0 (Polygon) nVertices=11, nParts=2
Bounds:( -1.113, 0.146, 0, 0)
to ( -0.493, 0.723, 0, 0)
( -1.087, 0.151, 0, 0) Ring
( -1.113, 0.686, 0, 0)
( -1.109, 0.712, 0, 0)
( -0.506, 0.723, 0, 0)
( -0.493, 0.146, 0, 0)
( -1.087, 0.151, 0, 0)
+ ( -0.896, 0.346, 0, 0) Ring
( -0.682, 0.339, 0, 0)
( -0.684, 0.556, 0, 0)
( -0.909, 0.558, 0, 0)
( -0.896, 0.346, 0, 0)
Shape:1 (Polygon) nVertices=19, nParts=4
Bounds:( -0.127, 0.159, 0, 0)
to ( 0.751, 0.783, 0, 0)
( 0.402, 0.783, 0, 0) Ring
( 0.751, 0.159, 0, 0)
( -0.127, 0.198, 0, 0)
( 0.402, 0.783, 0, 0)
+ ( 0.298, 0.539, 0, 0) Ring
( 0.417, 0.523, 0, 0)
( 0.426, 0.612, 0, 0)
( 0.339, 0.625, 0, 0)
( 0.298, 0.539, 0, 0)
+ ( 0.083, 0.263, 0, 0) Ring
( 0.248, 0.252, 0, 0)
( 0.213, 0.396, 0, 0)
( 0.107, 0.367, 0, 0)
( 0.083, 0.263, 0, 0)
+ ( 0.361, 0.281, 0, 0) Ring
( 0.530, 0.265, 0, 0)
( 0.435, 0.413, 0, 0)
( 0.330, 0.435, 0, 0)
( 0.361, 0.281, 0, 0)
Shape:2 (Polygon) nVertices=6, nParts=1
Bounds:( -0.389, 0.681, 0, 0)
to ( 0.141, 1.101, 0, 0)
( -0.383, 1.005, 0, 0) Ring
( -0.136, 1.101, 0, 0)
( 0.141, 0.964, 0, 0)
( -0.154, 0.681, 0, 0)
( -0.389, 0.809, 0, 0)
( -0.383, 1.005, 0, 0)
$ shpdump simple_shapes_with_holes_flattened.shp
Shapefile Type: Polygon # of Shapes: 7
File Bounds: ( -1.113, 0.146,0,0)
to ( 0.751, 1.101,0,0)
Shape:0 (Polygon) nVertices=6, nParts=1
Bounds:( -1.113, 0.146, 0, 0)
to ( -0.493, 0.723, 0, 0)
( -1.087, 0.151, 0, 0) Ring
( -1.113, 0.686, 0, 0)
( -1.109, 0.712, 0, 0)
( -0.506, 0.723, 0, 0)
( -0.493, 0.146, 0, 0)
( -1.087, 0.151, 0, 0)
Shape:1 (Polygon) nVertices=5, nParts=1
Bounds:( -0.909, 0.339, 0, 0)
to ( -0.682, 0.558, 0, 0)
( -0.896, 0.346, 0, 0) Ring
( -0.909, 0.558, 0, 0)
( -0.684, 0.556, 0, 0)
( -0.682, 0.339, 0, 0)
( -0.896, 0.346, 0, 0)
Shape:2 (Polygon) nVertices=4, nParts=1
Bounds:( -0.127, 0.159, 0, 0)
to ( 0.751, 0.783, 0, 0)
( 0.402, 0.783, 0, 0) Ring
( 0.751, 0.159, 0, 0)
( -0.127, 0.198, 0, 0)
( 0.402, 0.783, 0, 0)
Shape:3 (Polygon) nVertices=5, nParts=1
Bounds:( 0.298, 0.523, 0, 0)
to ( 0.426, 0.625, 0, 0)
( 0.298, 0.539, 0, 0) Ring
( 0.339, 0.625, 0, 0)
( 0.426, 0.612, 0, 0)
( 0.417, 0.523, 0, 0)
( 0.298, 0.539, 0, 0)
Shape:4 (Polygon) nVertices=5, nParts=1
Bounds:( 0.083, 0.252, 0, 0)
to ( 0.248, 0.396, 0, 0)
( 0.083, 0.263, 0, 0) Ring
( 0.107, 0.367, 0, 0)
( 0.213, 0.396, 0, 0)
( 0.248, 0.252, 0, 0)
( 0.083, 0.263, 0, 0)
Shape:5 (Polygon) nVertices=5, nParts=1
Bounds:( 0.330, 0.265, 0, 0)
to ( 0.530, 0.435, 0, 0)
( 0.361, 0.281, 0, 0) Ring
( 0.330, 0.435, 0, 0)
( 0.435, 0.413, 0, 0)
( 0.530, 0.265, 0, 0)
( 0.361, 0.281, 0, 0)
Shape:6 (Polygon) nVertices=6, nParts=1
Bounds:( -0.389, 0.681, 0, 0)
to ( 0.141, 1.101, 0, 0)
( -0.383, 1.005, 0, 0) Ring
( -0.136, 1.101, 0, 0)
( 0.141, 0.964, 0, 0)
( -0.154, 0.681, 0, 0)
( -0.389, 0.809, 0, 0)
( -0.383, 1.005, 0, 0)
I have a shapefile for which this doesn't work, however, and I suspect it's
perhaps a difference in the way the inner rings are represented. In this
real shapefile, the top part of the output from shpdump looks like
Shapefile Type: Polygon # of Shapes: 1
File Bounds: ( 462611.523, 3894164.223,0,0)
to ( 470338.347, 3901974.217,0,0)
Shape:0 (Polygon) nVertices=1575, nParts=14
Bounds:( 462611.523, 3894164.223, 0, 0)
to ( 470338.347, 3901974.217, 0, 0)
( 467768.987, 3900844.546, 0, 0) Ring
( 467744.058, 3900868.748, 0, 0)
( 467702.110, 3900879.841, 0, 0)
( 467670.108, 3900921.799, 0, 0)
( 467632.530, 3900960.784, 0, 0)
( 467689.198, 3900957.167, 0, 0)
( 467751.895, 3900922.201, 0, 0)
( 467785.655, 3900861.426, 0, 0)
( 467768.987, 3900844.546, 0, 0)
+ ( 467774.970, 3900834.757, 0, 0) Ring
( 467796.445, 3900858.379, 0, 0)
( 467850.763, 3900846.242, 0, 0)
( 467900.197, 3900841.419, 0, 0)
( 467927.928, 3900814.893, 0, 0)
( 467936.368, 3900767.871, 0, 0)
.
.
.
But my python code doesn't see all the parts.
$ python
Python 2.7 (r27:82500, Aug 07 2010, 16:54:59) [GCC] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from osgeo import ogr
>>> sf = ogr.Open('segments.shp')
>>> sl = sf.GetLayer(0)
>>> sl
<osgeo.ogr.Layer; proxy of <Swig Object of type 'OGRLayerShadow *' at
0xb7472980> >
>>> feat = sl[0]
>>> feat
<osgeo.ogr.Feature; proxy of <Swig Object of type 'OGRFeatureShadow *' at
0xb7472620> >
>>> sg = feat.GetGeometryRef()
>>> sg
<osgeo.ogr.Geometry; proxy of <Swig Object of type 'OGRGeometryShadow *' at
0xb749e770> >
>>> sg.GetGeometryCount()
2
>>>
I would expect the output of sg.GetGeometryCount() to be 15. Am I doing
something wrong?
The complete script is at
http://cires.colorado.edu/~braup/OGR/flatten_shapefile.py
The testfile is at
http://cires.colorado.edu/~braup/OGR/simple_shapes_with_holes.zip
The "real" testfile on which the script fails is at:
http://cires.colorado.edu/~braup/OGR/nanga_segments.zip
Any pointers would be greatly appreciated.
Cheers,
Bruce Raup
--
Bruce Raup
http://cires.colorado.edu/~braup/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140515/43a657cc/attachment.html>
More information about the gdal-dev
mailing list