10

I have tried several examples of code using libraries such as shapefile, fiona, and ogr to attempt to check whether a point (x, y) falls within the boundaries of a multipolygon created with ArcMap (and thus in shapefile format). However none of the examples work well with multipolygons, although they do fine with regular, single polygon shapefiles. Some snippets I tried are below:

# First example using shapefile and shapely:
from shapely.geometry import Polygon, Point, MultiPolygon
import shapefile

polygon = shapefile.Reader('shapefile.shp') 
polygon = polygon.shapes()  
shpfilePoints = []
for shape in polygon:
    shpfilePoints = shape.points 
polygon = shpfilePoints 
poly = Polygon(poly)

point = Point(x, y)
# point in polygon test
if polygon.contains(point):
    print 'inside'
else:
    print 'OUT'


# Second example using ogr and shapely:
from shapely.geometry import Polygon, Point, MultiPolygon
from osgeo import ogr, gdal

driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open("shapefile.shp", 0)

layer = dataset.GetLayer()
for index in xrange(layer.GetFeatureCount()):
    feature = layer.GetFeature(index)
    geometry = feature.GetGeometryRef()

polygon = Polygon(geometry)
print 'polygon points =', polygon  # this prints 'multipoint' + all the points fine

point = Point(x, y)
# point in polygon test
if polygon.contains(point):
    print 'inside'
else:
    print 'OUT'

The first example works fine with a single polygon at a time, but when I input a point within one of the shapes in my multipolygon shapefile it returns "out" even though it does fall inside one of the many parts.

For for the second example I get an error "object of type 'Geometry' has no len()" which I assume is because the geometry field can't be read as a normal, indexed list/array.

I additionally tried to replace the actual point in polygon code as suggested here to make sure it wasn't that part of the code that didn't work. And while that link's examples work fine with simple polygon shapefiles I can't get my complex multipolygon to test properly.

So I can't think of any other way to test whether a point falls within a multipolygon shapefile via python... Maybe there are other libraries out there I'm missing?

  • your second example looks like it might be coercing multipolygon to polygon? It might only be checking the point against the first part of the multipolygon. Try moving the point around to different parts and see if the check ever succeeds. – obrl_soil Aug 26 '16 at 22:44
  • @obrl_soil Thank you for your suggestion. However, the second example never works because of the error message I described above (object of type 'Geometry' has no len())" whether I try MultiPolygon (geometry) or simply Polygon(geometry). I have tried many points in the first example and only those within the main polygon work. Hope this clarification helps. – spartmar Aug 26 '16 at 22:48
  • Yeah, I think you need to replace polygon = Polygon(geometry) with some kind of try loop where it switches to polygon = MultiPolygon(geometry) if that error occurs. – obrl_soil Aug 27 '16 at 6:08
  • The problem in your first example is in the first loop. – xunilk Aug 27 '16 at 6:58
24

Shapefiles have no type MultiPolygon (type = Polygon), but they support them anyway (all rings are stored in one feature = list of polygons, look at Converting huge multipolygon to polygons)

The problem

enter image description here

If I open a MultiPolygon shapefile, the geometry is 'Polygon'

multipolys = fiona.open("multipol.shp")
multipolys.schema
{'geometry': 'Polygon', 'properties': OrderedDict([(u'id', 'int:10')])}
len(multipolys)
1

Solution 1 with Fiona

import fiona
from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon
multipol = fiona.open("multipol.shp")
multi= multipol.next() # only one feature in the shapefile
print multi
{'geometry': {'type': 'MultiPolygon', 'coordinates': [[[(-0.5275288092189501, 0.5569782330345711), (-0.11779769526248396, 0.29065300896286816), (-0.25608194622279135, 0.01920614596670933), (-0.709346991037132, -0.08834827144686286), (-0.8629961587708066, 0.18309859154929575), (-0.734955185659411, 0.39820742637644047), (-0.5275288092189501, 0.5569782330345711)]], [[(0.19974391805377723, 0.060179257362355965), (0.5480153649167734, 0.1293213828425096), (0.729833546734955, 0.03969270166453265), (0.8143405889884763, -0.13956466069142115), (0.701664532650448, -0.38540332906530095), (0.4763124199743918, -0.5006402048655569), (0.26888604353393086, -0.4238156209987196), (0.18950064020486557, -0.2291933418693981), (0.19974391805377723, 0.060179257362355965)]], [[(-0.3764404609475033, -0.295774647887324), (-0.11523687580025621, -0.3597951344430217), (-0.033290653008962945, -0.5800256081946222), (-0.11523687580025621, -0.7413572343149808), (-0.3072983354673495, -0.8591549295774648), (-0.58898847631242, -0.6927016645326505), (-0.6555697823303457, -0.4750320102432779), (-0.3764404609475033, -0.295774647887324)]]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', 1)])}

Fiona interprets the feature as a MultiPolygon and you can apply the solution presented in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (1)

points= ([pt for pt  in fiona.open("points.shp")])
for i, pt in enumerate(points):
    point = shape(pt['geometry'])
    if point.within(shape(multi['geometry'])):
         print i, shape(points[i]['geometry'])
1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.4993597951344431 -0.06017925736235585)
5 POINT (-0.3764404609475033 -0.4750320102432779)
6 POINT (-0.3098591549295775 -0.6312419974391805)

Solution 2 with pyshp (shapefile) and the geo_interface (GeoJSON like) protocol

This is a supplement to the answer of xulnik.

import shapefile
pts = shapefile.Reader("points.shp")
polys = shapefile.Reader("multipol.shp")
points = [pt.shape.__geo_interface__ for pt in pts.shapeRecords()]
multi = shape(polys.shapeRecords()[0].shape.__geo_interface__) # 1 polygon
print multi
MULTIPOLYGON (((-0.5275288092189501 0.5569782330345711, -0.117797695262484 0.2906530089628682, -0.2560819462227913 0.01920614596670933, -0.7093469910371319 -0.08834827144686286, -0.8629961587708066 0.1830985915492958, -0.734955185659411 0.3982074263764405, -0.5275288092189501 0.5569782330345711)), ((0.1997439180537772 0.06017925736235596, 0.5480153649167734 0.1293213828425096, 0.729833546734955 0.03969270166453265, 0.8143405889884763 -0.1395646606914211, 0.701664532650448 -0.3854033290653009, 0.4763124199743918 -0.5006402048655569, 0.2688860435339309 -0.4238156209987196, 0.1895006402048656 -0.2291933418693981, 0.1997439180537772 0.06017925736235596)), ((-0.3764404609475033 -0.295774647887324, -0.1152368758002562 -0.3597951344430217, -0.03329065300896294 -0.5800256081946222, -0.1152368758002562 -0.7413572343149808, -0.3072983354673495 -0.8591549295774648, -0.58898847631242 -0.6927016645326505, -0.6555697823303457 -0.4750320102432779, -0.3764404609475033 -0.295774647887324)))
for i, pt in enumerate(points):
    point = shape(pt)
    if point.within(multi): 
        print i, shape(points[i])
1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.4993597951344431 -0.06017925736235585)
5 POINT (-0.3764404609475033 -0.4750320102432779)
6 POINT (-0.3098591549295775 -0.6312419974391805)

Solution 3 with ogr and the geo_interface protocol (Python Geo_interface applications)

from osgeo import ogr
import json
def records(file):  
    # generator 
    reader = ogr.Open(file)
    layer = reader.GetLayer(0)
    for i in range(layer.GetFeatureCount()):
        feature = layer.GetFeature(i)
        yield json.loads(feature.ExportToJson())

points  = [pt for pt in records("point_multi_contains.shp")]
multipol = records("multipol.shp")
multi = multipol.next() # 1 feature
for i, pt in enumerate(points):
     point = shape(pt['geometry'])
     if point.within(shape(multi['geometry'])):
          print i, shape(points[i]['geometry'])

1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.499359795134443 -0.060179257362356)
5 POINT (-0.376440460947503 -0.475032010243278)
6 POINT (-0.309859154929577 -0.631241997439181)

Solution 4 with GeoPandas as in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (2)

import geopandas
point = geopandas.GeoDataFrame.from_file('points.shp') 
poly  = geopandas.GeoDataFrame.from_file('multipol.shp')
from geopandas.tools import sjoin
pointInPolys = sjoin(point, poly, how='left')
grouped = pointInPolys.groupby('index_right')
list(grouped)
[(0.0,      geometry                               id_left  index_right id_right  

1  POINT (-0.58898847631242 0.17797695262484)       None      0.0        1.0 
3  POINT (0.4993597951344431 -0.06017925736235585)  None      0.0        1.0
5  POINT (-0.3764404609475033 -0.4750320102432779)  None      0.0        1.0 
6  POINT (-0.3098591549295775 -0.6312419974391805)  None      0.0        1.0 ]
print grouped.groups
{0.0: [1, 3, 5, 6]} 

The points 1,3,5,6 falls within the boundaries of the MultiPolygon

  • Slightly old thread here, but how do you call multi = shape(polys.shapeRecords()[0].shape.__geo_interface__) in Solution 2? I cannot get a call to a shape() method from shapefile.py. I have even tried shapefile.Shape(); there is a class for it but it doesnt work. – pstatix Aug 31 '17 at 23:31
  • Further, where do you get the within() method from? – pstatix Aug 31 '17 at 23:45
  • 1
    from Shapely (from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon) – gene Sep 4 '17 at 15:04
  • I get this error using Solution 4 :File "C:\WinPython\python-3.6.5.amd64\lib\site-packages\geopandas\tools\sjoin.py", line 43, in sjoin if left_df.crs != right_df.crs: AttributeError: 'MultiPolygon' object has no attribute 'crs' – Aaron Bramson Apr 8 at 5:59
5

The problem in your first example is in this loop:

...
shpfilePoints = []
for shape in polygon:
    shpfilePoints = shape.points
...

It only appends the last feature points. I tried out my approach with this shapefile:

enter image description here

I modified your code to:

from shapely.geometry import Polygon, Point, MultiPolygon
import shapefile 

path = '/home/zeito/pyqgis_data/polygon8.shp'

polygon = shapefile.Reader(path) 

polygon = polygon.shapes() 

shpfilePoints = [ shape.points for shape in polygon ]

print shpfilePoints

polygons = shpfilePoints

for polygon in polygons:
    poly = Polygon(polygon)
    print poly

Above code was run at the Python Console of QGIS and the result was:

enter image description here

It works perfectly and now, you can check whether a point (x, y) falls within the boundaries of each feature.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.