clip_grade10000.py

clip_grade.py

"""
clip_by_grade.py
────────────────
Recorta qualquer camada vetorial pela grade fixa (maps_grade10000),
adicionando o campo grade_id ao resultado.

Uso:
    python clip_by_grade.py

Requisitos:
    pip install geopandas shapely pyogrio
"""

import time
import geopandas as gpd
from shapely.geometry import MultiPolygon, Polygon


# ─── CONFIGURAÇÕES ────────────────────────────────────────────────────────────

GPKG_PATH = "/Users/marcellodebarrosfilho/Downloads/AMAZONIA_LEGAL_FINAL.gpkg"

GRADE_LAYER  = "maps_grade10000"
GRADE_ID_COL = "id"               # PK da grade → vira grade_id no resultado

# Camadas de entrada: cada entrada define a camada, as colunas desejadas
# e o nome da camada de saída no GeoPackage.
# "columns": None  →  mantém todas as colunas
CLIP_JOBS = [
    {
        "layer"   : "amazonia_legal",
        "columns" : None,                       # None = todas as colunas
        "output"  : "amazonia_legal_por_grade",
    },
    # Descomente e ajuste para adicionar mais camadas:
    # {
    #     "layer"   : "outra_camada",
    #     "columns" : ["nome", "area_ha"],       # lista de colunas desejadas
    #     "output"  : "outra_camada_por_grade",
    # },
]

# ─── CONSTANTES INTERNAS ──────────────────────────────────────────────────────

MIN_AREA    = 1e-10   # área mínima para filtrar artefatos de borda
GPKG_DRIVER = "GPKG"
ENGINE      = "pyogrio"


# ─── FUNÇÕES ──────────────────────────────────────────────────────────────────

def log(msg: str, indent: int = 0) -> None:
    prefix = "  " * indent
    print(f"{prefix}{msg}")


def separator(char: str = "─", width: int = 55) -> None:
    print(char * width)


def read_layer(gpkg: str, layer: str) -> gpd.GeoDataFrame:
    log(f"Lendo '{layer}'...", indent=1)
    gdf = gpd.read_file(gpkg, layer=layer, engine=ENGINE)
    log(f"{len(gdf)} feições | CRS: {gdf.crs}", indent=2)
    return gdf


def align_crs(gdf: gpd.GeoDataFrame, reference: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    if gdf.crs == reference.crs:
        return gdf
    log(f"Reprojetando {gdf.crs} → {reference.crs}", indent=2)
    return gdf.to_crs(reference.crs)


def select_columns(gdf: gpd.GeoDataFrame, columns: list | None) -> gpd.GeoDataFrame:
    if columns is None:
        return gdf
    missing = [c for c in columns if c not in gdf.columns]
    if missing:
        raise ValueError(f"Colunas não encontradas na camada: {missing}")
    return gdf[columns + ["geometry"]]


def to_multipolygon(geom):
    if geom.geom_type == "Polygon":
        return MultiPolygon([geom])
    return geom


def clean(gdf: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    gdf = gdf[~gdf.geometry.is_empty & gdf.geometry.notna()]
    gdf = gdf[gdf.geometry.area > MIN_AREA]
    gdf["geometry"] = gdf["geometry"].apply(to_multipolygon)
    return gdf


def intersect(
    input_gdf: gpd.GeoDataFrame,
    grade: gpd.GeoDataFrame,
) -> gpd.GeoDataFrame:
    log("Executando overlay (interseção)...", indent=1)
    result = gpd.overlay(
        input_gdf,
        grade,
        how="intersection",
        keep_geom_type=True,
        make_valid=True,
    )
    result["grade_id"] = result["grade_id"].astype("Int64")
    log(f"Feições geradas: {len(result)}", indent=2)
    return result


def save_layer(gdf: gpd.GeoDataFrame, gpkg: str, layer: str) -> None:
    log(f"Salvando '{layer}' no GeoPackage...", indent=1)
    gdf.to_file(gpkg, layer=layer, driver=GPKG_DRIVER, engine=ENGINE)
    log(f"{len(gdf)} feições salvas.", indent=2)


def prepare_grade(gpkg: str, layer: str, id_col: str) -> gpd.GeoDataFrame:
    gdf = read_layer(gpkg, layer)
    return gdf[[id_col, "geometry"]].rename(columns={id_col: "grade_id"})


def run_job(job: dict, grade: gpd.GeoDataFrame) -> None:
    layer   = job["layer"]
    columns = job.get("columns")
    output  = job["output"]

    separator()
    log(f"JOB: {layer} → {output}")
    separator()

    t0 = time.time()

    input_gdf = read_layer(GPKG_PATH, layer)
    input_gdf = align_crs(input_gdf, grade)
    input_gdf = select_columns(input_gdf, columns)

    result = intersect(input_gdf, grade)

    log("Limpando resultado...", indent=1)
    result = clean(result)
    log(f"Feições válidas: {len(result)}", indent=2)
    log(f"Campos: {[c for c in result.columns if c != 'geometry']}", indent=2)

    save_layer(result, GPKG_PATH, output)

    log(f"✅ Concluído em {time.time() - t0:.1f}s", indent=1)


# ─── MAIN ─────────────────────────────────────────────────────────────────────

def main() -> None:
    separator("=")
    log("clip_by_grade.py — Recorte por grade fixa")
    separator("=")

    log("\n[Carregando grade fixa...]")
    grade = prepare_grade(GPKG_PATH, GRADE_LAYER, GRADE_ID_COL)

    total_jobs = len(CLIP_JOBS)
    for i, job in enumerate(CLIP_JOBS, start=1):
        log(f"\nJob {i}/{total_jobs}")
        try:
            run_job(job, grade)
        except Exception as e:
            separator()
            log(f"❌ Erro no job '{job.get('layer')}': {e}")
            separator()

    separator("=")
    log(f"Todos os jobs finalizados. GeoPackage: {GPKG_PATH}")
    separator("=")


if __name__ == "__main__":
    main()