Uzyskaj wszystkie wierzchołki wielokąta za pomocą OGR i Pythona

18

Mam trochę problemów z interfejsem API Python OGR. Staram się uzyskać wszystkie współrzędne każdego wierzchołka zewnętrznego pierścienia wielokąta.

Oto co mam do tej pory:

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

Słyszałem, że możesz przejść forprzez region, ale to zwraca tylko pierścienie w wielokącie, a nie węzły.

Każdy, kto może pomóc.

python open-source-gis ogr

— Nathan W.
źródło

15

Zależy to trochę od formatu pliku i geometrii, ale w zasadzie kontynuacja może wyglądać tak.

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

— relet
źródło

Jest to jedno z wyników: (1.8565347032449642e-313, 7.1913896573768921e-305, 6.3952653423603306e-305) Masz pojęcie, co z tym wspólnego?

— Nathan W

Nie do końca, to potrójne współrzędne, choć trochę małe;) - jak wyglądają twoje dane wejściowe i projekcja? np. co ogrinfo -almówi?

— relet

Wygląda mi to na interpretację liczb zmiennoprzecinkowych jako podwójnych lub podobnych.

— MerseyViking,

5

Ta linijka powinna brzmieć: lon, lat, z = ring.GetPoint(p)Która działa dla mnie.

— MerseyViking,

Dzięki MerseyViking, że to zrobił ... nie mogę uwierzyć, że to obejrzałem.

— Nathan W

6

FYI - kompletny przykład kodu oparty na oryginalnym pytaniu oraz skrypt, którego można użyć od razu .. patrz https://github.com/spatialguru/NME/blob/master/ogr_explode.py

— 1tylermitchell
źródło

5

Właśnie natrafiłem na ten sam problem. Skończyłem używać funkcji ExportToJson w ogrodzie, a następnie czytałem ciąg Json w słowniku. Używając moich danych i zapisu z pierwotnego pytania, wygląda to tak:

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.
źródło

4

Jeśli patrzysz tylko na pliki kształtów, możesz także użyć 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
źródło