Get shapefile polygon attribute value at a specific point using Python (e.g. via PyShp)

Question

I have a shapefile comprising multiple polygons. These polygons have attributes. I am trying to write a script in Python that takes a coordinate point, determines in which polygon it falls, and returns the attributes for that polygon.

So far, using PyShp I've got this, which prints off all the records.

import shapefile
sf = shapefile.Reader("England.shp")
records = sf.records()
print records

So I want to pass in a coordinates for Buckingham Palace, which is inside a polygon which has the attribute 'Area' value of 'London' and have it return 'London'

Can anyone point me in the right direction?

• gene: It is a spatial join and the same problem as in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc. You can use Pyshp as GeospatialPython.com shows, but also Fiona, shapely with or without a bounding spatial index as in the answer.

Answer 1

At the simplest level, you must loop through each polygon in your shapefile and check to see if they point is inside it. You can stop the loop when you find the right polygon. You can speed this process up a little bit by checking if the point is inside the extent of the polygon. By checking the extent first you are quickly narrowing the possibilies by quickly checking the 4 points of a rectangle.

The following function checks if a point is inside a polygon. This function alone will answer your question. If you use the pyshp shapeRecords() method, it will keep the geometry and the records associated so when you find the right polygon, you'll have easy access to the attributes:

def point_in_poly(x,y,poly):
 # check if point is a vertex
 if (x,y) in poly: return "IN"
 # check if point is on a boundary
 for i in range(len(poly)):
 p1 = None
 p2 = None
 if i==0:
 p1 = poly[0]
 p2 = poly[1]
 else:
 p1 = poly[i-1]
 p2 = poly[i]
 if p1[1] == p2[1] and p1[1] == y and x > min(p1[0], p2[0]) and x < max(p1[0], p2[0]):
 return "IN"
 n = len(poly)
 inside = False
 p1x,p1y = poly[0]
 for i in range(n+1):
 p2x,p2y = poly[i % n]
 if y > min(p1y,p2y):
 if y <= max(p1y,p2y):
 if x <= max(p1x,p2x):
 if p1y != p2y:
 xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
 if p1x == p2x or x <= xints:
 inside = not inside
 p1x,p1y = p2x,p2y
 if inside: return "IN"
 else: return "OUT"
# Test a vertex for inclusion
poligono = [(-33.416032,-70.593016), (-33.415370,-70.589604),
(-33.417340,-70.589046), (-33.417949,-70.592351),
(-33.416032,-70.593016)]
lat= -33.416032
lon= -70.593016
print point_in_poly(lat, lon, poligono)
# test a boundary point for inclusion
poly2 = [(1,1), (5,1), (5,5), (1,5), (1,1)]
x = 3
y = 1
print point_in_poly(x, y, poly2)