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.gill...@gmail.com> > 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, > > > On Wed, Jan 8, 2014 at 2:29 AM, Conrad Jackisch <coja...@gmail.com> wrote: > >> 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 >> >> _______________________________________________ >> Community mailing list >> Community@lists.gispython.org >> http://lists.gispython.org/mailman/listinfo/community >> >> > > > -- > Sean Gillies > > _______________________________________________ > Community mailing list > Community@lists.gispython.org > http://lists.gispython.org/mailman/listinfo/community > >
_______________________________________________ Community mailing list Community@lists.gispython.org http://lists.gispython.org/mailman/listinfo/community
