Re: testing if a point is inside a polygon
- To: mathgroup at smc.vnet.net
- Subject: [mg96233] Re: [mg96189] testing if a point is inside a polygon
- From: Frank Scherbaum <Frank.Scherbaum at geo.uni-potsdam.de>
- Date: Tue, 10 Feb 2009 05:48:09 -0500 (EST)
- References: <200902091032.FAA12225@smc.vnet.net>
Mitch,
please find below what I use for his purpose. I hope it is useful.
Best regards,
Frank
Am Feb 9, 2009 um 11:32 AM schrieb Mitch Murphy:
>
> is there a way to test whether a point is inside a polygon? ie.
>
> PointInsidePolygonQ[point_,polygon_] -> True or False
>
> i'm trying to do something like ...
>
> ListContourPlot[table,RegionFunction-
> >CountryData["Canada","Polygon"]]
>
> to create what would be called a "clipping mask" in photoshop.
>
> cheers,
> Mitch
>
PointInPolygonQ::usage="PointInPolygonQ[pt, poly] uses the winding-
number algorithm (Godkin and Pulli, 1984) to check, if point pt is
inside the closed polygon poly, which is given as list of its
vertices."
(*
checks, if a point is inside a polygon
pt: point as {lat [deg], lon [deg]} to test
poly: list of polygon vertices coordinates
GODKIN,C.B. AND J.J.PULLI: APPLICATION OF THE "WINDING-NUMBER
ALGORITHM" TO THE SPATIAL SORTING OF CATALOGED EARTHQUAKE DATA.
Bull. Seismol. Soc. Am. 74, 5, PP. 1845-1848, OCTOBER 1984
RETURN VALUE: 0 IF POINT OUTSIDE
+/-1 IF POINT INSIDE
2 IF POINT IS ON AN EDGE OR VERTEX
*)
PointInPolygonQ[pt_,poly_] := Module[
{ i,n,isicr,inside,px,py,pxx,pyy,x0,y0 },
n = Length[poly];
(* ACCUMULATE SIGNED CROSSING NUMBERS WITH INSIDE *)
inside = 0;
{x0,y0}=pt;
For[i=1,i < n,i++,
(* PROCEED AROUND POLYGON CHECKING EACH SEGMENT TO SEE IF
NEGATIVE X-AXIS WAS CROSSED
TRANSLATE COORDINATES OF POLYGON TO PUT TEST POINT AT
ORIGIN *)
{px,py} = poly[[i]];
{pxx,pyy} = poly[[i+1]];
isicr = ksicr[px - x0, py - y0, pxx - x0, pyy - y0];
(* STOP COUNTING IF POINT IS ON EDGE *)
If[isicr == 4, Return[2]];
inside += isicr;
];
(* CHECK SEGMENT FROM LAST VERTEX TO FIRST VERTEX *)
{px,py} = poly[[n]];
{pxx,pyy} = poly[[1]];
isicr = ksicr[px - x0, py - y0, pxx - x0, pyy - y0];
If[isicr == 4, Return[2]];
inside = (inside + isicr)/2;
Return[inside];
];
(*
COMPUTE SIGNED CROSSING NUMBER
A "SIGNED CROSSING NUMBER" TELLS WETHER A SEGMENT
(IN THIS CASE THE SEGMENT FROM (X1,Y1) TO (X2,Y2))
CROSSES THE NEGATIVE X-AXIS OR GOES THROUGH THE ORIGIN
THE RETURN VALUES ARE:
+2 IF SEGMENT CROSSES FROM BELOW
+1 IF SEGMENT EITHER ENDS ON -X-AXIS FROM BELOW OR STARTS
UPWARDS FROM -X-AXIS ("HALF CROSSING")
0 IF THERE IS NO CROSSING
-1 IF SEGMENT EITHER ENDS ON -X-AXIS FROM ABOVE OR STARTS
DOWNWARDS FROM -X-AXIS ("HALF CROSSING")
-2 IF SEGMENT CROSSES FROM ABOVE
+4 IF SEGMENT CROSSES THROUGH THE ORIGIN
*)
ksicr[x1_,y1_,x2_,y2_] := Module[
{
},
(* IF BOTH POINTS ARE ON THE SAME SIDE OF X-AXIS, RETURN 0 *)
If[N[y1*y2 > 0.], Return[0] (* no crossing *)];
(* CHECK IF SEGMENT CROSSES THROUGH THE ORIGIN *)
If[x1*y2 != x2*y1 || x1*x2 > 0.,
If[y1 * y2 < 0,
(*
COMPLETE CROSSING OF -X-AXIS?
BREAK INTO CASES ACCORDING TO CROSSING DIRECTION
*)
If[y1 > 0,
(* CASE Y1 > 0 > Y2 *)
If[y1 * x2 >= x1 * y2, Return[0];, (* no crossing *)
Return[-2]; (* downward crossing *)
];
,
(* CASE Y1 < 0 < Y2 *)
If[x1 * y2 >= y1 * x2, Return[0];, (* no crossing *)
Return[2]; (* upward crossing *)
];
];
,
(*
HALF CROSSING?
ONE END OF SEGMENT TOUCHES X-AXIS! WHICH END?
*)
If[y2 == 0,
(* HERE Y2==0 CHECK IF SEGMENT TOUCHES +X-AXIS *)
If[y1 == 0 || x2 > 0,
Return[0]; (* no crossing *)
,
(* UPWARD OR DOWNWARD? *)
If[y1 > 0.,
Return[-1]; (* Downward half crossing *)
,
Return[1]; (* Upward half crossing *)
];
];
,
(* HERE Y1==0 CHECK IF SEGMENT TOUCHES +X-AXIS *)
If[x1 > 0,
Return[0];
,
(* UPWARD OR DOWNWARD? *)
If[y2 > 0,
Return[1]; (* Upward half crossing *)
,
Return[-1]; (* Downward half crossing *)
];
];
];
(* HERE Y1=0 CHECK IF SEGMENT TOUCHES +X-AXIS *)
If[x1 > 0,
Return[0]; (* no crossing *)
];
(* UPWARD OR DOWNWARD? *)
If[y2 > 0.,
Return[-1]; (* Downward half crossing *),
Return[1]; (* Upward half crossing *)
];
];
,
Return[4];
];
];
- References:
- testing if a point is inside a polygon
- From: Mitch Murphy <mitch@lemma.ca>
- testing if a point is inside a polygon