¿Cómo almacenar un archivo de forma vectorial usando ogr python?

10

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()
usuario1765941
fuente

Respuestas:

2

Dos opciones que se me ocurren son:

  1. Calcule el grado equivalente de 500 kilómetros. Luego puede ingresarlo la función Buffer (). Sin embargo, debe tener cuidado ya que un grado no tiene un equivalente métrico constante. Depende de en qué latitud estés. Puede consultar la fórmula de Haversine si desea seguir esta ruta.

  2. Otra opción sería volver a proyectar el archivo de forma a UTM. De esa manera, puedes usar 500 kilómetros directamente. Sin embargo, tendrá que encontrar la zona UTM para su área de interés. (Que debería ser la zona UTM 32S para Angola si no me equivoco)

RK
fuente
3
No existe un "grado equivalente" de 500 kilómetros, ni de ninguna otra distancia, excepto aproximadamente para distancias pequeñas cerca del ecuador: eso se debe a que la relación entre distancias y grados cambia con la demora y la latitud. Por lo tanto, la primera opción generalmente no funcionará correctamente.
whuber
0
  1. Si desea hacer un búfer con grados, considere que los grados tienen una distancia bastante diferente dependiendo de la dirección cuando no esté cerca del ecuador. La latitud se mantiene igual, pero la longitud de un grado es mucho menor en latitudes altas. A continuación se muestra mi tabla de cuadrados de 500 km en grados en diferentes latitudes. Supongo que para Angola el valor de 4.4 puede ser una buena suposición si no necesita una alta precisión.
  2. Puede reproyectar objetos en python ogr (hay una función Transformar para ello) durante la lectura, luego no es necesario realizar conversiones de archivos de forma.
500 km en lat 0.0 es 4.491576420597608 x 4.486983030705042 grados
500 km en lat 10.0 es 4.491576420597608 x 4.389054945583991 grados
500 km en lat 20.0 es 4.491576420597608 x 4.16093408959923 grados
500 km en lat 30.0 es 4.491576420597608 x 3.8117296267699388 grados
500 km en lat 40.0 es 4.491576420597608 x 3.3535548944407267 grados
500 km en lat 50.0 es 4.491576420597608 x 2.8010165014556634 grados
500 km en lat 60.0 es 4.491576420597608 x 2.170722673038327 grados
500 km en lat 70.0 es 4.491576420597608 x 1.4808232946314916 grados
500 km en lat 80.0 es 4.491576420597608 x 0.7505852760718597 grados
500 km en lat 84.0 es 4.491576420597608 x 0.4516575041056399 grados
JaakL
fuente
44
El usuario @Dave X señala que esta tabla es errónea: una distancia fija abarca un mayor número de grados en latitudes más altas, no menor. Parece que podría haberse construido haciendo una multiplicación donde se requiere una división. Aun así, eso no explica completamente las discrepancias: quedan errores del orden de varios por ciento. Exactamente, ¿cómo calculaste estos números?
whuber