import shapefile
from shapely.geometry import Point # Point class
from shapely.geometry import shape # shape() is a function to convert geo objects through the interface
pt = (-97.759615,30.258773) # an x,y tuple
shp = shapefile.Reader('/home/af/Downloads/cb_2016_us_ua10_500k/cb_2016_us_ua10_500k.shp') #open the shapefile
all_shapes = shp.shapes() # get all the polygons
all_records = shp.records()
def is_urban(pt):
result = False
for i in range(len(all_shapes)):
boundary = all_shapes[i] # get a boundary polygon
#name = all_records[i][3] + ', ' + all_records[i][4] # get the second field of the corresponding record
if Point(pt).within(shape(boundary)): # make a point and see if it's in the polygon
result = True
return result
result = is_urban(pt)
我最终使用了从下载的shapely和shapefile
https://www.census.gov/geo/maps-data/data/cbf/cbf_ua.html
,它有美国的城市区域,所以如果一个点不在这些区域内,它就是农村。
我测试了它,它达到了我的预期。