Día y noche, mezclando visible y otros canales de satélite con xarray y Python

Siempre me fascina ver nubes, desde abajo y desde arriba. Desde arriba las veo usando satélites ya que, de momento, no he podido viajar al espacio 🙁

Los satélites nos muestran las nubes de varias formas, usando distintos canales para distintas bandas espectrales. Juntando varios de estos canales podemos crear productos interesantes. En el caso del Meteosat Second Generation se puede crear un producto parecido a si sacásemos una fotografía desde el espacio juntando los canales NIR1.6, VIS0.8 y VIS0.6.

Son las imágenes más bonitas, en mi modesta opinión, pero solo se pueden obtener de día.

¿Pero yo he visto que a veces retuiteas vídeos de ciclones/huracanes/etc y salen imágenes durante más de 24h?, ¿cómo lo hacen? Pues se hace una mezcla de imágenes visibles durante el día y por la noche se mezclan con imágenes de otros canales. Normalmente, se usa algún canal infrarrojo que nos da información de la temperatura. Las nubes -el tope nuboso, más bien-, al estar en altura, suelen estar más frías y contra más altas, normalmente, más frías.

Pues me he liado la manta a la cabeza y he intentado replicar lo que muestra, por ejemplo, @weatherdak en tuiter.

¿Vamos a ello?

Como siempre, empezamos importando algunas cosas:

import requests
import datetime as dt
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import rioxarray as rxr
import ffmpeg

Una de las gemas aquí es rioxarray. Hace que rasterio y xarray se pongan, incluso, más interesantes añadiendo funcionalidad extra.

Lo siguiente es una pequeña función que genera un rango de fechas que nos servirá para generar las URLs para descargar las imágenes que usaremos.

def get_datetime_range():
    end = dt.datetime.utcnow()
    end -= dt.timedelta(hours=1)
    mm = (end.minute // 15) * 15
    end = dt.datetime(end.year, end.month, end.day, end.hour, mm)
    start = end - dt.timedelta(days=2)
    date_range = []
    while start <= end:
        date_range.append(start)
        start += dt.timedelta(minutes=15)
    return date_range

El diccionario siguiente me sirve para seleccionar lo que quiero descargar. En realidad solo usaré el IR10.8 y el producto Natural (que viene a ser como una fotografía de la Tierra):

layers = {
    'IR39': 'meteosat%3Amsg_ir039',
    'IR108': 'meteosat%3Amsg_ir108',
    'VIS06': 'meteosat%3Amsg_vis006',
    'WV062': 'meteosat%3Amsg_wv062',
    'Natural': 'meteosat%3Amsg_natural',
    'NaturalEnhanced': 'meteosat%3Amsg_naturalenhncd'
}

El siguiente código es el que descarga las imágenes. La llamada a get_datetime_range me genera un listado de fechas para descargar las imágenes (cada 15 minutos). Como he comentado anteriormente, descargamos imágenes de dos productos. ATENCIÓN, EL SIGUIENTE CÓDIGO DESCARGARÁ UNOS CUANTOS GIGABYTES DE DATOS. Las fechas de lo siguiente corresponden al día en que hice esto:

date_range = get_datetime_range()
keys = ['Natural', 'IR108']

for _dt in date_range:
    for key in keys:
        print(_dt)
        layer = layers[key]
        url = (
            'https://eumetview.eumetsat.int/geoserv/wms?SERVICE=WMS&REQUEST=GetMap&'
            'TRANSPARENT=TRUE&EXCEPTIONS=INIMAGE&VERSION=1.3.0&'
            f'LAYERS={layer}'
            '%2Coverlay%3Ane_10m_coastline%2Coverlay%3Ane_10m_admin_0_boundary_lines_land&'
            'STYLES=raster%2C%2C&SRS=EPSG%3A4326&WIDTH=1400&HEIGHT=1100&'
            'BBOX=26.0,-20.0,48.0,8.0&'
            'FORMAT=image%2Fgeotiff&'
            f'TIME={_dt:%Y-%m-%d}T{_dt:%H:%M}:00.000Z&'
        )
        with requests.get(url) as r:
            with open(f'{key}_{_dt:%Y%m%d_%H%M}.tif', 'wb') as f:
                f.write(r.content)

Si abro la imagen Natural_20210501_0600.tif, una imagen que usa datos de canales visibles, verás que parte de la imagen no tiene información:

da = xr.open_rasterio('Natural_20210501_0600.tif')
fig, ax = plt.subplots()
ax.imshow(np.moveaxis(da.values, 0, -1))
fig.show()

Lo anterior mostrará:

Lo que vamos a hacer es rellenar esas blancos de la imagen de color natural (la visible) usando información del canal Infrarrojo 10.8.

Con el siguiente código mezclaremos las dos buscando cuando en la imagen visible se cumple cierta condición. Todos los píxeles con esa condición serán rellenados con píxeles de la imagen Infrarroja. Por último, guardamos las imágenes llamándolas merged…:

files_ir = sorted(list(Path('.').glob('IR*.tif')))
files_nat = sorted(list(Path('.').glob('Natural*.tif')))
for fir, fnat in zip(files_ir, files_nat):
    da_ir = rxr.open_rasterio(fir)
    da_nat = rxr.open_rasterio(fnat)
    mask = da_nat.sel(band=4) == 0
    da = da_nat.where(~mask, da_ir)
    part1 = 'merged'
    part2 = str(fir).split('_')[1]
    part3 = str(fir).split('_')[2]
    name = f'{part1}_{part2}_{part3}'
    da.rio.to_raster(name)

Y, por último, vamos a hacer un vídeo de esos que os gustan. Pero esta vez voy a usar ffmpeg-python por probar una cosa distinta.

(
    ffmpeg
    .input('./merged*.tif', pattern_type='glob', framerate=25)
    .output('msg_natural_ir104.mp4')
    .run()
)

El resultado debería ser algo parecido a lo siguiente: