Does point fall within polygon? VBA Function

[sijpie]

This is not a question, just a post of some code to help find if a point falls within a (complex) polygon.

I found the original code written in C here:
Determining Whether A Point Is Inside A Complex Polygon

The page also explains the issues with complex polygons. But the code of course also works for simple polygons.

I have modified the code to run as a function in Excel, and so it can be called as any other function.

The function takes as arguments the point to be checked as a range with X & Y, and a range with all the points of the polygon, X & Y.

In my example I am checking if wells drilled in the past are within a certain area.

Excel 2010
	A	B	C	D	E	F	G	H	I
1	Stage 1 Polygon
2						Well Database	SurfaceCoordinates
3							Stereo70		In Stage1 Area
4		x	y
5		549926	313550			Well Name	EastingY[m]	NorthingX[m]	Yes/No
6		549357	313316			1 A Est	559,008.51	309,204.47	FALSE
7		549923	312176			1 Ar Est	550,000.30	313,130.16	TRUE
8		550650	312573			1 Ar Vest	543,011.50	313,106.58	FALSE
9		551164	313205			1 I Est	556,226.56	312,713.64	TRUE
10		549926	313550			1 K Vest-Blejesti	539,357.63	313,165.48	FALSE
11						2 A Est	559,388.41	309,234.32	FALSE
12						2 Ar Est	549,925.25	313,140.10	TRUE
13						2 Ar Vest	543,080.28	312,970.75	FALSE
14						2 I Est	554,430.45	313,551.15	FALSE
15						2 P Vest	545,170.00	310,055.00	FALSE
16						3 A Est	559,722.95	309,279.51	FALSE
Wells inSt1 (2)

Cell Formulas
Range	Formula
I6	=PointInPolygon(G6:H6,Stage1Poly)

Named Ranges
Name	Refers To	Cells
'Wells inSt1 (2)'!Stage1Poly	='Wells inSt1 (2)'!$B$5:$C$9

To use the function, copy the code below. In Excel press Alt-F11 to open the VBA editor. Rightclick on the workbook name in the top left panel, and select Insert/Module. In the right hand pane, paste the code. Now in the workbook you can use the function as shown above. The range holding the polygone nodes is range Stage1Poly, which is B5:C9

The function code is as follows:

<font face=Courier New><SPAN style="color:#00007F">Function</SPAN> PointInPolygon(rXY <SPAN style="color:#00007F">As</SPAN> Range, rpolyXY <SPAN style="color:#00007F">As</SPAN> Range) <SPAN style="color:#00007F">As</SPAN> <SPAN style="color:#00007F">Boolean</SPAN><br>    <SPAN style="color:#007F00">' Function checks if X,Y given in rXY falls within complex _<br>      polygon as defined by node list rpolyXY. _<br>      rXY to be 2 cell range with one X and one Y value _<br>      rpolyXY to be 2 column range with for each node on the polygon _<br>      the X and the Y point</SPAN><br>    <br>    <SPAN style="color:#00007F">Dim</SPAN> i <SPAN style="color:#00007F">As</SPAN> <SPAN style="color:#00007F">Integer</SPAN>, j <SPAN style="color:#00007F">As</SPAN> <SPAN style="color:#00007F">Integer</SPAN>, polySides <SPAN style="color:#00007F">As</SPAN> <SPAN style="color:#00007F">Integer</SPAN><br>    <SPAN style="color:#00007F">Dim</SPAN> oddNodes <SPAN style="color:#00007F">As</SPAN> <SPAN style="color:#00007F">Boolean</SPAN><br>    <SPAN style="color:#00007F">Dim</SPAN> x <SPAN style="color:#00007F">As</SPAN> <SPAN style="color:#00007F">Double</SPAN>, y <SPAN style="color:#00007F">As</SPAN> <SPAN style="color:#00007F">Double</SPAN><br>    <SPAN style="color:#00007F">Dim</SPAN> aXY <SPAN style="color:#00007F">As</SPAN> <SPAN style="color:#00007F">Variant</SPAN><br>    <br>    oddNodes = <SPAN style="color:#00007F">False</SPAN><br>    x = rXY.Cells.Value2(1, 1)<br>    y = rXY.Cells.Value2(1, 2)<br>    aXY = rpolyXY.Value<br>    <br>    polySides = rpolyXY.Rows.Count<br>    j = polySides - 1<br>    <br>    <SPAN style="color:#00007F">For</SPAN> i = 1 <SPAN style="color:#00007F">To</SPAN> polySides<br>        <SPAN style="color:#00007F">If</SPAN> (((aXY(i, 2) < y And aXY(j, 2) >= y) _<br>            <SPAN style="color:#00007F">Or</SPAN> (aXY(j, 2) < y And aXY(i, 2) >= y)) _<br>            And (aXY(i, 1) <= x <SPAN style="color:#00007F">Or</SPAN> aXY(j, 1) <= x)) <SPAN style="color:#00007F">Then</SPAN><br>            oddNodes = oddNodes Xor (aXY(i, 1) + (y - aXY(i, 2)) / (aXY(j, 2) - aXY(i, 2)) * (aXY(j, 1) - aXY(i, 1)) < x)<br>        <SPAN style="color:#00007F">End</SPAN> <SPAN style="color:#00007F">If</SPAN><br>        j = i<br>    <SPAN style="color:#00007F">Next</SPAN> i<br>    PointInPolygon = oddNodes<br><SPAN style="color:#00007F">End</SPAN> <SPAN style="color:#00007F">Function</SPAN><br></FONT>

 

[chrism216]

KMZ are Google Earth files
I work at a company which gives maintenance to telecommunications equipment. I use Google Earth extensively to visualize work loads, by mapping the towers we have to visit each month, so i can sectorize efficientlyj and schedule the work so they drive the least distance possible.
The size of the area each technician can cover depends on how many towers there are, so i always have to count. But it gets a bit tiresome since theres thousands of towers...
You can imagine how helpful this will be.

 

[Kyle123]

On the subject of KMZ files, I use google fusion tables for finding a point in a polygon since it's baked in and very fast - just another possible option

 

[chrism216]

I have no idea what KMZ layers are, nor am I sure whose function you found useful (mine or AussieChuck's) in dealing with them, but I am glad you found something in this thread helpful.


I am sure AussieChuck would join me in saying, "You are quite welcome!"

Thanks again Rick

Over the course of this morning I programmed a little something to help me turn what once was pure visualization with no other analysis than "it looks good", to real, accurate data.

Using your function over 130,000 times, iterating through about 2700 geographical locations, checking with over 50 polygons, the end result is a clean database telling me which points are in which polygons. And everytime I tweak a polygon, I can automatically run an analysis to determine the impact on work loads. Just incredible.

Thank you so much.

Chris.

 

[Rick Rothstein]

Thanks again Rick

Over the course of this morning I programmed a little something to help me turn what once was pure visualization with no other analysis than "it looks good", to real, accurate data.

Using your function over 130,000 times, iterating through about 2700 geographical locations, checking with over 50 polygons, the end result is a clean database telling me which points are in which polygons. And everytime I tweak a polygon, I can automatically run an analysis to determine the impact on work loads. Just incredible.

Thank you so much.

You are quite welcome... I am glad I was able to be of help to you.

 

[Trikgard]

Mr. Rothstein,

I've been working with your point in polygon macro for a couple of days trying to resolve an issue. While your macro works for a point (or series of points) for a single polygon, I am trying to compare a multitude of points to a various polygons and determine each polygon as point may fall in or identify that the point is outside of all the polygons. I have my worksheet set up as a matrix (points top-bottom in the first column) and the polygons left-right in two-column sets (x (latitude) = first column, y (longitude) = second column). The macro then runs for each intersection, comparing the point with each polygon left-right. Also, the project needs to be dynamic, as the points can change (and do) and the polygon vertices can change.

In the event that two polygons are next two each other (left-right spatially), if a point is declared 'true' for the first polygon - it also comes up as 'true' for the second polygon as long as the point is below the most northern vertice in the polygon. (Obviously this is a bad thing).

I figured out that the macro is based on the 'looking up to determine the number of lines crossed) - which would same a left-right boundary doesn't come into play. I've thought about using a paired-"large/small" function set to determine left-most / right-most points and make sure the point falls somewhere in between - not sure if that would work.

Any assistance you can provide would be greatly appreciated.

 

[Rick Rothstein]

Mr. Rothstein,

I've been working with your point in polygon macro for a couple of days trying to resolve an issue. While your macro works for a point (or series of points) for a single polygon, I am trying to compare a multitude of points to a various polygons and determine each polygon as point may fall in or identify that the point is outside of all the polygons. I have my worksheet set up as a matrix (points top-bottom in the first column) and the polygons left-right in two-column sets (x (latitude) = first column, y (longitude) = second column). The macro then runs for each intersection, comparing the point with each polygon left-right. Also, the project needs to be dynamic, as the points can change (and do) and the polygon vertices can change.

In the event that two polygons are next two each other (left-right spatially), if a point is declared 'true' for the first polygon - it also comes up as 'true' for the second polygon as long as the point is below the most northern vertice in the polygon. (Obviously this is a bad thing).

I figured out that the macro is based on the 'looking up to determine the number of lines crossed) - which would same a left-right boundary doesn't come into play. I've thought about using a paired-"large/small" function set to determine left-most / right-most points and make sure the point falls somewhere in between - not sure if that would work.

I am not 100% sure of the problem you are trying to describe, but it almost sounds like you are trying to pass multiple polygons into the function at the same time. If so, that won't work. The function can only test a single point against a single polygon. If you need to make some judgment regarding multiple polygons and/or multiple points, you would need to do them one at a time and accumulate the results as each polygon is processed.

 

[gonurvia]

Hi, Id like to know if its possible to use this function por a 3D polygon (sloping surface). I'd like to use it in order to know if the points I introduce lay straight above/below the surface with 3 coordinates. I'd appreciate any help!

 

[sijpie]

I am quite sure it is. I also suspect that Rick Rothstein's function will be more appropriate. but you will have to do some maths yourself.

I think that if you take any plane through your point and which intersects your polygon, you can then use Rick's function to establish in- or outside.

 

[gonurvia]

Hi sijpie, thnx for the answer. I managed to solve it earlier by taking the x,y coord to know if the point lies inside/outside the area. After just by doing some math, I take x,y,z coord of 3 points in the plane, from which I obtain the plane equation and afterwards just simply check by introducing the x,y coord. of my point, find the height of the plane on that point. Compare heights of plane and point in those x,y coord, and I may then find out if its above/below.

 

[Rick Rothstein]

Hi, Id like to know if its possible to use this function por a 3D polygon (sloping surface). I'd like to use it in order to know if the points I introduce lay straight above/below the surface with 3 coordinates. I'd appreciate any help!

The answer to your question kind of depends on the nature of your 3D polygon. You could use my code to determine if the point lies within the planar polygon by zeroing our the Z-axis (assuming an XYZ coordinate system) or zeroing out the elevation assuming a survey Northing, Easting, Elevation system... but that is not your question... you want to know if point is in the air or in the ground (assuming a survey type system). Introducing the third dimension add an enormous amount of complexity. For example, let's say your polygon is composed of 10 angle points... are the all at the same elevation, probably not... are the 10 points all guaranteed to lie on a flat plane that is tilted to the Z axis? If they were, what about the area encompassed by them... do we have a perfectly flat plane or an undulating enclosed surface. If perfectly flat, I'd say we could calculate what you asked for somewhat easily... however, if the surface undulated, then you would need to triangulate the enclosed area (a fairly complex procedure), test each triangle to find which one it is above and then calculate if the point is above or below that triangle. So then, what can you tell us about the nature of the 3D polygon?