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

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
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?

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
p2 = poly
else:
p1 = poly[i-1]
p2 = poly[i]
if p1 == p2 and p1 == y and x > min(p1, p2) and x < max(p1, p2):
return "IN"

n = len(poly)
inside = False

p1x,p1y = poly
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)