[Community] Uh...I am unable to extract interior polygon coordinates?

[Ari Simmons]

Hi all,

I am new to Shapely (but enthusiastic about it), and recently I've
discovered a bit of a road bump.

I have a polygon shapefile that I am reading in via Fiona. This shapefile
contains BOTH polygon and multipolygon items and I need to build an array
for each feature of all the coordinates within it (i.e. both exterior
and/or interior). Notably, two of the polygon items have interior rings
(and they are valid).

I seem to have no problem accessing the exterior coordinates of the
polygon(s)/multipolygon(s) ... but I am not pulling anything for the
interior coordinates.

Do I need to take a new approach here (i.e. LinearRings)...?

Here is my function if it helps...I can send my test file (it is
small)...if you want to see what I mean...

import fiona
from numpy import asarray
from shapely.geometry import shape, Polygon, MultiPolygon

def convert_polygons(inFile):

for polys in fiona.open(inFile):
myShape = shape(polys['geometry'])
exterior_poly = 0
interior_poly = 0
if isinstance(myShape, Polygon):
print "yes, I am a polygon"
# count how many points for each interior polygon
try:
interior_poly += len(myShape.interior.coords)
except:
pass
# count how many points for each exterior polygon
exterior_poly += len(myShape.exterior.coords)
geomArray = asarray(myShape.exterior)
print geomArray
print "number of interior points in polygon " +
str(interior_poly)
print "number of exterior points in polygon " +
str(exterior_poly)
elif isinstance(myShape, MultiPolygon):
print "yes, I am a MultiPolygon"
# count how many points for each interior polygon
try:
interior_poly += len(myShape.interior.coords)
except:
pass
try:
# count how many points for each exterior polygon
exterior_poly += len(myShape.exterior.coords)
except:
pass
try:
geomArray = asarray(myShape.interior)
except:
pass
try:
geomArray = asarray(myShape.exterior)
except:
pass
print geomArray
print "number of interior points in polygon " +
str(interior_poly)
print "number of exterior points in polygon " +
str(exterior_poly)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.gispython.org/pipermail/community/attachments/20140217/6e6be42e/attachment.htm>

[Oliver Tonnhofer]

Hi,

I am new to Shapely (but enthusiastic about it), and recently I've discovered a bit of a road bump.
I have a polygon shapefile that I am reading in via Fiona. This shapefile contains BOTH polygon and multipolygon items and I need to build an array for each feature of all the coordinates within it (i.e. both exterior and/or interior). Notably, two of the polygon items have interior rings (and they are valid).
I seem to have no problem accessing the exterior coordinates of the polygon(s)/multipolygon(s) ... but I am not pulling anything for the interior coordinates.

There are errors in your code and the program tells you that, but you aren't listening :)

try:
make_error()
except:
pass

This try/except block swallows all errors. This is bad coding practice and you should avoid it.

interior_poly += len(myShape.interior.coords)

interior is a list, because a polygon can contain multiple interiors.

print "yes, I am a MultiPolygon"
# count how many points for each interior polygon
interior_poly += len(myShape.interior.coords)

and a polygon contains multiple polygons. You can iterate over myShape to get the sub-polygons.



Gru?/Regards,
Oliver

[Mike Toews]

Hi Ari,

As a follow-up, I posted an answer here:
http://stackoverflow.com/a/21922058/327026

Interior and exterior rings are structured differently. For any
polygon, there is always 1 exterior ring with zero or more interior
rings.

So looking at the structure of a geometry, `exterior` is a
`LinearRing` object, and `interiors` is a *list* of zero or more
`LinearRing` objects. Any `LinearRing` object will have `coords`,
which you can slice to see a list of the coordinates with `coords[:]`.

The following is a function that returns a dict of lists of exterior
and interior coordinates:

def extract_poly_coords(geom):
if geom.type == 'Polygon':
exterior_coords = geom.exterior.coords[:]
interior_coords = []
for int in geom.interiors:
interior_coords += i.coords[:]
elif geom.type == 'MultiPolygon':
exterior_coords = []
interior_coords = []
for part in geom:
epc = extract_poly_coords(part) # Recursive call
exterior_coords += epc['exterior_coords']
interior_coords += epc['interior_coords']
else:
raise ValueError('Unhandled geometry type: ' + repr(geom.type))
return {'exterior_coords': exterior_coords,
'interior_coords': interior_coords}

E.g.:

extract_poly_coords(myShape)

-Mike