[Community] Follow-up: Q: py:shapely :: is there a way to get index of multipoint intersection w/ polygon

[Conrad Jackisch]

dear community,

of cause there is the simple solution using .intersects and looping over
all points...
but this really takes too long.

- - - 8< snip - - -
poly=Polygon(zip([0.0,0.3,0.3,0.2,0.0],[0.2,0.2,0.4,0.5,0.4]))
points = MultiPoint(zip(np.random.rand(100000),np.random.rand(100000)))

#intersection
pm = np.ma.array(points, mask=False)
for i in np.arange(len(points)):
pm.mask[i]=points[i].intersects(poly)

#plot example
#warning: takes a while to render...
ax = plt.figure(figsize=(8, 4))
ax = plt.subplot(121)

patch = PolygonPatch(poly)
ax.add_patch(patch)
for p in points:
ax.plot(p.x, p.y, 'o', color='#999999')

for p in pm:
ax.plot(p[0],p[1], 'o', color='#A0522D')

pyplot.show()

- - - snap >8 - - -

so any idea with regard to speed was fabulous.
cheers,
conrad

2014/1/8 Conrad Jackisch <cojacoo at gmail.com>

dear community,
i am struggling with a intersection of a large number of points with a
polygon. i use shapely which helps already but returns only the coordinates
of intersecting points of a MultiPoint object. is there anyone who could
help with a way to simply get the indexes of intersecting points?
- - - 8< snip - - -
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
from shapely.geometry import asMultiPoint
poly=Polygon(zip([0.0,0.3,0.3,0.2,0.0],[0.2,0.2,0.4,0.5,0.4]))
a = array([[0.1,0.1],[0.2,0.3],[0.15,0.45]])
points = asMultiPoint(a)
#plot example
ax = plt.figure(figsize=(8, 4))
ax = plt.subplot(121)
patch = PolygonPatch(poly)
ax.add_patch(patch)
ax.plot(p.x, p.y, 'o', color='#999999')
pyplot.show()
#intersection
np.asarray(points.intersection(poly))
#this gives me the coordinates of the values - but not the index, which i
need much more.
- - - snap >8 - - -
of cause this example is very small and it was possible to compare the
results with the given points. however, for large sets (about one million
of points) this will slow down everything...
if you have any idea, i was most happy.
thanks,
conrad

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.gispython.org/pipermail/community/attachments/20140108/7eef540b/attachment.htm>