On Mon, May 11, 2009 at 3:36 AM, Mario Ceresa <mario.cer...@torrescalla.it> wrote: > Hi everybody, > I'd like to remove all the pixeles from an image which are inside a > given polygon. The first idea that come into my mind was to generate > all the points which are in the bounding box and then check if they > actually are inside the polygon: > > b = poly.bounds > xl = b[2] - b[0] > yl = b[3] - b[1] > points = (Point(i[0]+b[0],i[1]+b[1]) for i in numpy.ndindex((xl,yl))) > ps = list((p.x,p.y) for p in points if poly.contains(p)) > > For this test b was (38.0, 1073.0, 679.0, 1977.0) which lead to 579464 > points to be checked. > > Actually this is quite slow so I cannot use it on larger images: > > > In [10]: %time a=list(iterops.contains(poly,points)) > CPU times: user 43.82 s, sys: 0.21 s, total: 44.03 s > Wall time: 44.10 s > > In [12]: %time a=list(iterops.disjoint(poly,points)) > CPU times: user 42.64 s, sys: 0.20 s, total: 42.83 s > Wall time: 42.89 s > > In [16]: %time a=[p for p in points if poly.contains(p)] > CPU times: user 42.45 s, sys: 0.16 s, total: 42.61 s > Wall time: 42.68 s > > In [18]: %time a=map(poly.contains,points) > CPU times: user 37.98 s, sys: 0.15 s, total: 38.12 s > Wall time: 38.21 s > > In [20]: %time a=map(poly.disjoint,points) > CPU times: user 37.91 s, sys: 0.20 s, total: 38.10 s > Wall time: 38.18 s > > Is there a better way to do the same? Do you happen to know if the new > prepared geometry could help speeding up a little? > > Thanks and regards, > > Mario > _______________________________________________ > Community mailing list > Community@lists.gispython.org > http://lists.gispython.org/mailman/listinfo/community >
hi, i'd be interested to see how prepared geometries could speed that up. but i've found this: http://matplotlib.sourceforge.net/faq/howto_faq.html?highlight=nxutils#test-whether-a-point-is-inside-a-polygon to be _extremely_ fast. you could also stick your points in an rtree (http://pypi.python.org/pypi/Rtree/) and grab the points falling within the bounds, then do the contains() test on the points that pass that. -brent _______________________________________________ Community mailing list Community@lists.gispython.org http://lists.gispython.org/mailman/listinfo/community
