[gdal-dev] Difficulty finding all parts of shapefile polygon using Python and osgeo/ogr

Jan Heckman jan.heckman at gmail.com
Wed May 21 12:44:20 PDT 2014


Hi Bruce,
Fine, I liked looking into it as a way of refamiliarizing myself, so no
time wasted as far as I'm concerned.
Something that stuck in my mind is that the shapefile looks like a
contour-plot and that the cut-off of the little shape0 might not be
desired, at least not with the sort of minimal breach between the parts we
see in your example: the split could be regarded as an artifact.
If this is broadly correct in your case and you see the contouring as an
issue, maybe we could look into steering the contouring algorithm using
(speaking loosely) a minimal breach parameter. I've done something similar,
but it's a bit of time ago.
The problem I'm viewing occurs with narrow ridge- (or canyon-)like
structures, where, in view of the discretization, you could imagine a gap
in the ridge or a nice continuous case.
(The difference would be quite noticeable when you try to climb along the
ridge [?]). This is known in contouring as an issue of interpretation and
cannot be decided analytically without additional information, leading to a
case for supplying an additional input parameter to influence the result.
Best regards,
Jan


On Tue, May 20, 2014 at 5:54 AM, Bruce Raup <bruceraup at gmail.com> wrote:

> Jan,
>
> Thank you very much for spending time on my problem.  I ended up figuring
> out that the shapefile's structure consisted of a MULTIPOLYGON that
> contained two POLYGONs, one simple one (the small one you mentioned) and
> another with "holes".  I modified my code to account for the possibility of
> this case and it works great now.
>
> Thanks again, and best regards,
> Bruce Raup
>
>
> On Sun, May 18, 2014 at 6:19 AM, Jan Heckman <jan.heckman at gmail.com>wrote:
>
>> Bruce,
>> Attached is the digging - in C++, but should help for python as well.
>> I  used C++ for easier debugging, because there just might have been a
>> bug - but there wasn't.
>> Contains c++ source for reading the geometry of the testfile, the output
>> the program generates and the original zip of the nanga test shapefile.
>> Complexity of the geometric model follows from OGC specs.
>> But then, shapefiles support e.g. a collection of geometries (to be used
>> with a single set of attributes), as well.
>> (Nanga? Nanga Parbat? Couldn't make it out[?]).
>>  Jan
>>
>>
>> On Sun, May 18, 2014 at 12:53 AM, Jan Heckman <jan.heckman at gmail.com>wrote:
>>
>>> 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
>>>>
>>>
>>>
>>
>
>
> --
> 720-383-4668
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140521/e6435520/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 330.gif
Type: image/gif
Size: 96 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140521/e6435520/attachment-0001.gif>


More information about the gdal-dev mailing list