run_goes16_analysis.py

#!/usr/bin/env python3
"""
Pipeline de Detecção de Queimadas GOES-16 (K-Means) para dados em cache.
Executa análise completa: carregamento → K-Means → métricas → plotagem.
"""
import os, sys, json, math
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
import pyproj
from datetime import datetime

# Config paths
BASE_DIR = "/Users/naubergois/sistemasatelitefunceme"
OUTPUT_DIR = "/Users/naubergois/gois/climate-reports"
os.makedirs(OUTPUT_DIR, exist_ok=True)

NC_B07 = os.path.join(BASE_DIR, "b07.nc")
NC_B13 = os.path.join(BASE_DIR, "b13.nc")

print("=" * 70)
print(f"🛰️  GOES-16 QUEIMADAS — Análise K-Means (cache local)")
print(f"📅  Execução: {datetime.now().strftime('%Y-%m-%d %H:%M:%S')} UTC")
print("=" * 70)

# 1. Verify cache files
for f in [NC_B07, NC_B13]:
    if not os.path.exists(f):
        print(f"❌ Arquivo não encontrado: {f}")
        sys.exit(1)
    size_mb = os.path.getsize(f) / 1e6
    print(f"✅ Cache: {os.path.basename(f)} ({size_mb:.1f} MB)")

# 2. Load NetCDF data
nc07 = Dataset(NC_B07, 'r')
nc13 = Dataset(NC_B13, 'r')

t07 = nc07.variables['CMI'][:]  # Band 07 (3.9µm, Shortwave IR)
t13 = nc13.variables['CMI'][:]  # Band 13 (10.3µm, Clean IR)

print(f"📊  Band 07 shape: {t07.shape}, dtype: {t07.dtype}")
print(f"📊  Band 13 shape: {t13.shape}, dtype: {t13.dtype}")

# 3. Extract timestamp from NC file attributes
try:
    time_coverage = getattr(nc07, 'time_coverage_start', None)
    if time_coverage:
        print(f"🕐  Cobertura temporal: {time_coverage}")
except:
    pass

# 4. Georeferencing
proj_var = nc07.variables['goes_imager_projection']
h = proj_var.perspective_point_height
lon_0 = proj_var.longitude_of_projection_origin
p = pyproj.Proj(proj='geos', h=h, lon_0=lon_0, sweep='x')
x = nc07.variables['x'][:] * h
y = nc07.variables['y'][:] * h
xx, yy = np.meshgrid(x, y)
lons, lats = p(xx, yy, inverse=True)

print(f"🌍  Projeção GOES: h={h}, lon_0={lon_0}")
print(f"🌍  Grid: {len(x)}x{len(y)} pixels")

# 5. Ceará mask
# Rough bounding box for Ceará state
LAT_MIN, LAT_MAX = -7.9, -2.7
LON_MIN, LON_MAX = -41.5, -37.0

mask_ceara = (lats >= LAT_MIN) & (lats <= LAT_MAX) & (lons >= LON_MIN) & (lons <= LON_MAX)
n_ceara = np.sum(mask_ceara)
print(f"📍  Pixels no Ceará: {n_ceara}")

if n_ceara < 100:
    print("❌ Muito poucos pixels no Ceará. Abortando.")
    sys.exit(1)

# 6. Feature extraction
features_b07 = t07[mask_ceara].flatten()
features_b13 = t13[mask_ceara].flatten()

valid_pixels = ~np.isnan(features_b07) & ~np.isnan(features_b13)
n_valid = np.sum(valid_pixels)

print(f"📊  Pixels válidos (não-NaN): {n_valid}")

if n_valid < 100:
    print("❌ Poucos pixels válidos. Abortando.")
    sys.exit(1)

# Feature matrix: B07, B13, B07-B13 difference
X_valid = np.column_stack((
    features_b07[valid_pixels],
    features_b13[valid_pixels],
    (features_b07 - features_b13)[valid_pixels]
))

# Statistics
b07_mean = np.mean(X_valid[:, 0])
b07_std = np.std(X_valid[:, 0])
b13_mean = np.mean(X_valid[:, 1])
b13_std = np.std(X_valid[:, 1])
diff_mean = np.mean(X_valid[:, 2])
diff_std = np.std(X_valid[:, 2])

print(f"\n📊  Estatísticas dos pixels do Ceará:")
print(f"     B07 (3.9µm):  μ={b07_mean:.2f}K, σ={b07_std:.2f}K")
print(f"     B13 (10.3µm): μ={b13_mean:.2f}K, σ={b13_std:.2f}K")
print(f"     B07-B13:      μ={diff_mean:.2f}K, σ={diff_std:.2f}K")

# 7. K-Means clustering
n_clusters = 4
print(f"\n🔬  Executando K-Means (k={n_clusters})...")

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_valid)

kmeans = KMeans(n_clusters=n_clusters, random_state=42, n_init=10)
labels = kmeans.fit_predict(X_scaled)

# 8. Identify fire cluster
cluster_stats = []
for i in range(n_clusters):
    mask = labels == i
    count = np.sum(mask)
    mean_b07 = np.mean(X_valid[mask, 0])
    mean_b13 = np.mean(X_valid[mask, 1])
    mean_diff = np.mean(X_valid[mask, 2])
    std_b07 = np.std(X_valid[mask, 0])
    cluster_stats.append({
        'cluster_id': i,
        'count': int(count),
        'pct': float(count / n_valid * 100),
        'mean_b07': float(mean_b07),
        'mean_b13': float(mean_b13),
        'mean_diff': float(mean_diff),
        'std_b07': float(std_b07)
    })
    print(f"     Cluster {i}: {count:>6} pixels ({count/n_valid*100:5.1f}%) | "
          f"B07 μ={mean_b07:.1f}K, σ={std_b07:.1f}K | "
          f"diff μ={mean_diff:.1f}K")

# Sort by B07 temperature (descending) — hottest cluster = fire candidate
cluster_stats.sort(key=lambda c: c['mean_b07'], reverse=True)
fire_cluster = cluster_stats[0]
fire_cluster_id = fire_cluster['cluster_id']
fire_temp_avg = fire_cluster['mean_b07']

print(f"\n🔥  Cluster candidato a fogo: ID={fire_cluster_id}")
print(f"🔥  Temperatura média B07: {fire_temp_avg:.2f}K")
print(f"🔥  Pixels candidatos: {fire_cluster['count']} ({fire_cluster['pct']:.1f}%)")

# 9. Refinement: apply fire threshold filters
# Classic GOES fire detection: B07 > 315K AND B07-B13 > 5K
fire_mask_simple = (X_valid[:, 0] > 315.0)

# K-Means + filter (our approach)
fire_mask_km_filter = (labels == fire_cluster_id) & (X_valid[:, 0] > 315.0)

# K-Means + filter + difference
fire_mask_full = (labels == fire_cluster_id) & (X_valid[:, 0] > 315.0) & (X_valid[:, 2] > 5.0)

n_simple = np.sum(fire_mask_simple)
n_km = np.sum(labels == fire_cluster_id)
n_km_filter = np.sum(fire_mask_km_filter)
n_full = np.sum(fire_mask_full)

print(f"\n📋  Métricas de detecção (Ceará):")
print(f"     Simple threshold (B07 > 315K):           {n_simple:>5} pixels")
print(f"     K-Means cluster only:                    {n_km:>5} pixels")
print(f"     K-Means + B07 > 315K:                    {n_km_filter:>5} pixels")
print(f"     K-Means + B07 > 315K + Diff > 5K:        {n_full:>5} pixels")

# 10. Analyze the hottest fire pixels
if n_full > 0:
    fire_pixels = X_valid[fire_mask_full]
    hottest_idx = np.argmax(fire_pixels[:, 0])
    hottest_b07 = fire_pixels[hottest_idx, 0]
    hottest_b13 = fire_pixels[hottest_idx, 1]
    hottest_diff = fire_pixels[hottest_idx, 2]
    print(f"\n🔥  Pixel mais quente: B07={hottest_b07:.1f}K, B13={hottest_b13:.1f}K, Δ={hottest_diff:.1f}K")
    print(f"🔥  Média dos pixels de fogo: B07 μ={np.mean(fire_pixels[:,0]):.1f}K")

# 11. Reconstruct fire mask on original image
final_fire_mask = np.zeros(t07.shape, dtype=bool)
ceara_labels = np.full(features_b07.shape, -1)
ceara_labels[valid_pixels] = labels

# Full approach
is_fire_flat = (ceara_labels == fire_cluster_id) & (features_b07 > 315.0) & ((features_b07 - features_b13) > 5.0)
is_fire_flat_full = is_fire_flat & valid_pixels  # only where valid

final_fire_mask[mask_ceara] = is_fire_flat_full

# 12. Generate visualizations
print("\n🎨  Gerando visualizações...")

# Find bounding box for Ceará crop
rows, cols = np.where(mask_ceara)
if len(rows) > 0:
    y_min, y_max = rows.min(), rows.max()
    x_min, x_max = cols.min(), cols.max()

    t07_crop = t07[y_min:y_max+1, x_min:x_max+1]
    x_crop = x[x_min:x_max+1]
    y_crop = y[y_min:y_max+1]
    fire_mask_crop = final_fire_mask[y_min:y_max+1, x_min:x_max+1]
    lons_crop = lons[y_min:y_max+1, x_min:x_max+1]
    lats_crop = lats[y_min:y_max+1, x_min:x_max+1]

    # --- Plot 1: Thermal + Fire detections ---
    plt.figure(figsize=(14, 10))
    plt.imshow(t07_crop, cmap='inferno', vmin=280, vmax=340,
              extent=[x_crop.min(), x_crop.max(), y_crop.min(), y_crop.max()])
    cb = plt.colorbar(label='Brightness Temperature (K)')

    if np.sum(fire_mask_crop) > 0:
        yy_fire, xx_fire = np.where(fire_mask_crop)
        xx_grid, yy_grid = np.meshgrid(x_crop, y_crop)
        plt.scatter(xx_grid[fire_mask_crop], yy_grid[fire_mask_crop],
                   c='cyan', marker='x', s=40, linewidths=0.8,
                   label=f'Fire Detected ({np.sum(fire_mask_crop)} pixels)')
        plt.legend(fontsize=12, loc='upper right')

    plt.title(f'GOES-16 K-Means Fire Detection — Ceará\n'
              f'Band 07 + 13 | k=4 | B07>315K | Δ>5K | {np.sum(fire_mask_crop)} fire pixels',
              fontsize=14, fontweight='bold')
    plt.xlabel('x (m)', fontsize=11)
    plt.ylabel('y (m)', fontsize=11)
    plt.tight_layout()
    thermal_path = os.path.join(OUTPUT_DIR, "goes16_thermal_fire.png")
    plt.savefig(thermal_path, dpi=150)
    plt.close()
    print(f"   ✅ Thermal plot: {thermal_path}")

    # --- Plot 2: Scatter plot B07 vs B13 colored by cluster ---
    plt.figure(figsize=(12, 8))
    colors = plt.cm.tab10(np.linspace(0, 1, n_clusters))
    for i in range(n_clusters):
        mask = labels == i
        label_str = f'Cluster {i}' + (' (FIRE)' if i == fire_cluster_id else '')
        alpha = 0.8 if i == fire_cluster_id else 0.3
        s = 15 if i == fire_cluster_id else 5
        plt.scatter(X_valid[mask, 1], X_valid[mask, 0],
                   c=[colors[i]], label=label_str, alpha=alpha, s=s, edgecolors='none')

    plt.axhline(y=315, color='red', linestyle='--', alpha=0.7, label='B07 > 315K (threshold)')
    plt.axhline(y=300, color='orange', linestyle=':', alpha=0.5, label='300K reference')
    plt.xlabel('B13 Brightness Temperature (K)', fontsize=12)
    plt.ylabel('B07 Brightness Temperature (K)', fontsize=12)
    plt.title(f'K-Means Clustering — B07 vs B13\nCeará | {n_valid} pixels | k={n_clusters}',
              fontsize=14, fontweight='bold')
    plt.grid(alpha=0.3)
    plt.legend(fontsize=10)
    plt.tight_layout()
    scatter_path = os.path.join(OUTPUT_DIR, "goes16_cluster_scatter.png")
    plt.savefig(scatter_path, dpi=150)
    plt.close()
    print(f"   ✅ Scatter plot: {scatter_path}")

    # --- Plot 3: B07-B13 difference histogram ---
    plt.figure(figsize=(12, 6))
    fire_mask_cl = labels == fire_cluster_id
    non_fire_mask = labels != fire_cluster_id

    plt.hist(X_valid[non_fire_mask, 2], bins=100, alpha=0.5, label='Non-fire clusters',
             color='gray', density=True)
    if np.sum(fire_mask_cl) > 0:
        plt.hist(X_valid[fire_mask_cl, 2], bins=50, alpha=0.7, label='Fire cluster',
                 color='red', density=True)
    plt.axvline(x=5, color='red', linestyle='--', alpha=0.7, label='Δ > 5K (fire criterion)')
    plt.xlabel('B07 - B13 Difference (K)', fontsize=12)
    plt.ylabel('Density', fontsize=12)
    plt.title('B07-B13 Temperature Difference Distribution by Cluster', fontsize=14, fontweight='bold')
    plt.legend(fontsize=11)
    plt.grid(alpha=0.3)
    plt.tight_layout()
    hist_path = os.path.join(OUTPUT_DIR, "goes16_diff_histogram.png")
    plt.savefig(hist_path, dpi=150)
    plt.close()
    print(f"   ✅ Histogram: {hist_path}")

else:
    print("⚠️  Could not crop Ceará region for plotting.")

# 13. Final metrics summary
print("\n" + "=" * 70)
print("📊  RELATÓRIO FINAL — GOES-16 K-MEANS QUEIMADAS CEARÁ")
print("=" * 70)

print(f"""
📅  Data da análise:       {datetime.now().strftime('%Y-%m-%d %H:%M:%S')} UTC
📁  Cache GOES-16:          b07.nc + b13.nc (Jan 2026)
📍  Região:                 Ceará (BBox: {LON_MIN}°W a {LON_MAX}°W, {LAT_MIN}°S a {LAT_MAX}°S)
🔬  Algoritmo:              K-Means (k={n_clusters}, random_state=42)
📐  Features:               B07, B13, B07-B13

📋  DISTRIBUIÇÃO DOS CLUSTERS:
""")

for c in sorted(cluster_stats, key=lambda x: x['cluster_id']):
    cid = c['cluster_id']
    tag = "🔥 FOGO" if cid == fire_cluster_id else "     "
    print(f"  {tag}  Cluster {cid}: {c['count']:>6} pixels ({c['pct']:5.1f}%)  "
          f"B07={c['mean_b07']:.1f}K  B13={c['mean_b13']:.1f}K  Δ={c['mean_diff']:.1f}K")

print(f"""
📋  DETECÇÃO DE FOGO (Ceará):
     Método                    |  Pixels de Fogo
     ─────────────────────────┼───────────────
     Simple (B07 > 315K)      |  {n_simple:>6}
     K-Means cluster          |  {n_km:>6}
     K-Means + B07 > 315K     |  {n_km_filter:>6}
     K-Means + 315K + Δ>5K    |  {n_full:>6}

📊  VALIDAÇÃO CRUZADA (INPE BDQueimadas):
     Fonte:                    INPE AQUA Tarde (satélite referência)
     Período:                  01/01/2026 a 04/06/2026
     Focos CE (anual):         600 focos
     Focos CE (últ. 48h):      1 foco
     Ranking CE no Brasil:     8º lugar

     ⚠️  Nota: GOES-16 (geoestacionário) e AQUA Tarde (polar)
         têm resoluções espaciais e temporais diferentes.
         GOES-16 detecta anomalias termais a cada 10 min,
         enquanto AQUA Tarde tem ~2 passagens/dia.
         A comparação direta é indicativa, não exata.
""")

# Save metrics JSON
metrics = {
    'timestamp': datetime.now().strftime('%Y-%m-%d %H:%M:%S'),
    'region': 'Ceará',
    'bbox': {'lon_min': LON_MIN, 'lon_max': LON_MAX, 'lat_min': LAT_MIN, 'lat_max': LAT_MAX},
    'algorithm': f'KMeans(k={n_clusters})',
    'features': ['B07', 'B13', 'B07-B13'],
    'n_pixels_valid': int(n_valid),
    'clusters': cluster_stats,
    'fire_cluster_id': int(fire_cluster_id),
    'fire_detections': {
        'simple_threshold_315K': int(n_simple),
        'kmeans_only': int(n_km),
        'kmeans_plus_315K': int(n_km_filter),
        'kmeans_plus_315K_plus_diff5K': int(n_full),
    },
    'inpe_reference': {
        'satellite': 'AQUA Tarde',
        'source': 'https://terrabrasilis.dpi.inpe.br/queimadas/portal/',
        'period_annual': '2026-01-01 to 2026-06-04',
        'focos_ce_annual': 600,
        'focos_ce_48h': 1,
        'ranking_ce_brasil': 8,
        'total_brasil_annual': 14995,
    },
    'hottest_fire_pixel': {
        'b07_k': float(hottest_b07) if n_full > 0 else None,
        'b13_k': float(hottest_b13) if n_full > 0 else None,
        'diff_k': float(hottest_diff) if n_full > 0 else None,
    } if n_full > 0 else None,
    'output_files': {
        'thermal_map': 'goes16_thermal_fire.png',
        'scatter_plot': 'goes16_cluster_scatter.png',
        'diff_histogram': 'goes16_diff_histogram.png',
    }
}

metrics_path = os.path.join(OUTPUT_DIR, "goes16_metrics.json")
with open(metrics_path, 'w') as f:
    json.dump(metrics, f, indent=2, ensure_ascii=False)

print(f"✅  Métricas salvas: {metrics_path}")
print("✅  Análise concluída com sucesso.")