On Tue, Oct 14, 2008 at 8:46 AM, Pierre GM <[EMAIL PROTECTED]> wrote:
>> 2008/10/13 Mathew Yeates <[EMAIL PROTECTED]>
>>
>> > Is there a routine in scipy for telling whether a point is inside a
>> > convex 4 sided polygon?
>
> Mathew,
> You could use OGR (www.gdal.org)
>
> Example
> -------------
> import osgeo.ogr as ogr
>
> vert = [(0,0),(0,1),(1,1),(1,0)]
> listvert = [" %s %s" % (x,y) for (x,y) in vert]
> listvert.append(listvert[0])
> geo = ogr.CreateGeometryFromWkt("POLYGON ((%s))" % ','.join(listvert))
>
> querypoint = (0.5, 0.5)
> qpt = ogr.CreateGeometryFromWkt("POINT(%s %s)" % querypoint)
>
> assert geo.Contains(qpt)
>
> querypoint = (0.5, 1.5)
> qpt = ogr.CreateGeometryFromWkt("POINT(%s %s)" % querypoint)
>
> assert not geo.Contains(qpt)
If all you really need is a point in convex polygon test, it's
probably a little too trivial to be worth dragging in a dependency on
anything. But if you may find yourself needing more geometric tests
then it might be a good idea to get a good geom lib now.
As for this test, all you need to do is check that the point is to the
left of each of the edges, taken counter-clockwise.
Here's some code I wrote a while back that does a more general
even-odd test, but should work for your case too:
def inside_shape(p, verts, edges=None):
"""Test whether the point p is inside the specified shape.
The shape is specified by 'verts' and 'edges'
Arguments:
p - the 2d point
verts - (N,2) array of points
edges - (N,2) array of vert indices indicating edges
If edges is None then assumed to be in order. I.e.
[[0,1], [1,2], [2,3] ... [N-1,0]]
Returns:
- True/False based on result of in/out test.
Uses the 'ray to infinity' even-odd test.
Let the ray be the horizontal ray starting at p and going to +inf in x.
"""
verts = npy.asarray(verts)
if edges is None:
N = verts.shape[0]
edges = npy.column_stack([npy.c_[0:N],npy.c_[1:N,0]])
inside = False
x,y=p[0],p[1]
for e in edges:
v0,v1 = verts[e[0]],verts[e[1]]
# Check if both verts to the left of ray
if v0[0]<x and v1[0]<x:
continue
# check if both on the same side of ray
if (v0[1]<y and v1[1]<y) or (v0[1]>y and v1[1]>y):
continue
#check for horizontal line - another horz line can't intersect it
if (v0[1]==v1[1]):
continue
# compute x intersection value
xisect = v0[0] + (v1[0]-v0[0])*((y-v0[1])/(v1[1]-v0[1]))
if xisect >= x:
inside = not inside
return inside
License: public domain.
--bb
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
