[gdal-dev] Reading S57 Extents from Python

[Mike Slinn]

Below is a heavily modified version of Frank's get_soundg.py sample 
program.  It dumps out rather more information than the original. A few 
questions:

   1. layer.GetExtent() does not seem to return anything useful (lines
      58-61).  Is GetExtent() not implemented?  Is there another way to
      efficiently obtain the extents of the visible entities on a layer?
   2. Do Points and Features need to be Destroy()'d (line 74-76)?  Those
      Python classes have no Destroy() methods.

Mike

------------------------------------------------------------------------

try:
    from osgeo import osr, ogr
except ImportError:
    import osr, ogr
import string, sys

def usage():
    print 'Usage: read_s57.py <s57file>\n'
    sys.exit(1)

def dumpFields(layerDefn):
	"""Dump all fields in layer"""	
	field_count = layerDefn.GetFieldCount()
	for fld_index in range(field_count):
		fd = layerDefn.GetFieldDefn(fld_index)
		print "      Field %s; type %d" % (fd.GetName(), fd.GetType())

def dumpS57(s57filename):
	ds = ogr.Open(s57filename)
	print str(ds.GetLayerCount()) + " layers found:"
	for index in range(ds.GetLayerCount()):
		layer = ds.GetLayer(index)
		extents = layer.GetExtent()
		layerDefn = layer.GetLayerDefn()
		print "   Layer %s is of type %d and has %d fields:" % \
			(layerDefn.GetName(), layerDefn.GetGeomType(), layerDefn.GetFieldCount())
		dumpFields(layerDefn)
		#print dir(extents) # doesn't show any methods ... is GetExtent() not implemented?
		#print "      Extents: (" + \
			#extents.MinX + ", " + extents.MinY + ") to (" + \
			#extents.MaxX + ", " + extents.MaxY + ") " #+ "; " + layer.GetLayerDefn()

	layer_soundg = ds.GetLayerByName('SOUNDG')
	print "=== Layer SOUNDG features ==="
	field_count = layer_soundg.GetLayerDefn().GetFieldCount()
	feat = layer_soundg.GetNextFeature()
	while feat is not None:
		multi_geom = feat.GetGeometryRef()
		for iPnt in range(multi_geom.GetGeometryCount()):
		    pnt = multi_geom.GetGeometryRef(iPnt)
		    print "   Lat, long: (%.3f, %.3f); depth: %.0f (units?)" % (pnt.GetX(), pnt.GetY(), pnt.GetZ())
		    for fld_index in range(field_count):
		        print "      Field #%d: %s" % (fld_index, feat.GetField(fld_index))
		    #if pnt is not None:
		    #    pnt.Destroy()  # causes error: about to Destroy unowned geometry
		#feat.Destroy()  # causes error: about to Destroy unowned geometry
		feat = layer_soundg.GetNextFeature()
	ds.Destroy()


if __name__ == '__main__':
	if len(sys.argv) != 2:
		usage()
	dumpS57(sys.argv[1])


-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20080416/e69681e7/attachment.html