North Beach Buildings and Topography with BLOSM + USGS

import bpy
import mathutils
import bmesh
from mathutils.bvhtree import BVHTree
 
# Get terrain and building collection
terrain = bpy.data.objects["Terrain"]
buildings = bpy.data.collections["all_buildings"].objects
 
depsgraph = bpy.context.evaluated_depsgraph_get()
terrain_eval = terrain.evaluated_get(depsgraph)
terrain_mesh = terrain_eval.to_mesh()
 
# Create bmesh and BVH tree from terrain
bm = bmesh.new()
bm.from_mesh(terrain_mesh)
bm.transform(terrain.matrix_world)  # Important: account for object transforms
bm.normal_update()
bvh = BVHTree.FromBMesh(bm)
 
# Raycast each object to terrain
for obj in buildings:
    if obj.type != 'MESH':
        continue
 
    # Get world-space bounding box corners
    bbox_world = [obj.matrix_world @ mathutils.Vector(corner) for corner in obj.bound_box]
    min_point = min(bbox_world, key=lambda v: v.z)
 
    ray_origin = min_point + mathutils.Vector((0, 0, 1000))
    ray_direction = mathutils.Vector((0, 0, -1))
 
    hit = bvh.ray_cast(ray_origin, ray_direction)
 
 
 
    hit = bvh.ray_cast(ray_origin, ray_direction)
 
    if hit is not None and hit[0] is not None:
        location, normal, index, distance = hit
        delta_z = location.z - min_point.z
        obj.location.z += delta_z
    else:
        print(f"[WARN] Raycast missed for object '{obj.name}' at XY = ({min_point.x:.2f}, {min_point.y:.2f})")
 
# Cleanup
bm.free()
terrain_eval.to_mesh_clear()