Holen Sie sich alle Eckpunkte eines Polygons mit OGR und Python

18

Ich habe ein kleines Problem mit der Python-OGR-API. Ich versuche, alle Koordinaten jedes Scheitelpunkts des äußeren Rings eines Polygons abzurufen.

Das habe ich bisher:

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

Ich habe gehört, dass Sie nur forüber die Region können, aber das gibt nur die Ringe im Polygon zurück, nicht die Knoten.

Jeder, der helfen kann.

python open-source-gis ogr Nathan W
quelle

Antworten:

15

Es hängt ein bisschen von Ihrem Dateiformat und Ihrer Geometrie ab, aber im Prinzip könnte die Fortsetzung so aussehen.

  for p in xrange(points):
        lon, lat, z = ring.GetPoint(p)
relet
quelle
Dies ist einer der folgenden Ausgänge: (1.8565347032449642e-313, 7.1913896573768921e-305, 6.3952653423603306e-305) Irgendeine Idee, was damit los ist?
Nathan W
Nicht genau, es ist ein Dreifach von Koordinaten, wenn auch ein bisschen winzig;) - wie sehen Ihre Eingabedaten und Projektionen aus? zB Was ogrinfo -alsagt das?
Relet
Das sieht für mich so aus, als würde es Floats als Double oder ähnlich interpretieren.
MerseyViking
5
Diese Zeile sollte lauten: lon, lat, z = ring.GetPoint(p)Welche für mich funktioniert.
MerseyViking
Danke MerseyViking, das hat es getan. Ich kann nicht glauben, dass ich darüber nachgesehen habe.
Nathan W
5

Ich bin auf dasselbe Problem gestoßen. Ich beendete die Verwendung der ExportToJson-Funktion in OGR und las dann den Json-String in ein Wörterbuch. Unter Verwendung meiner Daten und der Notation aus der ursprünglichen Frage sieht dies folgendermaßen aus:

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
quelle
4

Wenn Sie sich nur Shapefiles ansehen , können Sie auch pyshp verwenden .

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
quelle