Extracting IFC Wall Geometries to Shapely

# ifcopenshell>=0.8.0, shapely>=2.0.0, numpy>=1.24.0
import ifcopenshell
import ifcopenshell.geom
import ifcopenshell.util.shape
from shapely.geometry import Polygon, MultiPolygon
import shapely
import shapely.ops
import numpy as np
from typing import Literal


def extract_walls_to_shapely(
    ifc_path: str,
    projection_axis: Literal["X", "Y", "Z"] = "Z",
    min_area: float = 0.001,
    include_standard_case: bool = True,
) -> list[dict]:
    """
    Extract IfcWall (and optionally IfcWallStandardCase) geometries and
    convert them to validated Shapely Polygon / MultiPolygon objects.

    Returns a list of dicts with keys:
        id        (str)  GlobalId of the wall element
        name      (str)  Name attribute, may be None
        geometry  (Polygon | MultiPolygon | None)
        status    (str)  "success" | "no_faces" | "error"
        message   (str)  present only when status == "error"
    """
    model = ifcopenshell.open(ifc_path)

    # Collect wall types — IfcWallStandardCase is a subtype of IfcWall in
    # IFC2x3; IFC4 merges them, but both strings remain valid to query.
    entity_types = ["IfcWall"]
    if include_standard_case:
        entity_types.append("IfcWallStandardCase")

    walls = []
    seen_ids: set[str] = set()
    for etype in entity_types:
        for w in model.by_type(etype):
            if w.GlobalId not in seen_ids:
                walls.append(w)
                seen_ids.add(w.GlobalId)

    # Configure geometry kernel: world coordinates, mesh-only output.
    settings = ifcopenshell.geom.settings()
    settings.set(settings.USE_WORLD_COORDS, True)
    settings.set(settings.EXCLUDE_SOLIDS_AND_SURFACES, False)
    settings.set(settings.INCLUDE_CURVES, False)
    settings.set(settings.APPLY_DEFAULT_MATERIALS, False)

    axis_map = {"X": 0, "Y": 1, "Z": 2}
    proj_idx = axis_map[projection_axis.upper()]
    other_axes = [i for i in range(3) if i != proj_idx]

    results: list[dict] = []

    for wall in walls:
        try:
            shape = ifcopenshell.geom.create_shape(settings, wall)
            # get_vertices returns a flat list; reshape to (N, 3)
            verts = np.array(
                ifcopenshell.util.shape.get_vertices(shape)
            ).reshape(-1, 3)
            faces = ifcopenshell.util.shape.get_faces(shape)

            polygons: list[Polygon] = []
            for face in faces:
                # Project each triangle vertex to the 2D plane
                coords = [
                    (float(verts[i][other_axes[0]]),
                     float(verts[i][other_axes[1]]))
                    for i in face
                ]
                try:
                    poly = Polygon(coords)
                    if poly.is_valid and poly.area > min_area:
                        polygons.append(poly)
                except Exception:
                    # Degenerate triangle (e.g. all vertices collinear)
                    continue

            if not polygons:
                results.append({
                    "id": wall.GlobalId,
                    "name": getattr(wall, "Name", None),
                    "geometry": None,
                    "status": "no_faces",
                })
                continue

            # Merge adjacent triangles into a coherent wall footprint
            merged = shapely.ops.unary_union(polygons)
            # Repair self-intersections produced by projection
            valid_geom = shapely.make_valid(merged)

            results.append({
                "id": wall.GlobalId,
                "name": getattr(wall, "Name", None),
                "geometry": valid_geom,
                "status": "success",
            })

        except Exception as exc:
            results.append({
                "id": wall.GlobalId,
                "name": getattr(wall, "Name", None),
                "geometry": None,
                "status": "error",
                "message": str(exc),
            })

    return results


# ── Usage ──────────────────────────────────────────────────────────────────
if __name__ == "__main__":
    import json

    walls = extract_walls_to_shapely("model.ifc", projection_axis="Z")

    success = [w for w in walls if w["status"] == "success"]
    errors  = [w for w in walls if w["status"] == "error"]

    print(f"Extracted {len(success)} walls, {len(errors)} errors")

    # Export to GeoJSON for PostGIS / QGIS ingestion
    features = []
    for w in success:
        features.append({
            "type": "Feature",
            "properties": {"globalId": w["id"], "name": w["name"]},
            "geometry": shapely.to_geojson(w["geometry"])
                        if hasattr(shapely, "to_geojson")
                        else w["geometry"].__geo_interface__,
        })

    with open("walls.geojson", "w") as f:
        json.dump({"type": "FeatureCollection", "features": features}, f)

from ifcopenshell.util.element import get_psets

def has_body_representation(wall) -> bool:
    if not wall.Representation:
        return False
    return any(
        rep.RepresentationIdentifier == "Body"
        for rep in wall.Representation.Representations
    )

walls = [w for w in model.by_type("IfcWall") if has_body_representation(w)]

def face_normal(verts, face):
    """Compute the unit normal of a triangular face."""
    v0, v1, v2 = [verts[i] for i in face]
    edge1 = v1 - v0
    edge2 = v2 - v0
    n = np.cross(edge1, edge2)
    length = np.linalg.norm(n)
    return n / length if length > 1e-10 else np.zeros(3)

# Keep only faces with a significant Z component (near-horizontal)
horizontal_faces = [
    f for f in faces
    if abs(face_normal(verts, f)[2]) > 0.5
]

from ifcopenshell.util.unit import calculate_unit_assignment

unit_scale = calculate_unit_assignment(model, "LENGTHUNIT")
# unit_scale == 0.001 when the file uses millimetres
verts = verts * unit_scale

from shapely.geometry import GeometryCollection

def extract_polygons(geom):
    if geom.is_empty:
        return None
    if isinstance(geom, (Polygon, MultiPolygon)):
        return geom
    if isinstance(geom, GeometryCollection):
        polys = [g for g in geom.geoms if isinstance(g, (Polygon, MultiPolygon))]
        return shapely.ops.unary_union(polys) if polys else None
    return None

BATCH_SIZE = 500

for i in range(0, len(walls), BATCH_SIZE):
    batch = walls[i : i + BATCH_SIZE]
    batch_results = [process_wall(w, settings, other_axes, min_area) for w in batch]
    # write batch_results to database here, then let the batch go out of scope