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

18

Tengo un pequeño problema con la API Python OGR. Lo que intento 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()

Escuché que puedes forsobre la región, pero eso solo devuelve los anillos en el polígono, no los nodos.

Cualquiera capaz de ayudar.

python open-source-gis ogr Nathan W
fuente

Respuestas:

15

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

  for p in xrange(points):
        lon, lat, z = ring.GetPoint(p)
realquilar
fuente
Esta es una de las salidas: (1.8565347032449642e-313, 7.1913896573768921e-305, 6.3952653423603306e-305) ¿Alguna idea de qué pasa con eso?
Nathan W
No exactamente, es un triple de coordenadas, aunque un poco pequeño;) - ¿cómo son sus datos de entrada y su proyección? Por ejemplo, ¿qué ogrinfo -aldice?
Relet
Me parece que está interpretando flotadores como dobles o similares.
MerseyViking 05 de
55
Esa línea debería leer: lon, lat, z = ring.GetPoint(p)Lo que funciona para mí.
MerseyViking 05 de
Gracias MerseyViking, eso lo hizo ... no puedo creer que haya revisado eso.
Nathan W
5

Me encontré con 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 así:

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'}
David M
fuente
4

Si solo estás mirando archivos de forma, 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
Marc Pfister
fuente