Since someone dragged in Canada...

Here is the code from a MathGroup post I had referenced. I have modified to compile to C and that speeds it further. The one-off preprocessing does take time but it seems not unreasonable. It takes a list of lists of polygons (so the "region" need not be connected). To account for this I slightly alter the setup from Mac's response.

Preprocessing the polygon:

```
getSegsC =
Compile[{{j, _Integer}, {minx, _Real}, {len, _Real}, {eps, _Real}, \
{segs, _Real, 3}}, Module[{lo, hi}, lo = minx + (j - 1)*len - eps;
hi = minx + j*len + eps;
Select[segs,
Module[{xlo, xhi}, {xlo, xhi} = Sort[{#[[1, 1]], #[[2, 1]]}];
lo <= xlo <= hi ||
lo <= xhi <= hi || (xlo <= lo && xhi >= hi)] &]]];
polyToSegmentList[poly_, nbins_] :=
Module[{xvals, yvals, minx, maxx, miny, maxy, segments, flatsegments,
segmentbins, xrange, len, eps}, {xvals, yvals} =
Transpose[Flatten[poly, 1]];
{minx, maxx} = {Min[xvals], Max[xvals]};
{miny, maxy} = {Min[yvals], Max[yvals]};
segments = Map[Partition[#, 2, 1, {1, 1}] &, poly];
flatsegments = Flatten[segments, 1];
xrange = maxx - minx;
eps = 1/nbins*len;
len = xrange/nbins;
segmentbins =
Table[getSegsC[j, minx, len, eps, flatsegments], {j, nbins}];
{{minx, maxx}, {miny, maxy}, segmentbins}]
```

The actual in-or-out code.

```
pointInPolygon[{x_, y_}, bins_, xmin_, xmax_, ymin_, ymax_] :=
Catch[Module[{nbins = Length[bins], bin},
If[x < xmin || x > xmax || y < ymin || y > ymax, Throw[False]];
bin = Ceiling[nbins*(x - xmin)/(xmax - xmin)];
If[EvenQ[countIntersectionsC[bins[[bin]], x, y, ymin - 1.]], False,
True]]]
countIntersectionsC =
Compile[{{segs, _Real, 3}, {x, _Real}, {yhi, _Real}, {ylo, _Real}},
Module[{tally = 0, yval, xlo, xhi, y1, y2},
Do[{{xlo, y1}, {xhi, y2}} = segs[[j]];
If[(x < xlo && x < xhi) || (x > xlo && x > xhi), Continue[]];
yval = y1 + (x - xlo)/(xhi - xlo)*(y2 - y1);
If[ylo < yval < yhi, tally++];, {j, Length[segs]}];
tally]];
```

The mainland of Canada will be the test again. As in Mac's example I rescale so coordinates are all between -1 and 1. This means I really don't need the x/ymin/max stuff but I opted to keep that in.

```
p = CountryData["Canada", "Polygon"][[1, 1]];
poly = {Transpose[{Rescale[
p[[All, 1]], {Min@#, Max@#} &@p[[All, 1]], {-1, 1}],
Rescale[p[[All, 2]], {Min@#, Max@#} &@p[[All, 2]], {-1, 1}]}]};
```

I'll use 1000 bins and do the preprocessing.

```
nbins = 1000;
Timing[{{xmin, xmax}, {ymin, ymax}, segmentbins} =
polyToSegmentList[poly, nbins];]
(* Out[369]= {5.15, Null} *)
```

For testing I'll take 10000 points.

```
npts = 10000;
pts = Partition[RandomReal[{-1, 1}, 2*npts], 2];
Timing[
inout = Map[pointInPolygon[#, segmentbins, xmin, xmax, ymin, ymax] &,
pts];]
(* Out[402]= {0.37, Null} *)
```

Visual check:

```
ListPlot[Pick[pts, inout], Joined -> False]
```

The northeast reminds me a bit of the duck's head seen here. But then...I've always found the Baffin...to be bafflin'.

Notebook from Mathematica-users.org PointInsidePolygon.nb contains additional solutions.

– Alexey Popkov – 2013-10-11T00:14:01.3001

I don't believe Mathematica has a build-in function for this. You could build your own in which case this is a good starting point: point in polygon

– jVincent – 2012-08-14T08:39:00.740Perhaps this answer by Heike?

– kglr – 2012-08-14T08:49:56.780@jVincent, I was hoping someone was willing to share their Mathematica code. – nilo de roock – 2012-08-14T08:52:17.610

@kguler Interesting stuff in the answers. Thanks, I'll have a closer look. – nilo de roock – 2012-08-14T08:53:31.913

2It's pretty depressing that you already accepted an answer before I posted mine, but I posted anyway. :-/ – Mr.Wizard – 2012-08-14T10:02:55.690

Anyway: there are a number of methods to choose from here (also here). Now, adapting them to

– J. M.'s ennui – 2012-08-14T13:27:58.140Mathematicais a different kettle of fish...4The best answer to this question depends on the intended use and other constraints. The most important determinants are (1) whether this will be a one-off test or if many points will be tested for a given polygon; (2) where the points are likely to fall; and (3) whether the test needs to be absolutely accurate. Except for the one-off accurate test, by

farthe fastest method--and one not yet offered in any answer--is to rasterize the polygon's interior, probe the raster at the point's location (a $O(1)$ operation), and revert to a more expensive test only if the probe is inconclusive. – whuber – 2012-08-14T17:52:25.6534

The various responses are quite good. That said, if you have to test many points and the polygon has many vertices, the method indicated at this link should be fairly efficient. http://forums.wolfram.com/mathgroup/archive/2009/Feb/msg00519.html

– Daniel Lichtblau – 2012-08-14T21:10:18.463