<div dir="ltr">Hi Bruce,<div>Fine, I liked looking into it as a way of refamiliarizing myself, so no time wasted as far as I'm concerned.</div><div>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.</div>
<div>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.</div>
<div>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.</div><div>(The difference would be quite noticeable when you try to climb along the ridge <img src="cid:330@goomoji.gmail" goomoji="330" style="margin: 0px 0.2ex; vertical-align: middle;">). 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.</div>
<div>Best regards,</div><div>Jan</div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, May 20, 2014 at 5:54 AM, Bruce Raup <span dir="ltr"><<a href="mailto:bruceraup@gmail.com" target="_blank">bruceraup@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Jan,<br><br>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.<br>
<br></div>Thanks again, and best regards,<br>Bruce Raup<br></div><div class="gmail_extra"><div><div class="h5"><br><br><div class="gmail_quote">On Sun, May 18, 2014 at 6:19 AM, Jan Heckman <span dir="ltr"><<a href="mailto:jan.heckman@gmail.com" target="_blank">jan.heckman@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div>Bruce,</div>Attached is the digging - in C++, but should help for python as well.<div>I used C++ for easier debugging, because there just might have been a bug - but there wasn't.<br>
<div>Contains c++ source for reading the geometry of the testfile, the output the program generates and the original zip of the nanga test shapefile.</div>
<div>Complexity of the geometric model follows from OGC specs.</div><div>But then, shapefiles support e.g. a collection of geometries (to be used with a single set of attributes), as well.</div><div>(Nanga? Nanga Parbat? Couldn't make it out<img src="cid:330@goomoji.gmail" goomoji="330" style="margin:0px 0.2ex;vertical-align:middle">).</div>
<span><font color="#888888">
<div>Jan</div></font></span></div></div><div><div><div class="gmail_extra"><br><br><div class="gmail_quote">On Sun, May 18, 2014 at 12:53 AM, Jan Heckman <span dir="ltr"><<a href="mailto:jan.heckman@gmail.com" target="_blank">jan.heckman@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hello Bruce,<div>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.</div>
<div>See <a href="http://www.gdal.org/ogr/ogr_apitut.html" target="_blank">ogr_apitutorial</a> - all assuming you are using 1.11: <i><span style="color:rgb(0,0,0);font-family:'Lucida Grande',Verdana,Geneva,Arial,sans-serif;font-size:12px">Starting with OGR 1.11, </span><a href="http://trac.osgeo.org/gdal/wiki/rfc41_multiple_geometry_fields" style="color:rgb(70,101,162);text-decoration:none;font-family:'Lucida Grande',Verdana,Geneva,Arial,sans-serif;font-size:12px" target="_blank">several geometry fields</a></i><span style="color:rgb(0,0,0);font-family:'Lucida Grande',Verdana,Geneva,Arial,sans-serif;font-size:12px"><i> can be associated to a feature</i>.</span></div>
<div>The little hidden piece starts at (approx) 467767;3900845. It happens to be the first part of the shape.</div><div>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..</div>
<div>Good luck,</div><div>Jan</div></div><div class="gmail_extra"><br><br><div class="gmail_quote"><div><div>On Thu, May 15, 2014 at 6:13 PM, Bruce Raup <span dir="ltr"><<a href="mailto:bruceraup@gmail.com" target="_blank">bruceraup@gmail.com</a>></span> wrote:<br>
</div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><div><div dir="ltr"><div>Hello,<br><br></div>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:<br>
<br>$ shpdump simple_shapes_with_holes.shp <br>Shapefile Type: Polygon # of Shapes: 3<br><br>File Bounds: ( -1.113, 0.146,0,0)<br> to ( 0.751, 1.101,0,0)<br><br>Shape:0 (Polygon) nVertices=11, nParts=2<br>
Bounds:( -1.113, 0.146, 0, 0)<br> to ( -0.493, 0.723, 0, 0)<br> ( -1.087, 0.151, 0, 0) Ring <br> ( -1.113, 0.686, 0, 0) <br> ( -1.109, 0.712, 0, 0) <br>
( -0.506, 0.723, 0, 0) <br> ( -0.493, 0.146, 0, 0) <br> ( -1.087, 0.151, 0, 0) <br> + ( -0.896, 0.346, 0, 0) Ring <br> ( -0.682, 0.339, 0, 0) <br>
( -0.684, 0.556, 0, 0) <br> ( -0.909, 0.558, 0, 0) <br> ( -0.896, 0.346, 0, 0) <br><br>Shape:1 (Polygon) nVertices=19, nParts=4<br> Bounds:( -0.127, 0.159, 0, 0)<br>
to ( 0.751, 0.783, 0, 0)<br> ( 0.402, 0.783, 0, 0) Ring <br> ( 0.751, 0.159, 0, 0) <br> ( -0.127, 0.198, 0, 0) <br> ( 0.402, 0.783, 0, 0) <br>
+ ( 0.298, 0.539, 0, 0) Ring <br> ( 0.417, 0.523, 0, 0) <br> ( 0.426, 0.612, 0, 0) <br> ( 0.339, 0.625, 0, 0) <br> ( 0.298, 0.539, 0, 0) <br>
+ ( 0.083, 0.263, 0, 0) Ring <br> ( 0.248, 0.252, 0, 0) <br> ( 0.213, 0.396, 0, 0) <br> ( 0.107, 0.367, 0, 0) <br> ( 0.083, 0.263, 0, 0) <br>
+ ( 0.361, 0.281, 0, 0) Ring <br> ( 0.530, 0.265, 0, 0) <br> ( 0.435, 0.413, 0, 0) <br> ( 0.330, 0.435, 0, 0) <br> ( 0.361, 0.281, 0, 0) <br>
<br>Shape:2 (Polygon) nVertices=6, nParts=1<br> Bounds:( -0.389, 0.681, 0, 0)<br> to ( 0.141, 1.101, 0, 0)<br> ( -0.383, 1.005, 0, 0) Ring <br> ( -0.136, 1.101, 0, 0) <br>
( 0.141, 0.964, 0, 0) <br> ( -0.154, 0.681, 0, 0) <br> ( -0.389, 0.809, 0, 0) <br> ( -0.383, 1.005, 0, 0) <br><br clear="all"><div><div><br><br><br>$ shpdump simple_shapes_with_holes_flattened.shp <br>
Shapefile Type: Polygon # of Shapes: 7<br><br>File Bounds: ( -1.113, 0.146,0,0)<br> to ( 0.751, 1.101,0,0)<br><br>Shape:0 (Polygon) nVertices=6, nParts=1<br> Bounds:( -1.113, 0.146, 0, 0)<br>
to ( -0.493, 0.723, 0, 0)<br> ( -1.087, 0.151, 0, 0) Ring <br> ( -1.113, 0.686, 0, 0) <br> ( -1.109, 0.712, 0, 0) <br> ( -0.506, 0.723, 0, 0) <br>
( -0.493, 0.146, 0, 0) <br> ( -1.087, 0.151, 0, 0) <br><br>Shape:1 (Polygon) nVertices=5, nParts=1<br> Bounds:( -0.909, 0.339, 0, 0)<br> to ( -0.682, 0.558, 0, 0)<br>
( -0.896, 0.346, 0, 0) Ring <br> ( -0.909, 0.558, 0, 0) <br> ( -0.684, 0.556, 0, 0) <br> ( -0.682, 0.339, 0, 0) <br> ( -0.896, 0.346, 0, 0) <br>
<br>Shape:2 (Polygon) nVertices=4, nParts=1<br> Bounds:( -0.127, 0.159, 0, 0)<br> to ( 0.751, 0.783, 0, 0)<br> ( 0.402, 0.783, 0, 0) Ring <br> ( 0.751, 0.159, 0, 0) <br>
( -0.127, 0.198, 0, 0) <br> ( 0.402, 0.783, 0, 0) <br><br>Shape:3 (Polygon) nVertices=5, nParts=1<br> Bounds:( 0.298, 0.523, 0, 0)<br> to ( 0.426, 0.625, 0, 0)<br>
( 0.298, 0.539, 0, 0) Ring <br> ( 0.339, 0.625, 0, 0) <br> ( 0.426, 0.612, 0, 0) <br> ( 0.417, 0.523, 0, 0) <br> ( 0.298, 0.539, 0, 0) <br>
<br>Shape:4 (Polygon) nVertices=5, nParts=1<br> Bounds:( 0.083, 0.252, 0, 0)<br> to ( 0.248, 0.396, 0, 0)<br> ( 0.083, 0.263, 0, 0) Ring <br> ( 0.107, 0.367, 0, 0) <br>
( 0.213, 0.396, 0, 0) <br> ( 0.248, 0.252, 0, 0) <br> ( 0.083, 0.263, 0, 0) <br><br>Shape:5 (Polygon) nVertices=5, nParts=1<br> Bounds:( 0.330, 0.265, 0, 0)<br>
to ( 0.530, 0.435, 0, 0)<br> ( 0.361, 0.281, 0, 0) Ring <br> ( 0.330, 0.435, 0, 0) <br> ( 0.435, 0.413, 0, 0) <br> ( 0.530, 0.265, 0, 0) <br>
( 0.361, 0.281, 0, 0) <br><br>Shape:6 (Polygon) nVertices=6, nParts=1<br> Bounds:( -0.389, 0.681, 0, 0)<br> to ( 0.141, 1.101, 0, 0)<br> ( -0.383, 1.005, 0, 0) Ring <br>
( -0.136, 1.101, 0, 0) <br> ( 0.141, 0.964, 0, 0) <br> ( -0.154, 0.681, 0, 0) <br> ( -0.389, 0.809, 0, 0) <br> ( -0.383, 1.005, 0, 0) <br>
<br>
<br></div><div>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<br>
<br><br>Shapefile Type: Polygon # of Shapes: 1<br><br>File Bounds: ( 462611.523, 3894164.223,0,0)<br> to ( 470338.347, 3901974.217,0,0)<br><br>Shape:0 (Polygon) nVertices=1575, nParts=14<br> Bounds:( 462611.523, 3894164.223, 0, 0)<br>
to ( 470338.347, 3901974.217, 0, 0)<br> ( 467768.987, 3900844.546, 0, 0) Ring <br> ( 467744.058, 3900868.748, 0, 0) <br> ( 467702.110, 3900879.841, 0, 0) <br> ( 467670.108, 3900921.799, 0, 0) <br>
( 467632.530, 3900960.784, 0, 0) <br> ( 467689.198, 3900957.167, 0, 0) <br> ( 467751.895, 3900922.201, 0, 0) <br> ( 467785.655, 3900861.426, 0, 0) <br> ( 467768.987, 3900844.546, 0, 0) <br>
+ ( 467774.970, 3900834.757, 0, 0) Ring <br>
( 467796.445, 3900858.379, 0, 0) <br> ( 467850.763, 3900846.242, 0, 0) <br> ( 467900.197, 3900841.419, 0, 0) <br> ( 467927.928, 3900814.893, 0, 0) <br> ( 467936.368, 3900767.871, 0, 0) <br>
.<br>
.<br>.<br><br></div><div>But my python code doesn't see all the parts. <br><br><br>$ python<br>Python 2.7 (r27:82500, Aug 07 2010, 16:54:59) [GCC] on linux2<br>Type "help", "copyright", "credits" or "license" for more information.<br>
>>> from osgeo import ogr<br>>>> sf = ogr.Open('segments.shp')<br>>>> sl = sf.GetLayer(0)<br>>>> sl<br><osgeo.ogr.Layer; proxy of <Swig Object of type 'OGRLayerShadow *' at 0xb7472980> ><br>
>>> feat = sl[0]<br>>>> feat<br><osgeo.ogr.Feature; proxy of <Swig Object of type 'OGRFeatureShadow *' at 0xb7472620> ><br>>>> sg = feat.GetGeometryRef()<br>>>> sg<br>
<osgeo.ogr.Geometry; proxy of <Swig Object of type 'OGRGeometryShadow *' at 0xb749e770> ><br>>>> sg.GetGeometryCount()<br>2<br>>>> <br><br><br></div><div>I would expect the output of sg.GetGeometryCount() to be 15. Am I doing something wrong?<br>
<br></div><div>The complete script is at<br><br><a href="http://cires.colorado.edu/~braup/OGR/flatten_shapefile.py" target="_blank">http://cires.colorado.edu/~braup/OGR/flatten_shapefile.py</a><br><br></div><div>The testfile is at<br>
<br>
<a href="http://cires.colorado.edu/~braup/OGR/simple_shapes_with_holes.zip" target="_blank">http://cires.colorado.edu/~braup/OGR/simple_shapes_with_holes.zip</a><br><br></div><div>The "real" testfile on which the script fails is at:<br>
<br><a href="http://cires.colorado.edu/~braup/OGR/nanga_segments.zip" target="_blank">http://cires.colorado.edu/~braup/OGR/nanga_segments.zip</a><br><br></div><div>Any pointers would be greatly appreciated.<br><br></div>
<div>Cheers,<br>Bruce Raup<span><font color="#888888"><br>
</font></span></div><span><font color="#888888"><div><br>-- <br>Bruce Raup<br><a href="http://cires.colorado.edu/~braup/" target="_blank">http://cires.colorado.edu/~braup/</a>
</div></font></span></div></div>
<br></div></div>_______________________________________________<br>
gdal-dev mailing list<br>
<a href="mailto:gdal-dev@lists.osgeo.org" target="_blank">gdal-dev@lists.osgeo.org</a><br>
<a href="http://lists.osgeo.org/mailman/listinfo/gdal-dev" target="_blank">http://lists.osgeo.org/mailman/listinfo/gdal-dev</a><br></blockquote></div><br></div>
</blockquote></div><br></div>
</div></div></blockquote></div><br><br clear="all"><br></div></div><span class="HOEnZb"><font color="#888888">-- <br><a href="tel:720-383-4668" value="+17203834668" target="_blank">720-383-4668</a>
</font></span></div>
</blockquote></div><br></div>