AnsweredAssumed Answered

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

Question asked by stephan.marc.garland on Nov 10, 2016
Latest reply on Nov 10, 2016 by CDow-esristaff

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))
#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.