Estoy tratando de aprender a usar ogr en python usando el país y lugares poblados de conjuntos de datos de http://www.naturalearthdata.com/downloads/50m-cultural-vectors/. Estoy tratando de usar filtros y buffers para encontrar puntos (ne_50m_populated_places.shp) dentro de un búfer especificado de nombre de país (filtrado de la clase de entidad de ADMINISTRACIÓN en ne_50m_admin_0_countries.shp). El problema parece ser que no entiendo lo de las unidades utilizadas para el buffer(). En el script simplemente he utilizado un valor arbitrario de 10 para probar si el script funciona. La secuencia de comandos se ejecuta, pero devuelve lugares poblados de todo el Caribe la región de país llamado 'Angola'. Idealmente, quiero ser capaz de especificar una distancia intermedia, digamos 500 km, pero no se puede averiguar cómo hacerlo ya que en mi entendimiento es de búfer() es el uso de las unidades de los países.shp que será en wgs84 (lat/long formato. Asesoramiento sobre el método para lograr esto sería muy apreciada.
# 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()
