[gdal-dev] Difficulty finding all parts of shapefile polygon using Python and osgeo/ogr
Jan Heckman
jan.heckman at gmail.com
Sat May 17 15:53:23 PDT 2014
Hello Bruce,
The point is that the geometrycount identifies two polygons, one a small
little piece sort of inside a bay of your shape, the other the rest, along
with the holes.
See ogr_apitutorial <http://www.gdal.org/ogr/ogr_apitut.html> - all
assuming you are using 1.11: *Starting with OGR 1.11, several geometry
fields <http://trac.osgeo.org/gdal/wiki/rfc41_multiple_geometry_fields>** can
be associated to a feature*.
The little hidden piece starts at (approx) 467767;3900845. It happens to be
the first part of the shape.
The example code in the link above should put you closer to all this;
however, retrieving the polygon(part)s will take a bit more digging..
Good luck,
Jan
On Thu, May 15, 2014 at 6:13 PM, Bruce Raup <bruceraup at gmail.com> wrote:
> 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/
>
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140518/75f16837/attachment-0001.html>
More information about the gdal-dev
mailing list