Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I'm having a little trouble with the Python OGR API. What I am trying to do is get all the coordinates of each vertex of the outer ring of a polygon.

This is what I have so far:

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

I have heard that you can just for over the region but that only returns the rings in the polygon, not the nodes.

Anyone able to help.

share|improve this question
add comment

3 Answers 3

up vote 8 down vote accepted

It depends a bit on your file format and geometry, but in principle the continuation could look like this.

  for p in xrange(points):
        lon, lat, z = ring.GetPoint(p)
share|improve this answer
    
This is one of the outputs: (1.8565347032449642e-313, 7.1913896573768921e-305, 6.3952653423603306e-305) Any idea what's up with that? –  Nathan W Apr 5 '11 at 13:06
    
Not exactly, it's a triple of coordinates, if a bit tiny ;) - what does your input data and projection look like? e.g. What does ogrinfo -al say? –  relet Apr 5 '11 at 13:11
    
That looks to me like it's interpreting floats as doubles or similar. –  MerseyViking Apr 5 '11 at 14:56
5  
That line should read: lon, lat, z = ring.GetPoint(p) Which works for me. –  MerseyViking Apr 5 '11 at 15:47
    
Thanks MerseyViking, that did it..can't believe I looked over that. –  Nathan W Apr 5 '11 at 20:05
show 1 more comment

FYI - for a complete code example based on the original question, and a script you can use right away .. see https://github.com/spatialguru/NME/blob/master/ogr_explode.py

share|improve this answer
add comment

If you're looking just at shapefiles, you could also use 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
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.