Estoy tratando de aprender cómo usar ogr en python usando los conjuntos de datos de países y lugares poblados de http://www.naturalearthdata.com/downloads/50m-cultural-vectors/. Estoy intentando usar filtros y zonas de influencia para encontrar puntos (ne_50m_populated_places.shp) dentro de una zona de influencia específica de un país con nombre (filtrado de la clase de entidad ADMIN en ne_50m_admin_0_countries.shp). El problema parece ser que no entiendo qué unidades usar para buffer (). En el script simplemente he usado un valor arbitrario de 10 para probar si el script funciona. El script se ejecuta pero devuelve lugares poblados de toda la región del Caribe para el país llamado 'Angola'. Idealmente, quiero poder especificar una distancia de búfer, digamos 500 km, pero no puedo entender cómo hacerlo, ya que entiendo que buffer () está usando las unidades de países.shp que estará en formato wgs84 lat / long . Se agradecería mucho el asesoramiento sobre el método para lograrlo.
# import modules
import ogr, os, sys
## data source
os.chdir('C:/data/naturalearth/50m_cultural')
# get the shapefile driver
driver = ogr.GetDriverByName('ESRI Shapefile')
# open ne_50m_admin_0_countries.shp and get the layer
admin = driver.Open('ne_50m_admin_0_countries.shp')
if admin is None:
print 'Could not open ne_50m_admin_0_countries.shp'
sys.exit(1)
adminLayer = admin.GetLayer()
# open ne_50m_populated_places.shp and get the layer
pop = driver.Open('ne_50m_populated_places.shp')
if pop is None:
print 'could not open ne_50m_populated_places.shp'
sys.exit(1)
popLayer = pop.GetLayer()
# use an attribute filter to restrict ne_50m_admin_0_countries.shp to "Angola"
adminLayer.SetAttributeFilter("ADMIN = ANGOLA")
# get the Angola geometry and buffer it by 10 units
adminFeature = adminLayer.GetFeature(0)
adminGeom = adminFeature.GetGeometryRef()
bufferGeom = adminGeom.Buffer(10)
# use bufferGeom as a spatial filter on ne_50m_populated_places.shp to get all places
# within 10 units of Angola
popLayer.SetSpatialFilter(bufferGeom)
# loop through the remaining features in ne_50m_populated_places.shp and print their
# id values
popFeature = popLayer.GetNextFeature()
while popFeature:
print popFeature.GetField('NAME')
popFeature.Destroy()
popFeature = popLayer.GetNextFeature()
# close the shapefiles
admin.Destroy()
pop.Destroy()
fuente