Obtener todos los vértices de un polígono usando OGR y Python

Question

Estoy teniendo un pequeño problema con la API OGR de Python. Lo que estoy tratando de hacer es obtener todas las coordenadas de cada vértice del anillo exterior de un polígono.

Esto es lo que tengo hasta ahora:

import osgeo.ogr
import glob

path = "/home/woo/maps/"
out = path + 'output.txt'

file = open(out,'w')
for filename in glob.glob(path + "*.shp"):
    ds = osgeo.ogr.Open(filename)
    layer1 = ds.GetLayer(0)
    print layer1.GetExtent()    
    for feat in layer1:
        geom = feat.GetGeometryRef()
        ring = geom.GetGeometryRef(0)
        points = ring.GetPointCount()
        #Not sure what to do here


file.close()

He oído que sólo se puede %-% for la región, pero que sólo devuelve los anillos en el polígono, no los nodos.

Cualquiera que pueda ayudar.

Answer 1

Depende un poco de su formato de archivo y geometría, pero en principio la continuación podría verse así.

  for p in xrange(points):
        lon, lat, z = ring.GetPoint(p)

Answer 2

Para su información, para obtener un ejemplo de código completo basado en la pregunta original y un script que puede usar de inmediato, consulte https://github.com/spatialguru/NME/blob/master/ogr_explode.py

Answer 3

Acabo de tener el mismo problema. Terminé usando la función ExportToJson en ogr y luego leyendo la cadena Json en un diccionario. Usando mis datos y la notación de la pregunta original, esto se ve como:

import json
...
ring_dict = json.loads(ring.ExportToJson())
ring_dict

{'coordinates': [[-4.94237, 55.725449],
  [-4.941922, 55.725585],
  [-4.9420024, 55.7252119],
  [-4.9422001, 55.7250997],
  [-4.9423197, 55.7251789],
  [-4.9425472, 55.7253089],
  [-4.94237, 55.725449]],
 'type': 'LineString'}

Answer 4

Si solo estás viendo shapefiles, también puedes usar pyshp.

import shapefile
sf = shapefile.Reader("shapefiles/blockgroups")
shapes = sf.shapes()
for shape in shapes:
  for vertex in shape.points:
    #do something with the vertex