0

I am uniformly generating points within the contours of Manhattan, using this geojson file.

I plotted the result using Bokeh to make sure that my output was looking good, and this is what I got:

enter image description here

...so there is no east side?

Here is the full code.

The only ideas I have about this are:

  1. My point-in-polygon calculations are somehow incorrect, in this weirdly specific way.
  2. The GeoJSON definition file crosses itself over somehow and causes an invisible hole in itself (kind of like what you'd get with an SVG).

I need help debugging this, it's such a weirdly specific error and I'm not sure what to blame.

Here is a MWE: https://gist.github.com/ResidentMario/d5e67b6c94988fdc168f.

Here is another point cloud using a much simpler GeoJSON representation of Manhattan:

enter image description here

Here are the diagrams, from the MWE:

enter image description here

Aleksey Bilogur
  • 3,686
  • 3
  • 30
  • 57
  • Welcome to StackOverflow. Please read and follow the posting guidelines in the help documentation. [Minimal, complete, verifiable example](http://stackoverflow.com/help/mcve) applies here. – Prune Feb 08 '16 at 23:39
  • Hmm good point. I'll grab specifically the bad bits, hold on... – Aleksey Bilogur Feb 08 '16 at 23:41
  • Also, reduce that large example to the minimal set necessary to reproduce the problem. – Prune Feb 08 '16 at 23:43
  • I've posted a MWE. I'm not sure that simplifying the shape would help too much, since if it does I still don't know if the problem is in the method or the underlying geojson. – Aleksey Bilogur Feb 08 '16 at 23:59
  • Looks like the East River to me. – stark Feb 09 '16 at 00:07
  • Please post your code here. – Barmar Feb 09 '16 at 00:36
  • All of my code is available, both the full code above and a minimum working example below (GitHub and Gist respectively). I think it's a problem with the underlying mapping. I constructed my own map using [geojson](http://geojson.io/) by hand it had no issues (bonus, performance was an order of magnitude better---fewer points). I still don't know what the issue is with this one map in particular, but I've been able to work around it. – Aleksey Bilogur Feb 09 '16 at 00:58
  • I'd put `if p1y != p2y` at the top of the `if` chain. Where it's at now, you're using a leftover value of `xints`. And people are complaining about this question because you shouldn't have to go off site to see the relevant code. – Mark Ransom Feb 09 '16 at 03:14
  • I *believe* that I have discovered what the issue is: point-in-polygon calculations will not work when there are multiple polygons defined using a single node list. So in this case the presence of linked-but-separate polygons for all of those small polygons messed up the calculation. This is an artifact of the way that the data was constructed: if the data was broken up into several polygons explicitly using a "FeatureSet", instead of represented monolithically as a single "Feature", I suspect the problem would disappear. However, I have not confirmed this fact. – Aleksey Bilogur Feb 10 '16 at 05:57
  • Putting multiple polygons into a single list shouldn't affect the point-in-polygon function. The line that connects one polygon to another is cancelled by the line that goes to the same points in the opposite direction. I think the problem is two adjacent points in the polygon with the same Y coordinate. – Mark Ransom Feb 10 '16 at 18:08

0 Answers0