Failing function to take union of polygons and write to a vector

union.py

"""Module for preprocessing operations."""

# Including all import statements from the original file, in case that matters
# Original file available (raw) at http://adept.invest-natcap.googlecode.com/hg/adept/preprocessing.py
import os
import logging
import tempfile
import shutil
import math
import json
import glob
import urllib
import zipfile


from osgeo import ogr
from osgeo import gdal
from osgeo import osr
import shapely
import shapely.wkb
import shapely.ops
import shapely.speedups
import shapely.geos
import shapely.prepared
import shapely.validation
import shapely.geometry
import rtree

LOGGER = logging.getLogger('adept.preprocessing')

def rm_shapefile(uri):
    """Delete all files associated with the user-defined ESRI Shapefile.

    uri - a URI to an ESRI shapefile on disk.

    Returns nothing."""
    shapefile_base = os.path.splitext(uri)[0]

    for filename in glob.glob(shapefile_base + '.*'):
        LOGGER.debug('Removing %s', filename)
        os.remove(filename)

def build_shapely_polygon(ogr_feature, prep=False, fix=False):
    geometry = ogr_feature.GetGeometryRef()
    try:
        polygon = shapely.wkb.loads(geometry.ExportToWkb())
    except shapely.geos.ReadingError:
        LOGGER.debug('Attempting to close geometry rings')
        # If the geometry does not form a closed circle, try to connect the
        # rings with the OGR function.
        geometry.CloseRings()
        polygon = shapely.wkb.loads(geometry.ExportToWkb())
        LOGGER.debug('Geometry fixed')

    if fix:
        polygon = polygon.buffer(0)
    if prep:
        polygon = shapely.prepared.prep(polygon)
    return polygon

def union_of_vectors(input_vector_uris, out_vector_uri):
    """Find the union of all features in the provided input vectors and write
    them out to a single geometry in out_vector_uri.  All vectors are assumed
    to be in the same projection.

        input_vector_uris - a list of python string URIs to OGR vectors on disk
        out_vector_uri - a python string URI.  The output OGR vector will be
            written to this location.

    Returns nothing."""

    source_vector_uri = input_vector_uris[0]
    LOGGER.debug('Opening input vector[0] %s', os.path.abspath(source_vector_uri))
    in_vector = ogr.Open(source_vector_uri)
    in_layer = in_vector.GetLayer()
    in_layer_srs = in_layer.GetSpatialRef()

    if os.path.exists(out_vector_uri):
        rm_shapefile(out_vector_uri)

    LOGGER.debug('Creating a new vector at %s', out_vector_uri)
    out_driver = ogr.GetDriverByName('ESRI Shapefile')
    out_vector = out_driver.CreateDataSource(out_vector_uri)
    # Cast to str here because the ESRI shapefile expects it.
    layer_name = str(os.path.basename(os.path.splitext(out_vector_uri)[0]))
    new_srs = osr.SpatialReference(wkt=in_layer_srs.ExportToWkt())
    out_layer = out_vector.CreateLayer(layer_name, new_srs, ogr.wkbPolygon)
    ogr.DataSource.__swig_destroy__(in_vector)

    LOGGER.debug('Building shapely geometries')
    all_features = []

    for in_vector in input_vector_uris:
        LOGGER.debug('Processing %s', in_vector)
        vector = ogr.Open(in_vector)
        layer = vector.GetLayer()

        for feature in layer:
            polygon = build_shapely_polygon(feature)
            all_features.append(polygon)

        # clean up the input vector.
        ogr.DataSource.__swig_destroy__(vector)

    sorted_features = sorted(all_features, key=lambda x: x.area, reverse=True)

    biggest_feature = sorted_features[0]
    features_to_test = sorted_features[1:]

    union_features = [biggest_feature]
    for feature in features_to_test:
        if biggest_feature.contains(feature):
            continue
        else:
            union_features.append(feature)

    # USUALLY FAILS HERE
    LOGGER.debug('Taking union of %s features', len(union_features))
    biggest_feature = shapely.ops.cascaded_union(union_features)

    union_polygon = biggest_feature

    LOGGER.debug('Creating new feature from empty definition')
    new_feature = ogr.Feature(ogr.FeatureDefn())
    new_geometry = ogr.CreateGeometryFromWkb(union_polygon.wkb)
    LOGGER.debug('Setting geometry from union_polygon with len(wkt)=%s',
        len(union_polygon.wkt))
    new_feature.SetGeometry(new_geometry)

    LOGGER.debug('Creating new feature')
    out_layer.CreateFeature(new_feature)

    LOGGER.debug('Cleaning up')
    new_feature = None
    out_layer = None
    out_vector = None
    out_driver = None
    LOGGER.debug('Finished cleanup of new vector')
    
if __name__ == '__main__':
    if not os.path.exists('vectors.zip'):
        urllib.urlretrieve ("http://ncp-yamato.stanford.edu/~jdouglass/shapely_sample_vectors.zip", "vectors.zip")
        with zipfile.ZipFile('vectors.zip') as archive:
            archive.extractall()
    
    out_vector = 'out.shp'
    in_vectors = glob.glob('sample_vectors/*.shp')
    union_of_vectors(in_vectors, out_vector)

Running this with python union.py causes the following crash:

Assertion failed: (!"should never be reached"), function itemsTree, file AbstractSTRtree.cpp, line 371.
Abort trap: 6