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

[Conrad Jackisch]

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)
for p in points:
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/62047500/attachment.htm>

[Sean Gillies]

Hi Conrad,

The key thing to understand about Shapely's operations (intersection,
union, etc) is that they involve sets of points in space and not sets of
Python objects. A.intersection(B) returns a Python object that represents
all the points in the intersection of A's points and B's points.

If you want just a point-in-polygon test and to track your points
precisely, you could do something like this:

from shapely.prepared import prep

poly = prep(poly)

for xy in a:
if poly.intersects(Point(xy)):
# plot it

Yours,

Post by Conrad Jackisch
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
_______________________________________________
Community mailing list
Community at lists.gispython.org
http://lists.gispython.org/mailman/listinfo/community

[Conrad Jackisch]

Dear Sean,

thanks.
i fear i wont get much faster than this:

from shapely.prepared import prep
from shapely.geometry import Point
from shapely.geometry import Polygon
from shapely.geometry import MultiPoint
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)))

poly_prep = prep(poly)
points_prep = prep(points)

hits = map(poly_prep.contains, points)
#from timeit import Timer
#t = Timer(lambda: map(poly_prep.contains, points))
#print t.timeit(number=1)

#t = Timer(lambda: points.intersection(poly))
#print t.timeit(number=1)

?which is about a second on my machine.
i will look for some other approaches instead.

thanks a million for your reply,
cheers.


2014/1/8 Sean Gillies <sean.gillies at gmail.com>

Post by Sean Gillies
Hi Conrad,
The key thing to understand about Shapely's operations (intersection,
union, etc) is that they involve sets of points in space and not sets of
Python objects. A.intersection(B) returns a Python object that represents
all the points in the intersection of A's points and B's points.
If you want just a point-in-polygon test and to track your points
from shapely.prepared import prep
poly = prep(poly)
# plot it
Yours,

Post by Conrad Jackisch
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
_______________________________________________
Community mailing list
Community at lists.gispython.org
http://lists.gispython.org/mailman/listinfo/community

--
Sean Gillies
_______________________________________________
Community mailing list
Community at lists.gispython.org
http://lists.gispython.org/mailman/listinfo/community

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.gispython.org/pipermail/community/attachments/20140109/8501e944/attachment.htm>