> On 26 October 2014 00:18, Hartmut Kaiser <hartmut.kai...@gmail.com> wrote:
> At this point we assume, that polys[0] is a linear ring to be interpreted
> as
> a polygon exterior and polys[1:] are the corresponding interiors for
> polys[0].
>
> Here are our questions:
>
> Is this assumption correct?
> Is there any detailed documentation describing the structure of the
> returned
> geometries?
> Are the linear rings supposed to be correctly oriented (area > 0 for
> exteriors and area < 0 for the interiors)?
>
> Hello Hartmut,
> In brief, the answers are no, no and yes.
> In more detail, assuming polys is not empty then it will contain one or
> more polygon exteriors and zero or more interiors, and they can be in any
> order. Here is a simple example where polys[0] is an interior and
> polys[1] an exterior:
>
> x = [0, 0, 1, 1, 0.5]
> y = [0, 1, 0, 1, 0.5]
> z = [0.5, 0.5, 0.5, 0.5, 0]
> triang = tri.Triangulation(x, y)
> contour = plt.tricontourf(triang, z, levels=[0.2, 0.4])
> The returned geometries are purposefully not documented. They are an
> 'implementation detail' and not considered part of the public interface.
> and as such they could change at any time and hence should not be relied
> upon. Of course you can choose to access them if you wish, as I do myself
> sometimes, but we make no promises about what the order of the polygons
> is, or that it won't change tomorrow.
> In reality the order of the polygons is chosen to be something that is
> easy for both the contour generation routines to create and for the
> various backends to render. If you were to look at the output generated
> by contourf, you will see that it is organised differently from that
> produced by tricontourf and is more like you would like it to be, i.e. one
> or more groups of an exterior polygon followed by zero or more
> interiors. This is historical as the contourf code dates from before all
> of the backends were able to render arbitrary groups of exterior and
> interior polygons, and so the contourf code has to calculate the order for
> the backends. When the tricontourf code was written the backends were all
> able to calculate the exterior/interior containment themselves, so there
> was no need for tricontourf to do it as well.
Thanks Ian! Your detailed answer is much appreciated.
As you might have already guessed, we have quite some problems creating clean
geometries from the generated contour data. I have tried to put together one
(reasonably) small test case illustrating at least one of our issues. I
apologize for the lengthy code attached.
The two attached figures demonstrate what we see. Matplotlib.png (generated by
the library) does not really look ok. Also, the attached shape.png shows that
there are a lot of geometries generated which are self-intersecting
(highlighted in dark blue), and we already skip polygons with less than 3
points. BTW, there are many polygons stacked with the same geometries.
Anything we can do about this?
Thanks!
Regards Hartmut
---------------
http://boost-spirit.com
http://stellar.cct.lsu.edu
import sys
import matplotlib.pyplot as plt
import matplotlib.tri as tri
from matplotlib import path
from shapely.geometry import geo, Polygon, Point
import numpy
# create contours for test case
def create_test_contours():
x = numpy.array( [-76.34244,-76.36786,-76.35223,-76.40118,-76.3905,
-76.41039,-76.42416,-76.43492,-76.4437,-76.44031,
-76.47977,-76.47267,-76.45914,-76.50938,-76.50371,
-76.50534,-76.53763,-76.53808,-76.53674,-76.57217,
-76.57292,-76.56992,-76.61163,-76.60898,-76.59344,
-76.65379,-76.65218,-76.62014] )
y = numpy.array( [37.00411,36.94697,36.91855,36.97744,36.89212,36.94417,
36.9803,36.90823,36.95995,36.99866,36.9427,36.98559,
37.02776,36.97269,37.01178,37.05462,36.99284,37.03721,
37.08237,37.02015,37.06149,37.09983,37.04387,37.08551,
37.1161,37.08103,37.11076,37.15342] )
z = numpy.ma.masked_equal(
[0.57833663768,0.59585130282,0.57725037291,0.61816105718,
-99999,0.62808432966,0.62889707758,0.63214800973,0.62361256066,
0.63571681755,0.63464181838,0.63445566205,-99999,0.63305024240,
0.63890298185,-99999,-99999,0.0014958446745,-99999,
0.0035448925787,0.00056665199692,-99999,0.0022350042433,0.0013006928908,
-99999,0.0035717841844,0.0036461797716,-99999], -99999 )
ele = [(1,0,3),(2,1,4),(4,1,5),(5,1,3),(5,3,6),(7,4,5),(7,5,8),(8,5,6),
(6,9,8),(7,8,10),(10,8,11),(8,9,11),(9,12,11),(10,11,13),(13,11,14),(11,12,14),
(14,12,15),(13,14,16),(16,14,17),(14,15,17),(17,15,18),(16,17,19),(19,17,20),(17,18,20),
(20,18,21),(19,20,22),(22,20,23),(20,21,23),(23,21,24),(22,23,25),(25,23,26),(23,24,26),
(24,27,26)]
triang = tri.Triangulation(x, y, triangles=ele)
levels = numpy.linspace(0, numpy.max(z), num=31)
contours = plt.tricontourf(triang, z, levels=levels)
plt.tricontour(triang, z, levels=levels)
return contours
# classify polygons based on their area
def classify_polygons(polys):
def signed_area(ring):
v2 = numpy.roll(ring, -1, axis=0)
return numpy.cross(ring, v2).sum() / 2.0
exteriors = []
interiors = []
for p in polys:
if signed_area(p) >= 0:
exteriors.append(p)
else:
interiors.append(p)
return exteriors, interiors
# returns mask for points which are inside given polygon
def points_inside_poly(points, polygon):
p = path.Path(polygon)
return p.contains_points(points)
# extract shapely Polygon objects from the contours
def extract_geometries(contours):
geoms = []
for coll in contours.collections:
for p in coll.get_paths():
polys = [p for p in p.to_polygons() if p.shape[0] >= 3]
exteriors, interiors = classify_polygons(polys)
if len(interiors) > 0:
interior_points = [pts[0] for pts, a in interiors]
for exterior in exteriors:
if len(interiors) > 0:
inout = points_inside_poly(interior_points, exterior)
my_interiors = [g for i, g in enumerate(interiors) if inout[i]]
poly = Polygon(exterior, my_interiors)
else:
poly = Polygon(exterior)
geoms.append(poly)
return geoms
def main(argv = None):
if argv is None:
argv = sys.argv
contours = create_test_contours()
geoms = extract_geometries(contours)
plt.show()
if __name__ == "__main__":
sys.exit(main())
------------------------------------------------------------------------------
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
