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

Outcomes