Find point in Polygon/MultiLineStrings with multiple interior areas, using Shapely and JSON

[StephanGarland]

X-posted from StackExchange.

I've been working on this for a few days, adding in features, and fixing the problems as they crop up.

The goal is to read in a CSV containing addresses for locates, parse/normalize them (that functionality will come later), geocode them, and compare them to our service territory, which I currently have as both a KML and GeoJSON (I used Ogre to conver the KML). I've uploaded the JSON to PasteBin. I've made a very simplified program that only reads in a point and compares it to the JSON for now. I've also tried a simple program using fiona and a .shp file, same result.

Our service area is bounded on the outside, and within it, there are several other areas that we do not service, with the exception of two - those two I've added placemarkers, which I assume I could check for. My first thought was to use Shapely's contains function, but after trying, and reading the documentation, I don't know if it would deal with the MultiLineStrings making up my boundaries.

The only way I get a positive hit out of this program is if I put in one of the two Points I made. I've tried inside the main boundary, inside a smaller bounding box (non-serviced area) inside that, and inside the boxes that we do service.

import json, shapely, time

from shapely.geometry import Point, shape#, MultiLineString

jsonFile = 'P:/Territory.json'

point = Point(-86.76642, 40.67597)

with open(jsonFile, 'r') as jsonFile:
js = json.load(jsonFile)

for feature in js['features']:
#polygon = MultiLineString(feature['geometry'])
polygon = shape(feature['geometry'])
if polygon.contains(point):
print("In our territory." + str(feature))
time.sleep(2)
#else:
#print("Not in our territory.")

The MultiLineString was throwing a KeyError: 0 (I can include the entire traceback if needed), so I commented it out and tried polygon. The else is commented out because as written, it checks every feature for a match, and generates a huge list of Nopes. I can fix that later; right now I just want to get the point in polygon working correctly. It's also worth mentioning that I'm happy to use any Python library and file format to accomplish this, arcpy included.

[ClintonDow1]

In the JSON the 'geometry' key's value is not consistent. Some of them are points (why you're getting the key error), so you'll want to handle it conditionally. Also many of the geometries are only a line between 2 points, you would need to group the lines that form a closed polygon into a collection, then pass all of their points to the polygon constructor. 

[StephanGarland]

Thank you, I'll start working on grouping them into a collection. That'll be awhile...

[ClintonDow1]

The 'object.touches()' predicate is useful for that, a polygon can have a minimum of three edges. So you can use touches to check if the third line and beyond touches both the first line and the prior line. If so that line closes the polygon. 

[StephanGarland]

I started using Dissolve in ArcMap, which then led me to discover that I have a few instances where the lines don't actually meet, so I'm now fixing those. I think Dissolve will handle it once the lines are closed, though.

[ClintonDow1]

Yes that will definitely work too. There is an environment setting XYTolerance which defines the distance before points 'snap together' that might make it easier as well.