Entrada 8: Proyecto BermejaPi — Índice de Vegetación Visible (VARI / ExG) y Visión Multiespectral

En esta octava entrega del proyecto BermejaPi, ampliamos la carga útil científica de nuestro script de vuelo introduciendo el análisis de la biomasa terrestre desde la Estación Espacial Internacional (ISS).

Mediante el uso de índices de vegetación en el espectro visible (VARI y ExG), evaluamos la densidad vegetal del planeta utilizando únicamente la cámara RGB estándar de la Astro Pi.

1. ¿Por qué medir la vegetación sin cámara infrarroja (NoIR)?

Tradicionalmente, los satélites de observación terrestre emplean el índice NDVI (Normalized Difference Vegetation Index), el cual requiere medir la reflectancia en el infrarrojo cercano (NIR). Sin embargo, la cámara visible estándar de la ISS opera exclusivamente en la banda del rojo, verde y azul (RGB).

Para superar esta limitación técnica, aplicamos métricas fotométricas avanzadas en el espectro visible:

• Índice VARI (Visible Atmospherically Resistant Index): Estima la fracción de cubierta vegetal verde minimizando el impacto del volumen atmosférico.
• Índice ExG (Excess Green Index): Resalta la presencia del pigmento de la clorofila contrastando el canal verde frente al rojo y azul.

2. Fundamento Matemático

Dado un fotograma en formato digital, extraemos las matrices de píxeles correspondientes a los canales Rojo ($R$), Verde ($G$) y Azul ($B$) convertidas a precisión de punto flotante (float32) para evitar desbordamientos de memoria (overflow).

Formulación de los Índices

• Fórmula VARI:
  $$\text{VARI} = \frac{G – R}{G + R – B}$$
• Fórmula ExG:
  $$\text{ExG} = 2G – R – B$$

Criterio de Interpretación

• $\text{VARI} > 0.10$ o $\text{ExG} > 0$: Presencia activa de vegetación y clorofila.
• $\text{VARI} \le 0$ o $\text{ExG} \le 0$: Masas de agua, suelo desnudo, desiertos o hielo.

3. Código Fuente Definitivo (main.py)

Esta versión incluye la integración completa del algoritmo VARI, el cálculo de albedo (HSV), la velocidad orbital (ORB), la lectura de sensores del Sense HAT y la regla estricta de purga de almacenamiento asignada por la ESA (máximo 42 imágenes o 240 MB):

Python

import csv
import time
from datetime import datetime, timedelta
from pathlib import Path
import cv2
import numpy as np

# Cargar módulo de cámara oficial
try:
    from picamera2 import Picamera2
except ImportError:
    Picamera2 = None

# Cargar Sense HAT (hardware real o emulador)
try:
    from sense_hat import SenseHat
    sense = SenseHat()
except Exception:
    try:
        from sense_emu import SenseHat
        sense = SenseHat()
    except Exception:
        sense = None

# ============================================
# CONSTANTES DE MISIÓN Y REGLAS ESA
# ============================================
DIR_BASE = Path(__file__).parent.resolve()
FICHERO_CSV = DIR_BASE / "data.csv"

DURACION_MINUTOS = 10
TIEMPO_LIMITE_SEG = DURACION_MINUTOS * 60  # 600 segundos
INTERVALO_CAPTURA_SEG = 15                  # 600s / 15s = 40 capturas max.

MAX_IMAGENES_PERMITIDAS = 42
LIMITE_ALMACENAMIENTO_MB = 240              # Margen de seguridad sobre los 250 MB

ALTURA_ISS_M = 400000.0                     # Altura orbital promedio (400 km)
ANCHO_SENSOR_MM = 7.564                     # HQ Camera
FOCAL_MM = 6.0
ANCHO_IMAGEN_PX = 4056


def inicializar_csv():
    """Crea el archivo data.csv con las cabeceras requeridas si no existe."""
    if not FICHERO_CSV.exists():
        with open(FICHERO_CSV, mode="w", newline="", encoding="utf-8") as f:
            escritor = csv.writer(f)
            escritor.writerow([
                "timestamp", 
                "temp_c", 
                "presion_hpa", 
                "humedad_pct", 
                "delta_t_s", 
                "desplazamiento_px", 
                "velocidad_kms",
                "albedo_pct",
                "vari_index",
                "veg_pct"
            ])


def leer_telemetria_sensores():
    """Lee temperatura, presión y humedad del Sense HAT o genera valores simulados."""
    if sense is not None:
        try:
            t = sense.get_temperature()
            p = sense.get_pressure()
            h = sense.get_humidity()
            return round(t, 2), round(p, 2), round(h, 2)
        except Exception:
            pass
    return 21.5, 1013.25, 42.0


def capturar_imagen(ruta, camara=None):
    """Captura foto con la cámara o genera terreno sintético con nubes y vegetación para VM."""
    if camara is not None:
        camara.capture_file(str(ruta))
    else:
        h, w = 1080, 1920
        # Fondo de terreno sintético en BGR (3 canales)
        terreno_base = np.random.randint(60, 140, (h, w, 3), dtype=np.uint8)
        terreno = cv2.GaussianBlur(terreno_base, (31, 31), 0)
        
        # Simulación de nubes brillantes (blancas) para probar Albedo
        cv2.circle(terreno, (600, 400), 180, (245, 245, 245), -1)
        cv2.circle(terreno, (1200, 700), 250, (250, 250, 250), -1)
        
        # Simulación de masa vegetal (verde) para probar el índice VARI (BGR)
        cv2.rectangle(terreno, (200, 200), (700, 800), (35, 160, 45), -1)
        
        cv2.imwrite(str(ruta), terreno)


def calcular_albedo_y_nubes(ruta_imagen):
    """Calcula el porcentaje de cobertura nubosa (albedo estimado) usando el espacio HSV."""
    img = cv2.imread(str(ruta_imagen))
    if img is None:
        return 0.0

    # 1. Convertir a espacio de color HSV
    hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)

    # 2. Definir rango para nubes (baja saturación, alto brillo/valor)
    umbral_inferior = np.array([0, 0, 180])
    umbral_superior = np.array([180, 50, 255])

    # 3. Crear máscara binaria (Blanco = Nube, Negro = Terreno/Océano)
    mascara_nubes = cv2.inRange(hsv, umbral_inferior, umbral_superior)

    # 4. Calcular proporción de píxeles
    pixeles_nubes = cv2.countNonZero(mascara_nubes)
    pixeles_totales = img.shape[0] * img.shape[1]

    porcentaje_albedo = (pixeles_nubes / pixeles_totales) * 100.0
    return round(porcentaje_albedo, 2)


def calcular_indice_vegetacion(ruta_imagen):
    """Calcula el índice VARI y el porcentaje de vegetación visible usando OpenCV/NumPy."""
    img = cv2.imread(str(ruta_imagen))
    if img is None:
        return 0.0, 0.0

    # Convertir canales a float32 para evitar desbordamiento de enteros uint8 (BGR)
    b = img[:, :, 0].astype(np.float32)
    g = img[:, :, 1].astype(np.float32)
    r = img[:, :, 2].astype(np.float32)

    # Evitar división por cero añadiendo constante épsilon (1e-6)
    denominador = g + r - b
    denominador[denominador == 0] = 1e-6

    # 1. Cálculo de la matriz VARI = (G - R) / (G + R - B)
    matriz_vari = (g - r) / denominador

    # 2. Máscara de vegetación activa (píxeles con VARI positivo significativo > 0.10)
    mascara_vegetacion = matriz_vari > 0.10

    # 3. Métricas
    vari_medio = float(np.mean(matriz_vari))
    pixeles_vegetacion = int(np.sum(mascara_vegetacion))
    pixeles_totales = img.shape[0] * img.shape[1]

    porcentaje_vegetacion = (pixeles_vegetacion / pixeles_totales) * 100.0

    return round(vari_medio, 3), round(porcentaje_vegetacion, 2)


def calcular_velocidad(ruta1, ruta2, delta_t):
    """Calcula desplazamiento y velocidad orbital mediante ORB de OpenCV."""
    img1 = cv2.imread(str(ruta1), cv2.IMREAD_GRAYSCALE)
    img2 = cv2.imread(str(ruta2), cv2.IMREAD_GRAYSCALE)

    if img1 is None or img2 is None:
        return 0.0, 0.0

    orb = cv2.ORB_create(nfeatures=1000)
    kp1, des1 = orb.detectAndCompute(img1, None)
    kp2, des2 = orb.detectAndCompute(img2, None)

    if des1 is None or des2 is None:
        return 0.0, 0.0

    bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=True)
    coincidencias = bf.match(des1, des2)
    
    if not coincidencias:
        return 0.0, 0.0

    coincidencias = sorted(coincidencias, key=lambda x: x.distance)[:50]

    distancias_px = [
        np.sqrt((kp2[m.trainIdx].pt[0] - kp1[m.queryIdx].pt[0])**2 + 
                (kp2[m.trainIdx].pt[1] - kp1[m.queryIdx].pt[1])**2)
        for m in coincidencias
    ]

    desplazamiento_px = float(np.mean(distancias_px))
    gsd = (ALTURA_ISS_M * ANCHO_SENSOR_MM) / (FOCAL_MM * ANCHO_IMAGEN_PX)
    distancia_m = desplazamiento_px * gsd
    velocidad_kms = (distancia_m / delta_t) / 1000.0 if delta_t > 0 else 0.0

    return desplazamiento_px, velocidad_kms


def verificar_y_purgar_almacenamiento():
    """Garantiza no superar 42 fotos ni el límite de espacio asignado por la ESA."""
    imagenes_guardadas = sorted(list(DIR_BASE.glob("foto_*.jpg")))
    
    # 1. Purga por número máximo de imágenes (Regla ESA <= 42)
    while len(imagenes_guardadas) > MAX_IMAGENES_PERMITIDAS:
        foto_antigua = imagenes_guardadas.pop(0)
        try:
            foto_antigua.unlink()
            print(f"[Purga ESA] Eliminada imagen antigua: {foto_antigua.name}")
        except Exception as e:
            print(f"Error al eliminar {foto_antigua.name}: {e}")

    # 2. Control de espacio total en megabytes
    tamano_total_bytes = sum(f.stat().st_size for f in DIR_BASE.glob("*") if f.is_file())
    tamano_total_mb = tamano_total_bytes / (1024 * 1024)

    if tamano_total_mb > LIMITE_ALMACENAMIENTO_MB and imagenes_guardadas:
        foto_a_borrar = imagenes_guardadas.pop(0)
        try:
            foto_a_borrar.unlink()
            print(f"[Alerta Almacenamiento] Eliminada por espacio: {foto_a_borrar.name}")
        except Exception as e:
            print(f"Error al liberar espacio: {e}")


def registrar_telemetria(temp, presion, humedad, delta_t, px, v_kms, albedo_pct, vari_index, veg_pct):
    """Guarda la fila de datos completa en data.csv."""
    ahora = datetime.now().strftime("%Y-%m-%d %H:%M:%S")
    with open(FICHERO_CSV, mode="a", newline="", encoding="utf-8") as f:
        escritor = csv.writer(f)
        escritor.writerow([
            ahora, 
            temp, 
            presion, 
            humedad, 
            round(delta_t, 2), 
            round(px, 2), 
            round(v_kms, 2),
            albedo_pct,
            vari_index,
            veg_pct
        ])


def bucle_principal_mision():
    print("=== INICIANDO MISIÓN BERMEJAPI (ALBEDO, VELOCIDAD Y VEGETACIÓN VARI) ===")
    inicializar_csv()

    camara = None
    if Picamera2 is not None:
        try:
            camara = Picamera2()
            camara.start()
        except Exception as e:
            print(f"Modo simulación cámara activo: {e}")

    hora_inicio = datetime.now()
    hora_fin = hora_inicio + timedelta(seconds=TIEMPO_LIMITE_SEG - 10)

    contador_fotos = 1
    foto_anterior = DIR_BASE / f"foto_{contador_fotos:02d}.jpg"
    
    capturar_imagen(foto_anterior, camara)
    tiempo_anterior = time.time()

    ciclo = 1
    while datetime.now() < hora_fin:
        try:
            print(f"\n[Ciclo {ciclo}] Procesando telemetría, visión multiespectral y velocidad...")
            time.sleep(INTERVALO_CAPTURA_SEG)

            contador_fotos += 1
            foto_actual = DIR_BASE / f"foto_{contador_fotos:02d}.jpg"

            capturar_imagen(foto_actual, camara)
            tiempo_actual = time.time()
            delta_t = tiempo_actual - tiempo_anterior

            # Lectura de sensores del Sense HAT
            temp, presion, humedad = leer_telemetria_sensores()

            # Cálculo de velocidad orbital (ORB)
            px, v_kms = calcular_velocidad(foto_anterior, foto_actual, delta_t)

            # Cálculo de albedo/cobertura nubosa (HSV)
            albedo = calcular_albedo_y_nubes(foto_actual)

            # Cálculo del índice de vegetación visible (VARI)
            vari, veg_pct = calcular_indice_vegetacion(foto_actual)

            print(f"Sensores -> T: {temp}°C | P: {presion} hPa | H: {humedad}%")
            print(f"Visión   -> Shift: {px:.1f}px | V: {v_kms:.2f} km/s | Albedo: {albedo}% | VARI: {vari} | Veg: {veg_pct}%")

            # Registro de datos en CSV
            registrar_telemetria(temp, presion, humedad, delta_t, px, v_kms, albedo, vari, veg_pct)

            # Control de purga y almacenamiento según reglas ESA
            verificar_y_purgar_almacenamiento()

            foto_anterior = foto_actual
            tiempo_anterior = tiempo_actual
            ciclo += 1

        except Exception as e:
            print(f"Excepción controlada en ciclo {ciclo}: {e}")
            time.sleep(2)

    if camara is not None:
        camara.stop()

    verificar_y_purgar_almacenamiento()
    print("\n=== PRUEBA DE MISIÓN COMPLETADA CORRECTAMENTE ===")


if __name__ == "__main__":
    bucle_principal_mision()

4. Estructura del Fichero de Salida (data.csv)

Con esta actualización, cada ciclo de ejecución añade 10 variables clave al fichero de registro:

Columna	Descripción	Rango / Unidad
timestamp	Marca temporal UTC	YYYY-MM-DD HH:MM:SS
temp_c	Temperatura interna	$\text{°C}$
presion_hpa	Presión atmosférica	$\text{hPa}$
humedad_pct	Humedad relativa	$\%$
delta_t_s	Intervalo entre capturas	Segundos
desplazamiento_px	Desplazamiento medio de características	Píxeles
velocidad_kms	Velocidad orbital estimada	$\text{km/s}$
albedo_pct	Cobertura nubosa estimada (HSV)	$\%$
vari_index	Índice medio de vegetación visible	$[-1.0, 1.0]$
veg_pct	Porcentaje de cubierta vegetal fértil	$\%$