Entrada 9: Proyecto BermejaPi — Mapeo del Campo Magnético y Detección de la Anomalía del Atlántico Sur (SAA)

En esta novena entrega del proyecto BermejaPi, ampliamos los objetivos científicos de nuestra misión orbital para la Agencia Espacial Europea (ESA). Hasta ahora hemos analizado la Tierra mediante visión por computador (velocidad orbital, albedo y cobertura vegetal); en esta etapa integramos el estudio de la magnetosfera terrestre utilizando el magnetómetro triaxial del Sense HAT.

1. El Objetivo Científico: Geofísica Orbital y la SAA

El campo magnético terrestre (la magnetosfera) actúa como un escudo protector contra la radiación solar y los rayos cósmicos. Al orbitar a 400 km de altitud a bordo de la Estación Espacial Internacional (ISS), el experimento registra fluctuaciones continuas en la intensidad y dirección de dicho campo.

La Anomalía del Atlántico Sur (SAA)

La SAA (South Atlantic Anomaly) es una región geográfica donde el cinturón interior de radiación de Van Allen se aproxima más a la superficie terrestre. En esta zona, el campo magnético global se debilita sensiblemente, aumentando el flujo de partículas cargadas. Medir estas variaciones permite:

• Detectar desviaciones locales en la magnetosfera.
• Correlacionar descensos de intensidad magnética con áreas de mayor exposición radiológica para la instrumentación espacial.

2. Fundamento Matemático y Vectorial

El magnetómetro del Sense HAT mide la densidad de flujo magnético en tres ejes espaciales perpendiculares (Bx, By, Bz) expresada en microteslas (µT).

Para obtener el módulo de la intensidad total del campo magnético (B_total), aplicamos la norma euclídea del vector tridimensional:

B_total = √(Bx² + By² + Bz²)

• Ejes individualizados (Bx, By, Bz): Indican la orientación del vector magnético respecto a la orientación espacial de la cámara e ISS.
• Módulo total (B_total): Representa la fuerza neta del campo magnético independiente de la rotación del módulo. El rango típico en órbita baja terrestre (LEO) oscila entre 25 µT y 65 µT.

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

Esta versión final consolida todas las ramas científicas del proyecto: telemetría ambiental, magnetometría 3D, velocidad orbital (ORB), albedo (HSV), biomasa vegetal (VARI) y purga estricta de memoria según los requisitos de la ESA.

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 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", 
                "mag_x_ut",
                "mag_y_ut",
                "mag_z_ut",
                "mag_total_ut",
                "delta_t_s", 
                "desplazamiento_px", 
                "velocidad_kms",
                "albedo_pct",
                "vari_index",
                "veg_pct"
            ])


def leer_telemetria_sensores():
    """Lee temperatura, presión, humedad y campo magnético 3D del Sense HAT o simula datos."""
    if sense is not None:
        try:
            t = sense.get_temperature()
            p = sense.get_pressure()
            h = sense.get_humidity()
            
            # Obtener lectura del magnetómetro triaxial (microteslas)
            raw_mag = sense.get_compass_raw()
            mx = round(raw_mag['x'], 2)
            my = round(raw_mag['y'], 2)
            mz = round(raw_mag['z'], 2)
            m_total = round(float(np.sqrt(mx**2 + my**2 + mz**2)), 2)
            
            return round(t, 2), round(p, 2), round(h, 2), mx, my, mz, m_total
        except Exception:
            pass
            
    # Valores de simulación para entorno sin Sense HAT conectado
    mx_sim, my_sim, mz_sim = 18.45, -22.10, 35.80
    mtot_sim = round(float(np.sqrt(mx_sim**2 + my_sim**2 + mz_sim**2)), 2)
    return 21.5, 1013.25, 42.0, mx_sim, my_sim, mz_sim, mtot_sim


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
        # Terreno sintético base (BGR)
        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 para 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 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

    hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
    umbral_inferior = np.array([0, 0, 180])
    umbral_superior = np.array([180, 50, 255])

    mascara_nubes = cv2.inRange(hsv, umbral_inferior, umbral_superior)
    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

    b = img[:, :, 0].astype(np.float32)
    g = img[:, :, 1].astype(np.float32)
    r = img[:, :, 2].astype(np.float32)

    denominador = g + r - b
    denominador[denominador == 0] = 1e-6

    matriz_vari = (g - r) / denominador
    mascara_vegetacion = matriz_vari > 0.10

    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")))
    
    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}")

    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, mx, my, mz, mtot, 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,
            mx,
            my,
            mz,
            mtot,
            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 (VISIÓN MULTIESPECTRAL + MAGNETOMETRÍA) ===")
    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 ambiental, geofísica y visión...")
            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

            # Telemetría ambiental y magnetómetro 3D
            temp, presion, humedad, mx, my, mz, mtot = leer_telemetria_sensores()

            # Análisis de visión por computador
            px, v_kms = calcular_velocidad(foto_anterior, foto_actual, delta_t)
            albedo = calcular_albedo_y_nubes(foto_actual)
            vari, veg_pct = calcular_indice_vegetacion(foto_actual)

            print(f"Sensores -> T: {temp}°C | P: {presion} hPa | H: {humedad}%")
            print(f"Magnético-> Bx: {mx} | By: {my} | Bz: {mz} | B_total: {mtot} uT")
            print(f"Visión   -> Shift: {px:.1f}px | V: {v_kms:.2f} km/s | Albedo: {albedo}% | VARI: {vari} | Veg: {veg_pct}%")

            # Registro CSV completo
            registrar_telemetria(temp, presion, humedad, mx, my, mz, mtot, delta_t, px, v_kms, albedo, vari, veg_pct)

            # Purga de almacenamiento
            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)

La base de datos del proyecto integra un total de 14 variables correlacionadas temporalmente:

Columna	Descripción	Unidad
timestamp	Fecha y hora UTC	YYYY-MM-DD HH:MM:SS
temp_c	Temperatura ambiente interna	°C
presion_hpa	Presión atmosférica interna	hPa
humedad_pct	Humedad relativa	%
mag_x_ut	Vector de fuerza magnética en eje X	µT
mag_y_ut	Vector de fuerza magnética en eje Y	µT
mag_z_ut	Vector de fuerza magnética en eje Z	µT
mag_total_ut	Módulo total del campo magnético	µT
delta_t_s	Intervalo temporal entre tomas	Segundos
desplazamiento_px	Desplazamiento medio detectado por ORB	Píxeles
velocidad_kms	Estimación de velocidad orbital	km/s
albedo_pct	Porcentaje de nubes detectadas (HSV)	%
vari_index	Índice medio de vegetación visible	[-1.0, 1.0]
veg_pct	Porcentaje de cubierta vegetal activa	%