How to use Rtree - 10 common examples

continue
            elif hole.intersects(polygons[place_id]):
                polygons[place_id] = polygons[place_id].union(hole)
                changed = True
            if changed:
                break
        if not changed:
            unassigned.append(hole)
            retries += 1
    if len(unassigned) == len(holes):
        # give up
        break
    holes = unassigned
print >>sys.stderr, "%d retried, %d unassigned." % (retries, len(unassigned))

hoodIndex = Rtree()

print >>sys.stderr, "Buffering polygons."
for place_id, polygon in polygons.items():
    if type(polygon) is Polygon:
        polygon = Polygon(polygon.exterior.coords)
    else:
        bits = []
        for p in polygon.geoms:
            if type(p) is Polygon:
                bits.append(Polygon(p.exterior.coords))
        polygon = MultiPolygon(bits)
    polygons[place_id] = polygon.buffer(0)
    hoodIndex.insert(place_id, polygons[place_id].bounds)

print >>sys.stderr, "Retconning blocks to shapes."
cur.execute("""select geom, geoid10 FROM tabblock10 tb WHERE statefp10 = %s AND countyfp10 = %s AND blockce10 NOT LIKE '0%%'""", (statefp10, countyfp10))

if do_cache:
      print >>sys.stderr, "Caching points, blocks, and names ..."
      pickle.dump((names, blocks, places), file(areaid + ".cache", "w"), -1)
    blocks = map(asShape, blocks)

points = []
place_list = set()
count = 0
for place_id, pts in places.items():
    count += 1
    print >>sys.stderr, "\rPreparing %d of %d places..." % (count, len(places)),
    for pt in pts:
        place_list.add((len(points), pt+pt, None))
        points.append((place_id, Point(pt)))
print >>sys.stderr, "Indexing...",
index = Rtree(place_list)
print >>sys.stderr, "Done."

def score_block(polygon):
    centroid = polygon.centroid
    #prepared = prep(polygon)
    score = {}
    outside_samples = 0
    for item in index.nearest((centroid.x, centroid.y), num_results=SAMPLE_SIZE):
        place_id, point = points[item]
        score.setdefault(place_id, 0.0)
        #if prepared.contains(point):
        #    score[place_id] += 1.0
        #else:
        score[place_id] += 1.0 / math.sqrt(max(polygon.distance(point)*SCALE_FACTOR, 1.0))
        outside_samples += 1
    return list(reversed(sorted((sc, place_id) for place_id, sc in score.items())))

count = 0
    while True:
      line = f.readline()
      if len(line) == 0:
        break
      id, lat, lon = line.strip().split(',')
      lat = float(lat)
      lon = float(lon)
      count += 1
      if count % 1000000 == 0:
        print('%d...' % count)
      yield int(id), (lat, lon, lat, lon), None
    

t0 = time.time()
p = rtree.index.Property()
p.filename = 'theindex'
p.overwrite = True
p.storage = rtree.index.RT_Disk
p.dimension = 2

# 17.5 sec
#p.variant = rtree.index.RT_Star
#p.fill_factor = 0.4

# 17.4 sec:
#p.variant = rtree.index.RT_Quadratic
#p.fill_factor = 0.4

# 17.5 sec:
#p.variant = rtree.index.RT_Linear
#p.fill_factor = 0.4

pvt_offset = 0

    mrf_dom = build_mrf_dom(tile_matrices, target_extents, tile_size, proj)
    with open(os.path.join(output_path, mrf_prefix) + '.mrf', 'w+') as f:
        f.write(mrf_dom.toprettyxml())

    spatial_dbs = []
    source_schemas = []
    # Dump contents of shapefile into a mutable rtree spatial database for faster searching.
    for input_file in input_file_path:
        print 'Processing ' + input_file
        with fiona.open(input_file) as shapefile:
            try:
                spatial_db = rtree.index.Index(
                    rtree_index_generator(list(shapefile), filter_list))
            except rtree.core.RTreeError as e:
                print 'ERROR -- problem importing feature data. If you have filters configured, the source dataset may have no features that pass. Err: {0}'.format(
                    e)
                sys.exit()
            spatial_dbs.append(spatial_db)
            source_schema = shapefile.schema['geometry']
            source_schemas.append(source_schema)
            if debug:
                print 'Points to process: ' + str(
                    spatial_db.count(spatial_db.bounds))

    # Build tilematrix pyramid from the bottom (highest zoom) up. We generate tiles left-right, top-bottom and write them
    # successively to the MRF.
    for i, tile_matrix in enumerate(reversed(tile_matrices)):
        z = len(tile_matrices) - i - 1

        for idx, spatial_db in enumerate(spatial_dbs):

if len(line) == 0:
        break
      id, lat, lon = line.strip().split(',')
      lat = float(lat)
      lon = float(lon)
      count += 1
      if count % 1000000 == 0:
        print('%d...' % count)
      yield int(id), (lat, lon, lat, lon), None
    

t0 = time.time()
p = rtree.index.Property()
p.filename = 'theindex'
p.overwrite = True
p.storage = rtree.index.RT_Disk
p.dimension = 2

# 17.5 sec
#p.variant = rtree.index.RT_Star
#p.fill_factor = 0.4

# 17.4 sec:
#p.variant = rtree.index.RT_Quadratic
#p.fill_factor = 0.4

# 17.5 sec:
#p.variant = rtree.index.RT_Linear
#p.fill_factor = 0.4

# 34 sec:
#p.leaf_capacity = 1000

def cherche_osm_buildings_proches(code_departement, code_commune, osm, transform_to_osm, transform_from_osm):
    """ Cherche a intégrer les nœuds "addr:housenumber" du fichier
        d'entrée osm avec les building extraits de la base OSM.
    """
    sys.stdout.write((u"Intégration avec les buidings proches présent dans la base OSM.\n").encode("utf-8"))
    sys.stdout.write((u"Chargement des buidings\n").encode("utf-8"))
    sys.stdout.flush();
    buildings_osm = get_osm_buildings_and_barrier_ways(code_departement, code_commune)
    for node in itertools.chain.from_iterable(
            [itervalues(o.nodes) for o in [osm, buildings_osm]]):
        if not hasattr(node,'xy'):
            node.position = transform_from_osm((float(node.attrs["lon"]), float(node.attrs["lat"])))
    # créé un index spatial de tous les ways:
    ways_index = rtree.index.Index()
    for way in itervalues(buildings_osm.ways):
        if way.nodes[0] == way.nodes[-1]:
            way.shape = Polygon([buildings_osm.nodes[id].position for id in way.nodes])
        else:
            way.shape = LineString([buildings_osm.nodes[id].position for id in way.nodes])
        ways_index.insert(way.id(), way.shape.bounds, way.id())
    sys.stdout.write((u"Recherche des buiding proches\n").encode("utf-8"))
    sys.stdout.flush();
    for node in osm.nodes.values():
        if "addr:housenumber" in node.tags:
            x,y = node.position
            search_bounds = [x - MAX_BUILDING_DISTANCE_METERS, y - MAX_BUILDING_DISTANCE_METERS,
                             x + MAX_BUILDING_DISTANCE_METERS, y + MAX_BUILDING_DISTANCE_METERS]
            near_ways = [buildings_osm.ways[e.object] for e in ways_index.intersection(search_bounds, objects=True)]
            if  hasattr(node, 'limite_parcelle') and node.limite_parcelle != None:
                    #and node.liimite_parcelle.distance(node.position) < MAX_BUILDING_DISTANCE_METERS:

#'landparcel_id': line[0],
                                #'mistral_function': line[1],
                                'toid': line['toid'],
                                #'household_id': line[4],
                                #'res_count': line[6],
                                #'lad': line[13],
                                'oa': line['oa'],
                            }
                        }
                        yield (i, geom_point.bounds, feature)
                        i += 1
    output = []

    try:
        # create index from generator (see http://toblerity.org/rtree/performance.html#use-stream-loading)
        idx = index.Index(premises())

        for n in idx.intersection((shape(exchange_area['geometry']).bounds), objects=True):
            point = n.object['representative_point']
            if prepared_area.contains(point):
                del n.object['representative_point']
                output.append(n.object)

    except:
        print('{} failed'.format(exchange_area['properties']['id']))

    return output

def compute_buildings_polygons_and_rtree(osm_data, tolerance):
    buildings_rtree = rtree.index.Index()
    osm_data.buildings_rtree = buildings_rtree 
    for way in itervalues(osm_data.ways):
        if way.isBuilding:
            if len(way.nodes) >= 3:
               way.polygon = Polygon([osm_data.nodes[i].position for i in way.nodes])
            else:
               way.polygon = LineString([osm_data.nodes[i].position for i in way.nodes])
            way.bbox = way.polygon.bounds
            way.tolerance_polygon = way.polygon.buffer(tolerance)
            buildings_rtree.insert(way.id(), way.bbox, way.textid())
    for rel in itervalues(osm_data.relations):
        if rel.isBuilding:
            exterior = None
            interiors = []
            for rtype,rref,rrole in rel.itermembers():
                if rtype == "way":

def createRtreeIndex(pf):
    idx = rtree.index.Index()
    for point in pf:
        point_id = int(point['id'])
        point_bounds = shapely.geometry.shape(point['geometry']).bounds
        idx.insert(point_id, point_bounds)
    print('Success -- created rtree')
    return idx

def compute_buildings_polygons_and_rtree(osm_data, tolerance):
    buildings_rtree = rtree.index.Index()
    osm_data.buildings_rtree = buildings_rtree 
    for way in osm_data.ways.itervalues():
        if way.isBuilding:
            if len(way.nodes) >= 3:
               way.polygon = Polygon([osm_data.nodes[i].position for i in way.nodes])
            else:
               way.polygon = LineString([osm_data.nodes[i].position for i in way.nodes])
            way.bbox = way.polygon.bounds
            way.tolerance_polygon = way.polygon.buffer(tolerance)
            buildings_rtree.insert(way.id(), way.bbox, way.textid())
    for rel in osm_data.relations.itervalues():
        if rel.isBuilding:
            exterior = None
            interiors = []
            for rtype,rref,rrole in rel.itermembers():
                if rtype == "way":