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?

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)