317 lines
14 KiB
Python
317 lines
14 KiB
Python
#!/usr/bin/env python
|
|
# -*- coding: utf-8 -*-
|
|
"""
|
|
Estimacion de la irradiancia global REAL (correccion de la saturacion
|
|
del piranometro en 1024 W/m2) para una sesion de campo.
|
|
|
|
Ubicacion de las pruebas: Aguascalientes, Aguascalientes, Mexico.
|
|
lat = +21.8756, lon = -102.2854, alt = 1896 m, TZ = UTC-6 (sin DST)
|
|
|
|
Metodo:
|
|
1) Posicion solar por formulas estandar (declinacion de Cooper,
|
|
ecuacion de tiempo de Spencer, correccion de longitud respecto al
|
|
meridiano del huso UTC-6 a 90 W, masa de aire de Kasten-Young).
|
|
2) Modelo de cielo despejado G_cs = A*exp(-B*AM) + C (transmisividad
|
|
de Beer sobre la masa de aire + difuso de cielo claro), ajustado SOLO
|
|
a muestras de cielo claro SIN saturar (el clipping a 1024 W/m2
|
|
sesgaria el ajuste hacia abajo). Para anclar la forma de la curva,
|
|
el ajuste se toma sobre el sobre-ensvolvente de la medida (muestra
|
|
cercana al maximo movil de 30 min), porque las nubes solo bajan la
|
|
curva de cielo despejado, nunca la suben.
|
|
3) La correccion es MONOTONICA: a las muestras saturadas (G >= 1024)
|
|
se les asigna G_est = max(1024, G_cs) (la correccion nunca baja la
|
|
meda); al resto se les conserva la medida (nubes incluidas).
|
|
4) Se recalculan con la misma ventana que el articulo (G > 50 W/m2, y
|
|
G > 600 para las estadisticas de cielo claro): H, R_opt, P_ref,
|
|
eta_tr, mediana P/Pref, fraccion P >= 0.95 P_ref y t95.
|
|
|
|
Uso: py -3.12 clearsky.py <tag> [nofig]
|
|
Salida: articulo/resultados_<tag>/preview_irradiance_est.png
|
|
articulo/resultados_<tag>/irradiance_est_<tag>.json
|
|
"""
|
|
|
|
import json
|
|
import sys
|
|
import tempfile
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
|
|
import matplotlib
|
|
matplotlib.use("Agg")
|
|
import matplotlib.dates as mdates
|
|
import matplotlib.pyplot as plt
|
|
|
|
import parse_log as pl
|
|
|
|
ROOT = Path(__file__).resolve().parent
|
|
|
|
# Aguascalientes, Ags.
|
|
LAT = 21.8756
|
|
LON_W = 102.2854 # longitud oeste positiva
|
|
ALT_M = 1896.0
|
|
L_STD_W = 90.0 # meridiano del huso UTC-6
|
|
G_CLIP = 1023.9 # ADC a tope => G = 1024.0 W/m2
|
|
|
|
|
|
def solar_geometry(ts: pd.Series):
|
|
"""Elevacion solar y masa de aire (Kasten-Young) a partir de la hora local."""
|
|
n = (ts.dt.dayofyear.values * 86400.0
|
|
+ ts.dt.hour.values * 3600.0 + ts.dt.minute.values * 60.0
|
|
+ ts.dt.second.values) / 86400.0
|
|
rad = np.pi / 180.0
|
|
phi = LAT * rad
|
|
B = 2 * np.pi / 365.0 * (n - 81.0)
|
|
delta = 23.45 * rad * np.sin(2 * np.pi / 365.0 * (n + 284.0))
|
|
eot_min = 229.18 * (0.000075 + 0.001868 * np.cos(B)
|
|
- 0.032077 * np.sin(B)
|
|
- 0.014615 * np.cos(2 * B)
|
|
- 0.040849 * np.sin(2 * B))
|
|
clock_min = (ts.dt.hour.values * 3600 + ts.dt.minute.values * 60
|
|
+ ts.dt.second.values) / 60.0
|
|
l_corr = 4.0 * (LON_W - L_STD_W) # minutos (+ si esta al oeste)
|
|
hour_angle = 15.0 * ((clock_min - l_corr + eot_min) / 60.0 - 12.0)
|
|
sin_el = (np.sin(phi) * np.sin(delta)
|
|
+ np.cos(phi) * np.cos(delta) * np.cos(hour_angle * rad))
|
|
el = np.degrees(np.arcsin(np.clip(sin_el, -1, 1)))
|
|
z = 90.0 - el
|
|
m_air = np.where(z < 75.0,
|
|
1.0 / (np.cos(z * rad) + 0.14 * (93.85 - z) ** (-1.375)),
|
|
np.inf)
|
|
return el, m_air
|
|
|
|
|
|
def fit_clearsky(el_deg, am, G, cloud_mask, t_idx):
|
|
"""Ajusta el modelo G_cs = A*exp(-B*AM) + C a las muestras claras no saturadas.
|
|
|
|
Conjunto de ajuste: el sobre-ensvolvente de la medida (muestra >= 0.98
|
|
del maximo movil de 30 min), fuera de transitorios nubosos, sin
|
|
saturar (G < 1015) y con el sol por encima de 12 grados. Al anclarse
|
|
en la parte alta de la curva, la extrapolacion al pico (que cae en la
|
|
zona saturada del piranometro) es corta y confiable.
|
|
"""
|
|
roll = pd.Series(G, index=pd.DatetimeIndex(t_idx)).rolling(
|
|
"30min", center=True, min_periods=5).max().values
|
|
fitm = ((el_deg > 12.0) & (G > 500.0) & (G < 1015.0)
|
|
& (~cloud_mask) & (G >= 0.98 * roll))
|
|
Gs = G[fitm]
|
|
AMs = am[fitm]
|
|
n_fit = int(fitm.sum())
|
|
if n_fit < 20:
|
|
raise RuntimeError("Muy pocas muestras claras para ajustar el modelo "
|
|
f"(n={n_fit}).")
|
|
|
|
def sse_of(A, Bc, C):
|
|
return float(np.sum((Gs - (A * np.exp(-Bc * AMs) + C)) ** 2))
|
|
|
|
best = None
|
|
for A in range(900, 1501, 50):
|
|
for Bc in np.arange(0.05, 0.61, 0.02):
|
|
for C in range(0, 301, 10):
|
|
s = sse_of(A, Bc, C)
|
|
if best is None or s < best[0]:
|
|
best = (s, float(A), float(Bc), float(C))
|
|
# refinamiento local
|
|
_, A0, B0, C0 = best
|
|
for A in range(int(A0) - 50, int(A0) + 51, 10):
|
|
for Bc in np.arange(B0 - 0.02, B0 + 0.021, 0.005):
|
|
for C in range(int(C0) - 10, int(C0) + 11, 5):
|
|
if A <= 0 or Bc <= 0 or C < 0:
|
|
continue
|
|
s = sse_of(A, Bc, C)
|
|
if best is not None and s < best[0]:
|
|
best = (s, float(A), float(Bc), float(C))
|
|
sse, A, Bc, C = best
|
|
rms = float(np.sqrt(sse / n_fit))
|
|
am_finite = np.isfinite(am)
|
|
am_safe = np.where(am_finite, am, 0.0)
|
|
# La curva solo se evalua donde la masa de aire es finita (el>=15). Para
|
|
# 0<el<15 la formula de Kasten-Young diverge (am=inf); reportar
|
|
# A*exp(-B*0)=A ahi seria fisicamente incorrecto.
|
|
curve = np.where((el_deg > 0.0) & am_finite,
|
|
A * np.exp(-Bc * am_safe) + C, 0.0)
|
|
# Pico del modelo a mediodia (masa de aire minima) = irradiancia maxima
|
|
# de cielo despejado en el sitio.
|
|
if am_finite.any():
|
|
am_min = float(np.min(am[am_finite]))
|
|
else:
|
|
am_min = 1.0
|
|
model_peak = A * np.exp(-Bc * am_min) + C
|
|
return A, Bc, C, rms, n_fit, curve, model_peak
|
|
|
|
|
|
def run(tag, nofig=False):
|
|
from compare_sessions import normalize
|
|
src = ROOT / f"{tag}.txt"
|
|
if not src.exists():
|
|
raise SystemExit(f"Log {src.name} no encontrado en la raiz.")
|
|
tmp = Path(tempfile.gettempdir()) / (tag + "_clearsky.txt")
|
|
normalize(src, tmp)
|
|
df_s, df_e = pl.parse_log(tmp)
|
|
df, ph = pl.build_physical(df_s, df_e)
|
|
|
|
G = df["G_POA"].values
|
|
Voc = df["Voc"].values
|
|
Isc = df["Isc"].values
|
|
Pc = df["P_CPV"].values
|
|
t = df["timestamp"]
|
|
|
|
el, am = solar_geometry(df["timestamp"])
|
|
ces = pl.cloud_events(df)
|
|
cloud_mask = np.zeros(len(df), dtype=bool)
|
|
for c in ces:
|
|
t0 = pd.Timestamp(c["t_start"])
|
|
t1 = t0 + pd.Timedelta(seconds=c["dur_s"])
|
|
cloud_mask |= (df["timestamp"] >= t0) & (df["timestamp"] <= t1)
|
|
|
|
A, Bc, C, rms, n_fit, curve, model_peak = fit_clearsky(
|
|
el, am, G, cloud_mask, t.values)
|
|
clipped = G >= G_CLIP
|
|
n_clip = int(clipped.sum())
|
|
dt_med = float(np.median(df["dt_s"]))
|
|
|
|
# correccion monotona: lo saturado toma el modelo (nunca por debajo
|
|
# de los 1024 medidos); lo demas conserva la medida
|
|
G_est = np.where(clipped, np.maximum(1024.0, curve), G)
|
|
|
|
# R_opt y P_ref: con y sin correccion
|
|
lit = G > 600
|
|
R_now = float((Isc[lit] / G[lit]).max())
|
|
P_ref_now = df["P_ref"].values
|
|
litc = G_est > 600
|
|
R_corr = float((Isc[litc] / G_est[litc]).max())
|
|
P_ref_corr = R_corr * G_est * Voc * pl.FF
|
|
|
|
# ventanas identicas a las del articulo (metrics.json)
|
|
m_full = pl.window_metrics(df, "Full (G>50)")
|
|
m_clear = pl.window_metrics(df[df["G_POA"] > 600], "Clear (G>600)")
|
|
df_c = df.copy()
|
|
df_c["G_POA"] = G_est
|
|
df_c["P_ref"] = P_ref_corr
|
|
m_full_c = pl.window_metrics(df_c, "Full (G>50)")
|
|
m_clear_c = pl.window_metrics(df_c[df_c["G_POA"] > 600], "Clear (G>600)")
|
|
|
|
df_mod = df.copy()
|
|
df_mod["P_ref"] = P_ref_corr
|
|
t95_now = pl.recovery_t95(df, df_e)
|
|
t95_corr = pl.recovery_t95(df_mod, df_e)
|
|
|
|
# MAE con la referencia corregida (ventana de cielo claro del articulo)
|
|
Isc_ref_corr = R_corr * G_est
|
|
wc = (G > 600)
|
|
MAE_I_now = float(np.mean(np.abs((R_now * G) - Isc)) * 1000.0)
|
|
MAE_I_corr = float(np.mean(np.abs(Isc_ref_corr - Isc)[wc])) * 1000.0
|
|
MAE_P_now = m_clear["MAE_P_mW"]
|
|
MAE_P_corr = float(np.mean(np.abs(P_ref_corr - Pc)[wc])) * 1000.0
|
|
|
|
print("\n" + "=" * 72)
|
|
print(f" CORRECCION DE SATURACION (1024 W/m2) - {tag}")
|
|
print(f" Aguascalientes, Ags., Mexico (21.8756 N, 102.2854 W, 1896 m)")
|
|
print("=" * 72)
|
|
print(f" Modelo cielo despejado: G_cs = {A:.0f} * exp(-{Bc:.3f} * AM) + {C:.0f}")
|
|
print(f" rms sobre n={n_fit} muestras claras no saturadas: {rms:.0f} W/m2")
|
|
print(f" Pico: medido = {G.max():.0f} W/m2 (tope del ADC) | "
|
|
f"estimado = {G_est.max():.0f} W/m2")
|
|
print(f" Muestras saturadas: {n_clip} "
|
|
f"({100.0 * n_clip / len(df):.1f} % del log, ~{n_clip * dt_med / 60:.0f} min)")
|
|
print(f" H (kWh/m2, G>50): {m_full['H_kWh_m2']} -> "
|
|
f"{m_full_c['H_kWh_m2']} "
|
|
f"(+{100.0 * (m_full_c['H_kWh_m2'] / m_full['H_kWh_m2'] - 1):.1f} %)")
|
|
print(f" R_opt (uA/W/m2): {R_now * 1e6:.0f} -> {R_corr * 1e6:.0f}")
|
|
print(f" eta_tr (G>50): {m_full['eta_tr_pct']} % -> "
|
|
f"{m_full_c['eta_tr_pct']} %")
|
|
print(f" eta_tr (G>600): {m_clear['eta_tr_pct']} % -> "
|
|
f"{m_clear_c['eta_tr_pct']} %")
|
|
print(f" MAE_I (G>600): {MAE_I_now:.1f} mA -> "
|
|
f"{MAE_I_corr:.1f} mA")
|
|
print(f" MAE_P (G>600): {MAE_P_now} mW -> {MAE_P_corr:.1f} mW")
|
|
print(f" med P/Pref (G>600): {m_clear['ratio_P_Pref_median'] if 'ratio_P_Pref_median' in m_clear else float(np.median(Pc[wc] / P_ref_now[wc])):.3f} -> "
|
|
f"{float(np.median(Pc[wc] / P_ref_corr[wc])):.3f}")
|
|
if len(t95_now) == len(t95_corr) and len(t95_now):
|
|
print("\n t95 (s): antes -> corregido")
|
|
for a, b in zip(t95_now.itertuples(), t95_corr.itertuples()):
|
|
print(f" {a.trigger:>14} @ {str(a.t_trigger)[11:19]}: "
|
|
f"{a.t95_s} -> {b.t95_s}")
|
|
|
|
if not nofig:
|
|
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12.5, 4.4),
|
|
sharex=True)
|
|
ax1.plot(t, Pc * 1000, color="tab:red", lw=0.8, label="P_CPV")
|
|
ax1.plot(t, P_ref_now * 1000, color="k", lw=1.1,
|
|
label="P_ref (G medida)")
|
|
ax1.plot(t, P_ref_corr * 1000, color="tab:blue", ls="--", lw=1.1,
|
|
label="P_ref (G corregida)")
|
|
ax1.set_ylabel("Potencia (mW)")
|
|
ax1.set_xlabel("Hora local (UTC-6)")
|
|
ax1.legend(fontsize=8)
|
|
ax1.grid(alpha=0.3)
|
|
|
|
if n_clip:
|
|
idx = np.where(clipped)[0]
|
|
start = 0
|
|
for j in range(1, len(idx) + 1):
|
|
if j == len(idx) or idx[j] != idx[j - 1] + 1:
|
|
ax2.axvspan(t.iloc[idx[start]], t.iloc[idx[j - 1]],
|
|
color="tab:orange", alpha=0.18)
|
|
start = j
|
|
ax2.plot(t, G, color="tab:orange", lw=0.8, label="G medida")
|
|
ax2.plot(t, G_est, color="tab:orange", ls="--", lw=1.3,
|
|
label="G estimada (corregida)")
|
|
ax2.plot(t, curve, color="gray", ls=":", lw=0.8,
|
|
label="Modelo cielo despejado")
|
|
ax2.axhline(1024, color="gray", ls="-.", lw=0.8)
|
|
ax2.text(t.iloc[1], 1030, "saturacion 1024 W/m2",
|
|
fontsize=8, color="dimgray")
|
|
ax2.set_ylabel("Irradiancia global (W/m2)")
|
|
ax2.set_xlabel("Hora local (UTC-6)")
|
|
ax2.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M"))
|
|
ax2.legend(fontsize=8, loc="upper right")
|
|
ax2.grid(alpha=0.3)
|
|
fig.tight_layout()
|
|
|
|
outdir = ROOT / "articulo" / f"resultados_{tag}"
|
|
outdir.mkdir(parents=True, exist_ok=True)
|
|
fig.savefig(outdir / "preview_irradiance_est.png", dpi=150)
|
|
|
|
res = {
|
|
"tag": tag,
|
|
"model": "G_cs = A*exp(-B*AM) + C, Ajustado al sobre-ensvolvente "
|
|
"claro no saturado",
|
|
"model_scale_A_Wm2": round(A, 1), "model_B": round(Bc, 4), "model_C": C,
|
|
"model_rms_clear_Wm2": round(rms, 1), "n_fit_clear": n_fit,
|
|
"model_peak_Wm2": round(model_peak, 0),
|
|
"G_max_measured_Wm2": round(float(G.max()), 1),
|
|
"G_max_estimated_Wm2": round(float(G_est.max()), 1),
|
|
"n_clipped": n_clip,
|
|
"clipped_pct_log": round(100 * n_clip / len(df), 2),
|
|
"H_kWh_m2_full_measured": m_full["H_kWh_m2"],
|
|
"H_kWh_m2_full_corrected": m_full_c["H_kWh_m2"],
|
|
"R_opt_uA_per_Wm2_measured": round(R_now * 1e6, 1),
|
|
"R_opt_uA_per_Wm2_corrected": round(R_corr * 1e6, 1),
|
|
"eta_tr_pct_full_measured": m_full["eta_tr_pct"],
|
|
"eta_tr_pct_full_corrected": m_full_c["eta_tr_pct"],
|
|
"eta_tr_pct_clear_measured": m_clear["eta_tr_pct"],
|
|
"eta_tr_pct_clear_corrected": m_clear_c["eta_tr_pct"],
|
|
"MAE_I_mA_clear_measured": round(MAE_I_now, 2),
|
|
"MAE_I_mA_clear_corrected": round(MAE_I_corr, 2),
|
|
"MAE_P_mW_clear_measured": MAE_P_now,
|
|
"MAE_P_mW_clear_corrected": round(MAE_P_corr, 2),
|
|
"ratio_P_Pref_median_clear_measured":
|
|
round(float(np.median(Pc[wc] / P_ref_now[wc])), 3),
|
|
"ratio_P_Pref_median_clear_corrected":
|
|
round(float(np.median(Pc[wc] / P_ref_corr[wc])), 3),
|
|
"location": "Aguascalientes, Ags., Mexico "
|
|
"(21.8756 N, 102.2854 W, 1896 m, UTC-6)",
|
|
}
|
|
outdir = ROOT / "articulo" / f"resultados_{tag}"
|
|
outdir.mkdir(parents=True, exist_ok=True)
|
|
(outdir / f"irradiance_est_{tag}.json").write_text(
|
|
json.dumps(res, indent=2, ensure_ascii=False), encoding="utf-8")
|
|
print(f"\n[ok] JSON (y figura) en {outdir}")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
tag = sys.argv[1] if len(sys.argv) > 1 else "21092025"
|
|
nofig = len(sys.argv) > 2 and sys.argv[2] == "nofig"
|
|
run(tag, nofig=nofig) |