tag:blogger.com,1999:blog-2044555433406398156.post8764472954672783495..comments2017-09-24T16:49:45.719-05:00Comments on GeospatialPython.com: Point in PolygonJoel Lawhead, PMP, GISPhttp://www.blogger.com/profile/02508246012088522732noreply@blogger.comBlogger33125tag:blogger.com,1999:blog-2044555433406398156.post-55380870607325873182017-02-15T10:06:50.274-06:002017-02-15T10:06:50.274-06:00I know it has been years but just wanted to say th...I know it has been years but just wanted to say thank you !Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-60326495289387929122016-06-14T18:41:34.121-05:002016-06-14T18:41:34.121-05:00This worked really well for me and I want to appla...This worked really well for me and I want to applaud your efforts. My only wish at this point is that the program worked a bit faster. I'm dealing with a large file, though. Cheers and nice work!Matt Cremeenshttps://www.blogger.com/profile/14167148241424111480noreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-5370723645520270332015-03-29T18:26:30.223-05:002015-03-29T18:26:30.223-05:00Hi, Many thanks for this great tutorial. Regarding...Hi, Many thanks for this great tutorial. Regarding the second program in which you've generated some points on a shapefile, I need just the coordinates of the points. So I dropped the pngcanvas parts of your code and just stored the px and py values. But when I map them, it gives me straight line. What am I missing here? Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-4596588054445837132014-10-30T11:31:12.242-05:002014-10-30T11:31:12.242-05:00Hi,
thanks for providing the code. I presume that...Hi, <br />thanks for providing the code. I presume that in the function point_in_poly(x,y,poly), the loop could really start at 1 (omitting evaluation of a degenerate line).<br />Btw. point_in_poly(x,y,poly) should work both for clockwise and counter clockwise perimeters. Those who want to understand the algorithm should note that it relies on external point == (+infinity, y) for ray casting in the figure on the top of the pageAnonymousnoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-54752089288326287852013-11-21T11:59:59.804-06:002013-11-21T11:59:59.804-06:00Hi, I was aiming something like that, and I tried ...Hi, I was aiming something like that, and I tried your approach and it actually worked, thanks, but it took too long, around 0.5 secs to calculate in which of 60 polygons was 1 point (but i need to run it with a hundred million points)<br /><br />But soon, I realize, I needed to accelerate it, I was trying to do a point_on_polygon but for a list of more than 100.000 points and 60 (non overlapping) polygons (each of them with hundreds of vertices). So I decided to give it a try to shapely, I read my points and turned them into an array of points, and then did the same for the polygons, and then created a loop like:<br />```<br />points = Array of shapely.Points<br />polygons = Array of shapely.Polygons<br /><br />for p in points:<br /> for m in polygons:<br /> if p.within(m):<br /> whatever<br />```<br /><br />The performance for the point in polygon is spectacular, so if you are working on something like this, you should also try shapely.ahttps://www.blogger.com/profile/02673187612643118001noreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-48480935888454731982013-09-26T05:38:54.492-05:002013-09-26T05:38:54.492-05:00danken!!danken!!Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-21217400410319173702013-05-02T07:57:17.691-05:002013-05-02T07:57:17.691-05:00Ok, but let me show you another thing:
Look at th...Ok, but let me show you another thing:<br /><br />Look at this question, it's about this code:<br />http://stackoverflow.com/questions/16325720/point-in-convex-polygon<br /><br />I'm also trying to implement such a function in python, and this solution is what I've found on most of websites. I'm trying it for different shapes/points, but it's sometimes acting the opposite it should. (I'm not a programmer, just a beginner with python, so I don't know exactly what's wrong)Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-38272305549157016592013-05-01T15:15:24.185-05:002013-05-01T15:15:24.185-05:00You are correct. This algorithm does not account ...You are correct. This algorithm does not account for points along the polygon boundary. Please see the following post from a few months later. I'll put a note on this post to steer people in the right direction:<br /><br />http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.htmlJoel Lawheadhttps://www.blogger.com/profile/15645953215148786163noreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-75857820971120576912013-05-01T14:56:57.263-05:002013-05-01T14:56:57.263-05:00This code does not work. Just try your example for...This code does not work. Just try your example for x=9, y=0 (the polygon is the same: (0,0) (0,10) (10,0) (10,10) )<br /><br />It should give true, but it gives false.Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-41419663159542645292012-12-27T13:43:05.378-06:002012-12-27T13:43:05.378-06:00Is there any exception(s) when the procedure expla...Is there any exception(s) when the procedure explained above, to identify if a point is within or outside a polygon would fail to work? <br /><br />lizthornblom@gmail.comAnonymousnoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-4473996592535329572012-12-27T13:41:52.516-06:002012-12-27T13:41:52.516-06:00The operation “point in polygon” is used every tim...The operation “point in polygon” is used every time a GIS user clicks with the mouse somewhere on a themecontaining polygons. The GIS program must then calculate if the point is inside or outside a polygon in order to identify the correct polygon where the user clicked. Other times this operation is used are e.g. when we want totest which cities that are situated in a certain district, or the distribution of soil fertility samples on different land use polygons.<br /><br />By tracing an indefinitely long line in any direction from a point, if the number of times this line intersects the edges of a polygon is uneven (or odd) (1, 3, 5, 7, …) then the point is inside the polygon. On the other hand, if the number of intersections is even (0, 2, 4, 6, 8, …), the point is outside the polygon. This implies that if we are able to calculate intersections between lines we are also able to check if a point lies inside a specific polygon.<br />This method is always working if the line does not completely overlap a “border line” of the polygon. Then we have a special case, which is also possible to solve. However, this is very rare.<br /><br />Is any exception(s) when the procedure explained above, to identify if a point is within or outside a polygon would fail to work? <br /><br />lizthornblom@gmail.comAnonymousnoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-67072151047773410802012-07-14T02:35:25.641-05:002012-07-14T02:35:25.641-05:00Here is the simplest Python port of the C code I m...Here is the simplest Python port of the C code I mentioned above.<br /><br /># pnpoly_orig - the original pnpoly, the simplest port from the C-version. Written in "simple" Python.<br /><br />def pnpoly_orig(polyx, polyy, x, y):<br /><br /> nvert = len(polyx)<br /> c = False<br /> i = 0<br /> j = nvert-1<br /> while i < nvert:<br /> if ((polyy[i]>y) != (polyy[j]>y)) and (x < (polyx[j]-polyx[i]) / (polyy[j]-polyy[i]) * (y-polyy[i]) + polyx[i]):<br /> c = not c<br /> j=i<br /> i=i+1<br /> return c<br /><br /><br />I'm currently working on a really optimized version, vectorized in numpy and in which edge inverse slopes computation is done only once for each edge, even if the function is called multiple times (I keep them in a list).NicoLarvenoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-74881990892652000042012-07-14T02:29:39.832-05:002012-07-14T02:29:39.832-05:00I see shape.parts is never used in your code but m...I see shape.parts is never used in your code but many shapefiles actually use it to include nested polygons. Maybe this can explain the issue you describe here. The solution is quite simple: take care of the structure of the nested polygons (shape.points *must* be sliced according to the indexes in shape.parts). Also, if you have nested polygons, the same logic applies as if you where testing your point against edges: if (x,y) is inside an odd number of polygons then inside must be True, etc. (e.g. consider an island on a river and you're looking for a point no ground or water.)NicoLarvenoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-55903389904881508122012-05-31T15:12:24.832-05:002012-05-31T15:12:24.832-05:00I recently dig into the point in polygon problem. ...I recently dig into the point in polygon problem. I think I found two major references about this topic:<br /><br />http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html<br />7 lines of C, assume semi-open polygons (means a point one the edge will always match only one polygon if many adjacent ones are to be considered). Note that is does not require min/max tests at all, really smart.<br /><br />while the other discusses about several algorithms:<br />http://erich.realtimerendering.com/ptinpoly/NicoLarvenoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-20952420405941256182012-05-03T12:10:42.426-05:002012-05-03T12:10:42.426-05:00There is a standard method that works every time.
...There is a standard method that works every time.<br /><br />You the visual studio solution downloadable from http://www.inferentialpower.com/Bio/Ing/ObjectParser/.<br /><br />It contains an "InPolygon" method written in C#.<br /><br />The solution also contains an "Imprint" method for imprinting polygons onto bitmaps.David Inghttps://www.blogger.com/profile/02412043685055845045noreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-3554481866774207982012-02-08T09:15:01.496-06:002012-02-08T09:15:01.496-06:00Mike,
There's no conversion or measurement go...Mike,<br /><br />There's no conversion or measurement going on so I think this is fine for lat/long. It should be abstract enough and works on arbitrary polygons so I think it would be ok. Correct me if you discover otherwise.Joel Lawheadhttps://www.blogger.com/profile/15645953215148786163noreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-10477043564137755652012-02-07T22:27:39.068-06:002012-02-07T22:27:39.068-06:00Joel,
I would be happy to provide my script, but ...Joel,<br /><br />I would be happy to provide my script, but I want to make sure I have my head on straight. Don't you need to convert the lat/lon points into points on a Cartesian coordinate system for this point-in-polygon algorithm to work?Mikenoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-31623738166454054962012-02-06T14:13:00.270-06:002012-02-06T14:13:00.270-06:00Mike,
Sorry for the slow reply. I don't have...Mike,<br /><br />Sorry for the slow reply. I don't have anything coded up but here's some pretty straight forward math that may help:<br /><br />http://tinyurl.com/7cqzwdw<br /><br />If you implement this in Python please let me know. I'd love to share that with the community.Joel Lawheadhttps://www.blogger.com/profile/15645953215148786163noreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-26541020779889360932012-02-03T20:00:53.875-06:002012-02-03T20:00:53.875-06:00Joel: Do you have any suggestions for determining ...Joel: Do you have any suggestions for determining (with Python) if a convex or concave polygon runs clockwise? Your script works well, but the set of polygons I have are a mixture of clockwise and counterclockwise. Thank you in advance.Mikenoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-4317580064807962482011-11-07T20:09:05.651-06:002011-11-07T20:09:05.651-06:00Ah - sorry - no, the starting point doesn't ma...Ah - sorry - no, the starting point doesn't matter.Joel Lawheadhttps://www.blogger.com/profile/15645953215148786163noreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-69255897909324073232011-11-07T19:27:13.552-06:002011-11-07T19:27:13.552-06:00Sorry, I was a little ambiguous. What I meant was ...Sorry, I was a little ambiguous. What I meant was do the order of the vertexes of the polygon list matter? Or can they be specified in any clockwise manner. For example, if you had a simple square defined by 4 points, does it matter which vertex is listed first, given that they are listed clockwise?Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-22875266003318067742011-11-07T19:01:45.775-06:002011-11-07T19:01:45.775-06:00If I understand your question, no - as long as the...If I understand your question, no - as long as they end up clockwise before you save. When you add points for a polygon they go inside a nested list called "parts". A point Is just a list of [x,y] inside the parts list. You can create the parts list and add it to the poly method at any time.Joel Lawheadhttps://www.blogger.com/profile/15645953215148786163noreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-76484456760458278632011-11-07T18:10:35.537-06:002011-11-07T18:10:35.537-06:00Does it matter what order the points are put in, a...Does it matter what order the points are put in, as long as they are clockwise?Anonymousnoreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-1317197766329752032011-11-06T14:43:36.213-06:002011-11-06T14:43:36.213-06:00I keep a copy here:
http://pyshp.googlecode.com/f...I keep a copy here:<br /><br />http://pyshp.googlecode.com/files/shapefile.pdfJoel Lawheadhttps://www.blogger.com/profile/15645953215148786163noreply@blogger.comtag:blogger.com,1999:blog-2044555433406398156.post-77924350739787817912011-11-06T11:56:58.565-06:002011-11-06T11:56:58.565-06:00Thanks Joel! I was doing the points left to right,...Thanks Joel! I was doing the points left to right, and got the right answer most of the times but couldn't figure out why certain times the answer was wrong. Big help. When you said spec, where is this spec?Anonymousnoreply@blogger.com