22

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(polygon) point = Point(x, y) # point in polygon test if poly.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?

Improve this question
4
  • 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. Commented Aug 26, 2016 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. Commented Aug 26, 2016 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. Commented Aug 27, 2016 at 6:08
  • The problem in your first example is in the first loop. Commented Aug 27, 2016 at 6:58

3 Answers 3

Reset to default
35

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

Improve this answer
5
  • 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. Commented Aug 31, 2017 at 23:31
  • Further, where do you get the within() method from? Commented Aug 31, 2017 at 23:45
  • 1
    from Shapely (from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon) Commented Sep 4, 2017 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' Commented Apr 8, 2019 at 5:59
  • Reading this helped me to understand better gene's answer regarding geopandas's join usage: geopandas.org/en/stable/docs/user_guide/… Commented Mar 25, 2022 at 19:01
9

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.

Improve this answer
0
1

If you are trying to check a latitude, longitude point within a polygon, make sure you you have point object is created by the following:

from shapely.geometry.point import Point Point(LONGITUDE, LATITUDE) .. point.within(poly) # Returns true if the point within the polygon poly.contains(point) # Returns true if the polygon contains the point. Works as well.

Point takes longitude, then latitude in the argument. Not latitude first. You can call polygon_object.contains or polygon_object.within function to check if the point is within the shape.

Improve this answer
3
  • 1
    Any suggestion how I can do it fast and efficient when I have many points and many polygons? Looping over all is very costly. Commented Nov 8, 2020 at 6:05
  • You can use postgis to store your polygons and use one of the built-in functions of the server. Commented Nov 9, 2020 at 7:32
  • @mah65 If the distribution is very imbalanced (e.g. most points are in a few countries), ordering the polygons by order of most frequently seen can greatly increase performances. You can sample a number of random points from your data to find the most efficient ordering. Commented Apr 25, 2023 at 13:40

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.