Unions of polygonal boxes in GDAL

boxunion.py

#!/usr/bin/python

from __future__ import division
from numpy import *
from osgeo import ogr

def box(scale, offset):
  sample = offset + array([[-1,-1],[-1,1],[1,1],[1,-1]])
  sample = scale * sample
  ring = ogr.Geometry(ogr.wkbLinearRing)  # create Polygon geometry
  for (x, y) in sample:		      
    ring.AddPoint(float(x), float(y))	      # add each point
  ring.CloseRings()
  poly = ogr.Geometry(ogr.wkbPolygon)
  poly.AddGeometry(ring)
  return poly   # convex hull produces a POLYGON (ogr.wkbPolygon)


if __name__ == "__main__":
  poly1 = box(20, 0)
  poly2 = box(20, 1.5)
  poly3 = box(20, -1)
  poly  = poly1.Union(poly2)
  poly  = poly.Union(poly3)

  ring = poly.GetGeometryRef(0)    # with a single LINEARRING (ogr.wkbLinearRing)

  # produce postscript output
  print "100 500 translate 2 2 scale 0 0 moveto"
  print "/tick {moveto 0 2 rlineto 0 -4 rlineto 0 2 rlineto"
  print "              2 0 rlineto -4 0 rlineto 2 0 rlineto} def"
  
  for ref in range(poly.GetGeometryCount()):
    ring = poly.GetGeometryRef(ref)
    if "POLYGON" == ring.GetGeometryName():
      ring = ring.GetGeometryRef(0)
    print "stroke"
    x, y, _ = ring.GetPoint(0)
    print x, y, "moveto"
    points = ring.GetPointCount()
    for p in range(points)[1:]:
      x, y, _ = ring.GetPoint(p)
      print x, y, "lineto"
    print "closepath stroke"
  print "showpage"