[Community] speeding up "contains"

[Sean Gillies]

Hi Mario,

Masking is probably better left to image or raster processing  
software. If you can think of any way that Shapely might work better  
with Gimp, GRASS, or whatever, do let me know.

Sean

On May 11, 2009, at 4:36 AM, Mario Ceresa 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 at lists.gispython.org
> http://lists.gispython.org/mailman/listinfo/community