Realiza o clip das camadas vetoriais em função de uma grade

grade_clip.py

import os
import geopandas as gpd
from pathlib import Path
from shapely.geometry import MultiPolygon

# ============================================================================
# CONFIGURAÇÕES
# ============================================================================
# Caminhos de entrada
CAMINHO_INPUT = '/Users/marcellodebarrosfilho/code/curso_geopandas/dados/pb.gpkg'
CAMADA_INPUT = 'municipios'

CAMINHO_GRID = '/Users/marcellodebarrosfilho/Downloads/pluviosidade_media/Pluviosidade_Media_Anual_1991_2020.gpkg'
CAMADA_GRID = 'maps_grade10000'

# Caminho de saída (usando Path para compatibilidade multiplataforma)
PASTA_SAIDA = Path('/Users/marcellodebarrosfilho/code/curso_geopandas/dados/')
NOME_ARQUIVO = 'pb.gpkg'
NOME_CAMADA = 'municipios_grade'

# ============================================================================
# 1. CARREGAR DADOS
# ============================================================================
print("Carregando dados...")
input_layer = gpd.read_file(CAMINHO_INPUT, layer=CAMADA_INPUT)
print(f"  ✓ Camada de entrada carregada: {len(input_layer)} features")

# Carregar apenas colunas necessárias do grid
overlay_grid = gpd.read_file(
    CAMINHO_GRID,
    layer=CAMADA_GRID,
    columns=['id', 'geometry']
)
print(f"  ✓ Grid carregado: {len(overlay_grid)} features")

# ============================================================================
# 2. GARANTIR MESMO CRS
# ============================================================================
if input_layer.crs != overlay_grid.crs:
    print(f"Reprojetando de {input_layer.crs} para {overlay_grid.crs}...")
    input_layer = input_layer.to_crs(overlay_grid.crs)
    print("  ✓ Reprojeção concluída")

# ============================================================================
# 3. FILTRO ESPACIAL (OTIMIZAÇÃO)
# ============================================================================
print("\nAplicando filtro espacial...")
bbox_input = input_layer.total_bounds
overlay_grid_filtered = overlay_grid.cx[bbox_input[0]:bbox_input[2], 
                                        bbox_input[1]:bbox_input[3]]

# Filtro adicional com intersects para maior precisão
overlay_grid_filtered = overlay_grid_filtered[
    overlay_grid_filtered.intersects(input_layer.union_all())
]

print(f"  Grid original: {len(overlay_grid)} features")
print(f"  Grid filtrado: {len(overlay_grid_filtered)} features")
print(f"  Redução: {100 * (1 - len(overlay_grid_filtered)/len(overlay_grid)):.1f}%")

# ============================================================================
# 4. OVERLAY (INTERSEÇÃO)
# ============================================================================
print("\nExecutando overlay (interseção)...")
clipped = gpd.overlay(
    input_layer, 
    overlay_grid_filtered[['id', 'geometry']], 
    how='intersection',
    keep_geom_type=False
)
print(f"  ✓ Overlay concluído: {len(clipped)} features resultantes")

# ============================================================================
# 5. AJUSTES FINAIS
# ============================================================================
# Renomear coluna id do grid
clipped = clipped.rename(columns={'id': 'grade_id'})

# Validar geometrias
clipped = clipped[clipped.geometry.is_valid]
print(f"  ✓ Dados preparados: {len(clipped)} features válidas")

# Converter todas as geometrias para MultiPolygon
def to_multipolygon(geom):
    """Converte Polygon para MultiPolygon"""
    if geom.geom_type == 'Polygon':
        return MultiPolygon([geom])
    elif geom.geom_type == 'MultiPolygon':
        return geom
    else:
        # Para outros tipos (GeometryCollection, etc), tenta extrair polígonos
        return MultiPolygon([g for g in geom.geoms if g.geom_type == 'Polygon'])

print("\nConvertendo geometrias para MultiPolygon...")
clipped['geometry'] = clipped['geometry'].apply(to_multipolygon)
print(f"  ✓ Conversão concluída")
print(f"  📊 Tipos de geometria: {clipped.geometry.geom_type.value_counts().to_dict()}")

# ============================================================================
# 6. SALVAR RESULTADO
# ============================================================================
print(f"\nSalvando resultado...")

# Criar pasta de saída se não existir
try:
    PASTA_SAIDA.mkdir(parents=True, exist_ok=True)
    arquivo_saida = PASTA_SAIDA / NOME_ARQUIVO
    
    # Salvar com opções otimizadas
    clipped.to_file(
        arquivo_saida, 
        layer=NOME_CAMADA,
        driver='GPKG',
        engine='fiona'
    )
    
    print(f"  ✓ Arquivo salvo com sucesso!")
    print(f"  📁 Localização: {arquivo_saida}")
    print(f"  📊 Layer: {NOME_CAMADA}")
    print(f"  📐 Features: {len(clipped)}")
    print(f"  🗺️  CRS: {clipped.crs}")
    print(f"  🔷 Tipo de geometria: MultiPolygon")
    
except PermissionError as e:
    print(f"\n❌ ERRO DE PERMISSÃO:")
    print(f"   Não foi possível escrever em: {arquivo_saida}")
    print(f"   Detalhes: {e}")
    print("\n💡 Soluções:")
    print("   1. Feche o arquivo se estiver aberto no QGIS/ArcGIS")
    print("   2. Execute o script como administrador")
    print("   3. Escolha outro diretório de saída")
    
    # Tentar salvar em local alternativo
    pasta_alternativa = Path.home() / 'Desktop' / 'resultado_intersecao'
    pasta_alternativa.mkdir(parents=True, exist_ok=True)
    arquivo_alternativo = pasta_alternativa / NOME_ARQUIVO
    
    print(f"\n🔄 Tentando salvar em: {arquivo_alternativo}")
    clipped.to_file(arquivo_alternativo, layer=NOME_CAMADA, driver='GPKG')
    print(f"  ✓ Salvo com sucesso no local alternativo!")
    
except Exception as e:
    print(f"\n❌ ERRO INESPERADO:")
    print(f"   Tipo: {type(e).__name__}")
    print(f"   Mensagem: {e}")
    print("\n💡 Verifique:")
    print("   - Espaço em disco disponível")
    print("   - Nome do arquivo (caracteres especiais)")
    print("   - Permissões da pasta")

# ============================================================================
# 7. RESUMO FINAL
# ============================================================================
print("\n" + "="*70)
print("RESUMO DO PROCESSAMENTO")
print("="*70)
print(f"Features de entrada:     {len(input_layer)}")
print(f"Grid original:           {len(overlay_grid)}")
print(f"Grid após filtro:        {len(overlay_grid_filtered)}")
print(f"Features resultantes:    {len(clipped)}")
print(f"Colunas no resultado:    {len(clipped.columns)}")
print(f"Colunas: {', '.join(clipped.columns[:-1])}")  # Exceto geometry
print(f"Tipo de geometria:       MultiPolygon")
print("="*70)

qgis_grade_clip.py

import sys
from pathlib import Path
from qgis.core import (
    QgsVectorLayer,
    QgsVectorFileWriter,
    QgsCoordinateTransformContext,
    QgsProject
)
import processing

# ============================================================
# CONFIGURAÇÕES DE CAMINHOS E CAMADAS
# ============================================================
# Defina o caminho principal do arquivo de entrada e o nome da camada (se for GPKG)
# Exemplo GPKG:
CAMINHO_INPUT = Path("/Users/marcellodebarrosfilho/code/curso_pyqgis_udemy/dados/pb.gpkg")
CAMADA_INPUT_NOME = "estados_ne" # Preencha este campo se o CAMINHO_INPUT for um GPKG

# Exemplo Shapefile (Deixe CAMADA_INPUT_NOME = None ou "")
# CAMINHO_INPUT = Path("/caminho/para/meu_shapefile.shp")
# CAMADA_INPUT_NOME = None

# Definições do GRID (Pluviosidade)
CAMINHO_GRID = Path("/Users/marcellodebarrosfilho/Downloads/pluviosidade_media/Pluviosidade_Media_Anual_1991_2020.gpkg")
CAMADA_GRID = "maps_grade10000"

# --- Lógica de Flexibilidade de Entrada ---

# Define a string de conexão (URI) e o nome base para a saída
if CAMINHO_INPUT.suffix.lower() in ['.gpkg', '.dxf', '.sqlite'] and CAMADA_INPUT_NOME:
    # Se for um container (GPKG), usa a sintaxe layername
    CAMADA_INPUT_URI = f"{CAMINHO_INPUT.as_posix()}|layername={CAMADA_INPUT_NOME}"
    NOME_BASE = CAMADA_INPUT_NOME
else:
    # Para arquivos simples (Shapefile, KML, etc), o URI é o próprio caminho
    CAMADA_INPUT_URI = CAMINHO_INPUT.as_posix()
    NOME_BASE = CAMINHO_INPUT.stem

# Definição automática do caminho e nome da camada de saída
CAMINHO_SAIDA = CAMINHO_INPUT.parent / f"{NOME_BASE}_grade.gpkg"
NOME_CAMADA_SAIDA = f"{NOME_BASE}_grade"

print(f"Arquivo de entrada: {CAMINHO_INPUT} (Camada Base: {NOME_BASE})")
print(f"Arquivo de saída:   {CAMINHO_SAIDA}")

# ============================================================
# 1. CARREGAR CAMADAS
# ============================================================
print("\nCarregando camadas...")
input_layer = QgsVectorLayer(CAMADA_INPUT_URI, "input", "ogr")
grid_layer = QgsVectorLayer(f"{CAMINHO_GRID.as_posix()}|layername={CAMADA_GRID}", "grid", "ogr")

if not input_layer.isValid() or not grid_layer.isValid():
    sys.exit("ERRO CRÍTICO: Não foi possível carregar uma ou mais camadas.")
print("✓ Camadas carregadas")

# ============================================================
# 2. PRÉ-FILTRO DO GRID (Apenas BBox para otimização)
# ============================================================
print("Filtrando grid pela extensão (Bounding Box)...")
bbox = input_layer.extent()

grid_bbox = processing.run("native:extractbyextent", {
    "INPUT": grid_layer,
    "EXTENT": bbox,
    "CLIP": False,
    "OUTPUT": "TEMPORARY_OUTPUT"
})["OUTPUT"]
print("✓ Grid pré-filtrado (Extensão)")

# ============================================================
# 3. INTERSEÇÃO
# ============================================================
print("Executando interseção (Geometria precisa)...")
intersection = processing.run("native:intersection", {
    "INPUT": input_layer,
    "OVERLAY": grid_bbox,
    "INPUT_FIELDS": [],
    "OVERLAY_FIELDS": ["id"],
    "OVERLAY_FIELDS_PREFIX": "",
    "OUTPUT": "TEMPORARY_OUTPUT"
})["OUTPUT"]
print("✓ Interseção concluída")

# ============================================================
# 4. RENOMEAR CAMPO id → grade_id E REMOVER FID DUPLICADO
# ============================================================
print("Reestruturando campos (Renomeando 'id' e removendo 'fid')...")

fields_mapping = []

for field in intersection.fields():
    name = field.name()
    
    # 1. IGNORAR O FID (Obrigatório para evitar o erro UNIQUE constraint)
    if name.lower() == 'fid':
        continue
        
    # 2. Renomear 'id' original do grid para 'grade_id'
    elif name == 'id':
        fields_mapping.append({
            "expression": "id", 
            "name": "grade_id", 
            "type": 4, "length": 0, "precision": 0
        })
        
    # 3. Manter todos os outros campos (municipios, pluviometria, etc)
    else:
        fields_mapping.append({
            "expression": f'"{name}"', 
            "name": name, 
            "type": field.type(), 
            "length": field.length(), 
            "precision": field.precision()
        })

# Executar a refatoração
renamed = processing.run("native:refactorfields", {
    "INPUT": intersection,
    "FIELDS_MAPPING": fields_mapping,
    "OUTPUT": "TEMPORARY_OUTPUT"
})["OUTPUT"]
print("✓ Campos reestruturados (FID removido com sucesso)")

# ============================================================
# 5. CORRIGIR GEOMETRIAS
# ============================================================
print("Corrigindo geometrias...")
fixed = processing.run("native:fixgeometries", {
    "INPUT": renamed,
    "OUTPUT": "TEMPORARY_OUTPUT"
})["OUTPUT"]

# Validação final: Checa se alguma feição não pôde ser consertada
valid_result = processing.run("qgis:checkvalidity", {
    "INPUT_LAYER": fixed,
    "METHOD": 2,  # GEOS
    "IGNORE_RING_SELF_INTERSECTION": False,
    "VALID_OUTPUT": "TEMPORARY_OUTPUT",
    "INVALID_OUTPUT": "TEMPORARY_OUTPUT"
})

final_layer = valid_result["VALID_OUTPUT"]
invalid_count = valid_result["INVALID_OUTPUT"].featureCount()

if invalid_count > 0:
    print(f"⚠️ AVISO: {invalid_count} feições foram descartadas por geometria inválida!")
else:
    print("✓ Todas as geometrias estão válidas.")

# ============================================================
# 6. SALVAR
# ============================================================
print("Salvando arquivo final...")

# Deleta o arquivo existente, tratando erros de permissão
if CAMINHO_SAIDA.exists():
    try:
        CAMINHO_SAIDA.unlink()
        print("  Arquivo antigo removido.")
    except PermissionError:
        sys.exit("ERRO: O arquivo de saída está aberto em outro programa. Feche e tente novamente.")

# Configurações de escrita do GPKG
options = QgsVectorFileWriter.SaveVectorOptions()
options.driverName = "GPKG"
options.layerName = NOME_CAMADA_SAIDA
options.fileEncoding = "UTF-8"

# Usar contexto de transformação do projeto atual
transform_context = QgsProject.instance().transformContext()

error = QgsVectorFileWriter.writeAsVectorFormatV3(
    final_layer,
    CAMINHO_SAIDA.as_posix(),
    transform_context,
    options
)

if error[0] == QgsVectorFileWriter.NoError:
    print("=" * 60)
    print("🎉 PROCESSO CONCLUÍDO COM SUCESSO!")
    print(f"Arquivo criado: {CAMINHO_SAIDA}")
    print(f"Camada: {NOME_CAMADA_SAIDA}")
    print(f"Total de feições: {final_layer.featureCount()}")
    print("=" * 60)
else:
    print(f"✗ Erro ao salvar: {error[1]}")