How to check if a 2D point is in a polygon?

The undocumented Graphics`PolygonUtils`PointWindingNumber (if you're on versions < 10, use Graphics`Mesh`PointWindingNumber) does exactly this — it gives you the winding number of a point. A point lies inside the polygon if and only if its winding number is non-zero.

Using this, you can create a Boolean function to test if a point is inside the polygon

inPolyQ[poly_, pt_] := Graphics`PolygonUtils`PointWindingNumber[poly, pt] =!= 0

(* Examples *)
inPolyQ[{{-1, 0}, {0, 1}, {1, 0}}, {1/3, 1/3}]
(* True *)
inPolyQ[{{-1, 0}, {0, 1}, {1, 0}}, {1, 1}]
(* False *)

Or, you can use the aptly named Graphics`PolygonUtils`InPolygonQ which has the same 2-argument syntax and is a predicate.

Using the function winding from Heike's answer to a related question

 winding[poly_, pt_] := 
 Round[(Total@ Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly), 
  2 Pi, -Pi]/2/Pi)]

to modify the test function in this Wolfram Demonstration by R. Nowak to

testpoint[poly_, pt_] := 
Round[(Total@ Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly), 
    2 Pi, -Pi]/2/Pi)] != 0

gives

Update: Full code:

Manipulate[With[{p = Rest@pts, pt = First@pts},
   Graphics[{If[testpoint[p, pt], Pink, Orange], Polygon@p},
   PlotRange -> 3 {{-1, 1}, {-1, 1}},
   ImageSize -> {400, 475},
   PlotLabel -> Text[Style[If[testpoint[p, pt], "True ", "False"], Bold, Italic]]]],
 {{pts, {{0, 0}, {-2, -2}, {2, -2}, {0, 2}}}, 
 Sequence @@ (3 {{-1, -1}, {1, 1}}), Locator, LocatorAutoCreate -> {4, Infinity}},
 SaveDefinitions -> True,   
 Initialization :> {
 (* test if point pt inside polygon poly *)
    testpoint[poly_, pt_] := 
    Round[(Total@ Mod[(# - RotateRight[#]) &@(ArcTan @@ (pt - #) & /@ poly), 
       2 Pi, -Pi]/2/Pi)] != 0 } ]

Update 2: An alternative point-in-polygon test using yet another undocumented function:

 testpoint2[poly_, pt_] := Graphics`Mesh`InPolygonQ[poly, pt]

 testpoint2[{{-1, 0}, {0, 1}, {1, 0}}, {1/3, 1/3}]
 (*True*)
 testpoint2[{{-1, 0}, {0, 1}, {1, 0}}, {1, 1}]
 (*False*)

Sometimes speed is an issue if there are many polygons and or many points to check. There is an excellent reference on this issue under http://erich.realtimerendering.com/ptinpoly/ with the main conclusion that the angle summation algorithm should be avoided if speed is the objective.

Below is my Mathematica implementation of the point in polygon problem which appears to be roughly 5x faster than the inPolyQ[] algorithm posted above.

Test case - use triangle

poly = {{-1, 0}, {0, 1}, {1, 0}};

My code implementation

inPoly2[poly_, pt_] := Module[{c, nvert,i,j},
   nvert = Length[poly];
   c = False;
   For[i = 1, i <= nvert, i++,
    If[i != 1, j = i - 1, j = nvert];
    If[(
      ((poly[[i, 2]] > pt[[2]]) != (poly[[j, 2]] > pt[[2]])) && (pt[[
      1]] < (poly[[j, 1]] - 
         poly[[i, 1]])*(pt[[2]] - poly[[i, 2]])/(poly[[j, 2]] - 
          poly[[i, 2]]) + poly[[i, 1]])), c = ! c];
    ];
   c
   ];

An the timing output testing on point {0,0.99}

Timing[t1 = Table[inPolyQ[poly, 0, 0.99], {10000}];]
Timing[t2 = Table[inPoly2[poly, 0, 0.99], {10000}];]

Out[115]= {0.062, Null}
Out[116]= {0.016, Null}

Update Following a suggestion from ruebenko I've now investigated the actual performance of all the different point-in-polygon routines for two specific cases.

Test No1: Simple triangle polyon and testing using 5000 random test points

poly = {{-1, 0}, {0, 1}, {1, 0}};
pts = Partition[RandomReal[{-1, 1}, 10000], 2];
npts = Length@pts;
Print["inPoly2: ", 
 Timing[Table[inPoly2[poly, pts[[i]]], {i, npts}];][[1]]]
Print["testpoint: ", 
 Timing[Table[testpoint[poly, pts[[i]]], {i, npts}];][[1]]]
Print["testpoint2: ", 
 Timing[Table[testpoint2[poly, pts[[i]]], {i, npts}];][[1]]]
Print["inPolyQ: ", 
 Timing[Table[inPolyQ[poly, pts[[i]]], {i, npts}];][[1]]]
Print["InsidePolygonQ: ", 
 Timing[Table[InsidePolygonQ[poly, pts[[i]]], {i, npts}];] [[1]]]
Print["inPolyQ2: ", 
 Timing[Table[
     inPolyQ2[poly, pts[[i, 1]], pts[[i, 2]]], {i, npts}];][[1]]]

with the following results

inPoly2: 0.202
testpoint: 0.25
testpoint2: 0.016
inPolyQ: 0.015
InsidePolygonQ: 12.277
inPolyQ2: 0.032

Test No2: Very complicated polygon. The main CountryData[] polygon for Canada has over 10 000 vertices and a fairly complex shape. I've focused on the fastest routines and excluded the InsidePolygonQ[] routine in this case and used 200 test points.

p = CountryData["Canada", "Polygon"][[1, 1]];
poly = {Rescale[p[[All, 1]], {Min@#, Max@#} &@p[[All, 1]], {-1, 1}],
    Rescale[p[[All, 2]], {Min@#, Max@#} &@p[[All, 2]], {-1, 1}]} // 
   Transpose;
pts = Partition[RandomReal[{-1, 1}, 400], 2];
npts = Length@pts;
Print["inPoly2: ", 
 Timing[Table[inPoly2[poly, pts[[i]]], {i, npts}];][[1]]]
Print["testpoint: ", 
 Timing[Table[testpoint[poly, pts[[i]]], {i, npts}];][[1]]]
Print["testpoint2: ", 
 Timing[Table[testpoint2[poly, pts[[i]]], {i, npts}];][[1]]]
Print["inPolyQ: ", 
 Timing[Table[inPolyQ[poly, pts[[i]]], {i, npts}];][[1]]]
Print["inPolyQ2: ", 
 Timing[Table[
 inPolyQ2[poly, pts[[i, 1]], pts[[i, 2]]], {i, npts}];][[1]]]

with the following results

inPoly2: 8.237
testpoint: 11.295
testpoint2: 0.156
inPolyQ: 0.436
inPolyQ2: 0.078

My verdict: There is an astonishing 3 orders of magnitude difference in the performance of the different routines. InsidePolygonQ[] while mathematically elegant, is very slow. It pays to use either the undocumented routine for point in polygon in Mathematica, in this case testpoint2[] (with the usual caveats), or the compiled routine inPolyQ2[] which both had excellent performance for both simple and complex test polygons.