Acceso a pronósticos GFS en la nube

from arraylake import Client
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from datetime import datetime, timedelta

client = Client()
client.login()  # Solo necesario la primera vez por sesión

repo = client.get_repo("earthmover-public/gfs")

# Establecer vista de solo lectura del branch "main" inmutable
session = repo.readonly_session("main")

%%time
ds_gfs = xr.open_dataset(
    session.store,
    engine="zarr",
    consolidated=False,
    zarr_format=3,
    chunks={},
    group="solar"  # Grupo con variables de superficie
)

display(ds_gfs)

# Calcular fecha de mañana
tomorrow = datetime.now() + timedelta(days=1)
fecha_str = f"{tomorrow:%Y-%m-%d}"

print(f"Pronóstico para: {fecha_str}")

# Seleccionar temperatura para mañana al mediodía
t2m_forecast = ds_gfs.t2m.sel(time=f"{fecha_str} 12:00", method="nearest")

# Convertir de Kelvin a Celsius
t2m_celsius = t2m_forecast - 273.15
t2m_celsius

# Creamos una visualizacion 
fig, ax = plt.subplots(figsize=(12, 6))
t2m_celsius.isel(step=10).plot(ax=ax, x="longitude", y="latitude", cmap="RdYlBu_r", robust=True)
ax.set_title(f"Temperatura a 2m - Pronóstico GFS ({fecha_str} 12:00 UTC)", fontsize=14)
ax.set_xlabel("Longitud [°E]")
ax.set_ylabel("Latitud [°N]")
plt.tight_layout()
plt.show()

# Seleccionar región de Colombia
# Longitud: 275° a 300° (equivalente a -85° a -60°)
# Latitud: 14° a -5° (Norte a Sur)
gfs_colombia = ds_gfs.sel(
    longitude=slice(275, 300),
    latitude=slice(14, -5)).sel(
    time=f"{fecha_str} 12:00",
    method="nearest"
)

print(f"Shape Colombia: {gfs_colombia.t2m.shape}")

gfs_colombia

# Crear mapa para Colombia con proyección PlateCarree
fig, ax = plt.subplots(
    # figsize=(10, 8),
    subplot_kw={"projection": ccrs.PlateCarree()},
    dpi=150
)

# Temperatura en Celsius
t2m_col = gfs_colombia.t2m - 273.15

# Graficar temperatura
t2m_col.isel(step=0).plot.contourf(
    ax=ax,
    x="longitude", 
    y="latitude",
    cmap="RdYlBu_r",
    transform=ccrs.PlateCarree(),
    cbar_kwargs={
        "label": "Temperatura [°C]",
        "orientation": "horizontal",
        "shrink": 0.8
    }
)

# Añadir características geográficas
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.gridlines(draw_labels=True, linewidth=0.5, alpha=0.5)

ax.set_title(f"Pronóstico GFS - Colombia ({fecha_str} 12:00 UTC)", fontsize=12)
plt.show()

# Crear panel con 4 variables
fig, axes = plt.subplots(
    2, 2,
    figsize=(8, 6),
    subplot_kw={"projection": ccrs.PlateCarree()},
    dpi=120
)
axes = axes.flatten()

# Variable 1: Temperatura
(gfs_colombia.t2m.isel(step=0) - 273.15).plot(
    x="longitude", y="latitude",
    ax=axes[0], cmap="RdYlBu_r", transform=ccrs.PlateCarree(),
    cbar_kwargs={"label": "Temperatura [°C]", "shrink": 0.8}
)
axes[0].set_title("Temperatura 2m")

# Variable 2: Humedad relativa
gfs_colombia.r2.isel(step=0).plot(
    x="longitude", y="latitude",
    ax=axes[1], cmap="YlGnBu", transform=ccrs.PlateCarree(),
    cbar_kwargs={"label": "Humedad relativa [%]", "shrink": 0.8}
)
axes[1].set_title("Humedad relativa 2m")

# Variable 3: Cobertura de nubes
(gfs_colombia.isel(step=0).tcc * 100).plot(
    x="longitude", y="latitude",
    ax=axes[2], cmap="gray_r", transform=ccrs.PlateCarree(),
    cbar_kwargs={"label": "Cobertura [%]", "shrink": 0.8}
)
axes[2].set_title("Cobertura de nubes total")

# Variable 4: Precipitación
(gfs_colombia.isel(step=0).prate * 86400).plot(  # Convertir a mm/día
    x="longitude", y="latitude",
    ax=axes[3], cmap="Blues", transform=ccrs.PlateCarree(),
    vmin=0, vmax=100,
    cbar_kwargs={"label": "Precipitación [mm/día]", "shrink": 0.8}
)
axes[3].set_title("Tasa de precipitación")

# Añadir características a todos los mapas
for ax in axes:
    ax.coastlines()
    ax.add_feature(cfeature.BORDERS, linestyle=":")
    ax.gridlines(alpha=0.3)

plt.suptitle(f"Pronóstico GFS - Colombia ({fecha_str} 12:00 UTC)", fontsize=16, y=0.98)
plt.tight_layout()
plt.show()

# Coordenadas de Bogotá: 4.7°N, 285.9°E (-74.1°W)
today = datetime.now()
week_ahead = timedelta(hours=168)

t2m_bogota = ds_gfs.t2m.sel(
    latitude=4.7,
    longitude=285.9,
    time=f"{today}",
    method="nearest"
).sel(step=slice(None, week_ahead))

# Convertir a Celsius
t2m_bogota_celsius = t2m_bogota - 273.15

print(f"Pronóstico para Bogotá: {len(t2m_bogota_celsius)} pasos de tiempo")

t2m_bogota

fig, ax = plt.subplots(figsize=(12, 5))

t2m_bogota_celsius.plot(ax=ax, linewidth=1.5, color="steelblue")

ax.set_title("Pronóstico de temperatura - Bogotá (próximos 7 días)", fontsize=14)
ax.set_xlabel("Fecha")
ax.set_ylabel("Temperatura [°C]")
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Tu código aquí

# Tarea 1: Mapa de ráfagas de viento para Colombia
fig, ax = plt.subplots(
    figsize=(10, 8),
    subplot_kw={'projection': ccrs.PlateCarree()},
    dpi=120
)

# Seleccionar ráfagas de viento
gust_col = gfs_colombia.gust.isel(step=0)

# Graficar
gust_col.plot(
    x='longitude',
    y='latitude',
    ax=ax,
    cmap='YlOrRd',
    transform=ccrs.PlateCarree(),
    cbar_kwargs={'label': 'Ráfagas de viento [m/s]', 'shrink': 0.8}
)

ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.gridlines(draw_labels=True, linewidth=0.5, alpha=0.5)
ax.set_title(f'Pronóstico de ráfagas de viento - Colombia ({fecha_str} 12:00 UTC)')
plt.show()

# Tarea 2: Precipitación acumulada de 24 horas
# Seleccionar las primeras 24 horas de pronóstico
precip_24h = gfs_colombia.prate.isel(step=slice(0, 24))

# Convertir de kg/m²/s a mm/día y sumar sobre 24 horas
# 1 kg/m²/s * 3600 s/h = mm/h, luego sumar 24 horas
precip_acum_24h = (precip_24h * 3600).sum(dim='step')

# Visualizar
fig, ax = plt.subplots(
    figsize=(10, 8),
    subplot_kw={'projection': ccrs.PlateCarree()},
    dpi=120
)

precip_acum_24h.plot(
    x='longitude',
    y='latitude',
    ax=ax,
    cmap='Blues',
    transform=ccrs.PlateCarree(),
    cbar_kwargs={'label': 'Precipitación acumulada [mm]', 'shrink': 0.8}
)

ax.coastlines()
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.gridlines(draw_labels=True, linewidth=0.5, alpha=0.5)
ax.set_title(f'Precipitación acumulada 24h - Colombia ({fecha_str})')
plt.show()

print(f'Precipitación máxima en 24h: {float(precip_acum_24h.max().values):.1f} mm')

# Tarea 3: Evolución temporal de nubosidad para Medellín
# Coordenadas: lat=6.25, lon=284.5 (equivalente a -75.5°W)
week_ahead = timedelta(hours=168)

tcc_medellin = ds_gfs.tcc.sel(
    latitude=6.25,
    longitude=284.5,
    time=f"{today}",
    method='nearest'
).sel(step=slice(None, week_ahead))

# Convertir a porcentaje
tcc_medellin_pct = tcc_medellin * 100

# Visualizar
fig, ax = plt.subplots(figsize=(12, 5))

tcc_medellin_pct.plot(ax=ax, linewidth=1.5, color='gray')

ax.set_title('Pronóstico de cobertura de nubes - Medellín (próximos 7 días)', fontsize=14)
ax.set_xlabel('Paso de tiempo')
ax.set_ylabel('Cobertura de nubes [%]')
ax.set_ylim(0, 100)
ax.grid(True, alpha=0.3)
ax.axhline(y=50, color='orange', linestyle='--', alpha=0.5, label='50% cobertura')
ax.legend()
plt.tight_layout()
plt.show()

print(f'Cobertura promedio: {float(tcc_medellin_pct.mean().values):.1f}%')
print(f'Cobertura máxima: {float(tcc_medellin_pct.max().values):.1f}%')
print(f'Cobertura mínima: {float(tcc_medellin_pct.min().values):.1f}%')