Viendo la calima de Canarias con satpy en Python

Estos días (principalmente el 22 y 23 de febrero del 2020) hemos visto imágenes que han convertido a las islas Canarias en escenarios propios de Blade runner 2049.

Fotografía de Bengt Nyman. Foto original con licencia CC BY 2.0.

Esto es debido a la entrada de arena del Sahara y pueden llegar a distancias de miles de kilómetros:

Como comprenderás, respirar esto no es nada bueno pero hoy no voy a hablar de esto, principalmente porque no tengo ni idea. Voy a hablar de cómo crear una imagen animada usando imágenes del satélite meteorológico geostacionario Meteosat y Pillow y luego voy a mostrar cómo ver una imagen del satélite polar Sentinel-3 usando satpy.

Meteosat y `Pillow`

El código que voy a mostrar para descargar las imágenes del Meteosat está explicado más en detalle en este artículo. También voy a mostrar una forma alternativa a la que usé en este otro artículo de crear un GIF animado. En lugar de usar imageio voy a usar Pillow. También hable sobre cómo usar Pillow en este otro artículo. En las referencias más abajo vuelvo a dejar todos estos enlaces.

Bueno, a lo que iba, que me disperso. Pongo el código y comentaré la parte de Pillow, que es la novedosa.

import datetime as dt
from io import BytesIO 

from owslib.wms import WebMapService 
from PIL import Image

# Parámetros para descargar imágenes de un Web Map Service
wms = WebMapService(
    "https://eumetview.eumetsat.int/geoserv/wms", 
    version="1.3.0"
)
wms_params = { 
    "bbox": (-24, 24, -5, 35), 
    "size": (690, 400), 
    "srs": "EPSG:4326", 
    "layers": ["meteosat:msg_natural"], 
    "format": "image/jpeg"
}

# Descargamos una imagen cada 15 minutos y las guardamos en la lista `imgs`
# La lista `imgs` contiene objetos PIL.Image.Image
tdelta = dt.timedelta(minutes=15)
datetimes = [dt.datetime(2020, 2, 22, 7, 0) + tdelta * i for i in range(48)]
imgs = [] 
for datetime in datetimes: 
    wms_params["time"] = datetime.isoformat() + "Z" # ugly hack 
    imgfile = BytesIO()
    imgfile.write(wms.getmap(**wms_params).read())
    imgs.append(Image.open(imgfile))

# El objeto PIL.Image.Image tiene un método save que uso aquí
imgs[0].save(
    'msg_calima_canarias_20200222.gif',
    save_all=True,
    append_images=imgs[1:48],
    duration=100,
    loop=0
)

import datetime as dt

from io import BytesIO

from owslib.wms import WebMapService

from PIL import Image

# Parámetros para descargar imágenes de un Web Map Service

wms = WebMapService(

"https://eumetview.eumetsat.int/geoserv/wms",

version="1.3.0"

)

wms_params = {

"bbox": (-24, 24, -5, 35),

"size": (690, 400),

"srs": "EPSG:4326",

"layers": ["meteosat:msg_natural"],

"format": "image/jpeg"

}

# Descargamos una imagen cada 15 minutos y las guardamos en la lista `imgs`

# La lista `imgs` contiene objetos PIL.Image.Image

tdelta = dt.timedelta(minutes=15)

datetimes = [dt.datetime(2020, 2, 22, 7, 0) + tdelta * i for i in range(48)]

imgs = []

for datetime in datetimes:

wms_params["time"] = datetime.isoformat() + "Z" # ugly hack

imgfile = BytesIO()

imgfile.write(wms.getmap(**wms_params).read())

imgs.append(Image.open(imgfile))

# El objeto PIL.Image.Image tiene un método save que uso aquí

imgs[0].save(

'msg_calima_canarias_20200222.gif',

save_all=True,

append_images=imgs[1:48],

duration=100,

loop=0

)

La línea con imgs[0].save(...) es la que hace la magia para crear el GIF animado. Guardamos la imagen inicial usando el nombre msg_calima_canarias_20200222.gif. La opción save_all se usa para que se incluyan todos los frames al guardar. La opción append_images se usa para indicar qué imágenes incluir. La opción duration se usa para indicar la duración en que se muestra cada imagen. La opción loop permite definir si el loop se ve una sola vez, dos,…, o infinitas (usando 0),

El resultado lo puedes ver aquí:

Se ven cosas interesantes, la lengua de arena “comiéndose” las islas Canarias, convección, como algunas islas actúan de obstáculo (mirad cuando la lengua llega a Gran Canaria).

Sentinel-3 y `satpy`

El Meteosat, como está en órbita geoestacionaria, puede ver el mismo punto de la tierra todo el rato. Los satélites con órbita baja (LEO, Low Earth Orbit) están dando vueltas a la tierra y van haciendo barridos de la misma. En este caso uso una imagen del instrumento OLCI (Ocean and Land Colour Instrument) a bordo del satélite Sentinel-3.

Para acceder a las imágenes de Sentinel-3 el proceso es un poco más farragoso. Debes ir al servicio Copernicus Online Data Access (CODA), registrarte, si no tienes cuenta, y buscar la imagen y descargarla. Aquí te explican cómo (en inglés).

Para este ejemplo he descargado el fichero S3A_OL_1_EFR____20200222T110740_20200222T111040_20200222T131218_0179_055_151_2520_MAR_O_NR_002.SEN3. En realidad no es un fichero, es una carpeta con ficheros con formatos NetCDF y xml y pesa más de 600 Mb. Como esto no es un tutorial de satpy (quizá llegue uno en el futuro) solo voy a mostrar un código mínimo para crear la imagen (que además he copiado de aquí):

from datetime import datetime

from satpy import Scene, find_files_and_readers
files = find_files_and_readers(start_time=datetime(2020, 2, 22, 11, 7),
                               end_time=datetime(2020, 2, 22, 11, 10),
                               base_dir="./S3-data/OLCI",
                               reader="olci_l1b")
scn = Scene(filenames=files)
composite = 'true_color'
scn.load([composite])
scn.show(composite)
scn.save_dataset(composite, filename="sand.png")

from datetime import datetime

from satpy import Scene, find_files_and_readers

files = find_files_and_readers(start_time=datetime(2020, 2, 22, 11, 7),

end_time=datetime(2020, 2, 22, 11, 10),

base_dir="./S3-data/OLCI",

reader="olci_l1b")

scn = Scene(filenames=files)

composite = 'true_color'

scn.load([composite])

scn.show(composite)

scn.save_dataset(composite, filename="sand.png")

El anterior código crea una imagen que se llama sand.jpg que es la que se muestra a continuación:

La imagen es de las 11 UTC por lo que solo había llegado a las islas más orientales (Fuerteventura, Lanzarote, la Graciosa,…). Esta imagen tiene mucha más resolución espacial que el Meteosat por lo que se pueden ver más detalles. En las imágenes de abajo hago zoom sobre Lanzarote para resaltar algunos fenómenos físicos interesantes:

¿Qué veis ahí?

Al Nordeste de la Graciosa hay una isla que se llama Montaña Clara y tiene casi tanto de alto como de ancho. Debajo puedes ver una imagen extraida de Earth;

Para un flujo de aire puede resultar más fácil rodear la isla por los laterales que ascender. Se puede ver que en la imagen con zoom hay una estela (como una calle con menos arena) detrás de esta isla. También se pueden ver una especie de bandas a sotavento de Lanzarote. Esto podrían ser ondas de gravedad (no confundir con ondas gravitacionales) generadas por la propia isla (otro ejemplo de ondas de gravedad a sotavento en Lanzarote aquí).

Lo que comento está destacado en la imagen de abajo:

PRECIOSO, ¿VERDAD?

Referencias

Más imágenes interesantes vistas por satélites de la NOAA: https://www.nesdis.noaa.gov/content/noaa-20-watches-sandstorm-blow-over-canary-islands
Predicción de polvo del BSC (que lo hizo bastante bien durante este episodio): https://ess.bsc.es/bsc-dust-daily-forecast
Predicción de polvo del BSC desde la web de Aemet: https://dust.aemet.es/
Accediendo a Web Map Services con Python: https://pybonacci.org/2018/11/27/creacion-de-videos-al-vuelo-en-el-notebook/
Creando GIFs animados con imageio: https://pybonacci.org/2019/12/13/creacion-de-gif-animado-con-imageio-y-python/
Ejemplo de uso de Pillow: https://pybonacci.org/2019/12/31/crea-un-mosaico-de-imagenes-usando-pillow-pil-en-python/

2 comentarios en «Viendo la calima de Canarias con satpy en Python»

Pablo el 2020-02-27 a las 10:55

Muy interesante y didáctico. Gracias por publicar post como este.

Responder
1. Kiko el 2020-02-27 a las 21:44
  
  Muchas gracias por comentar. Me alegro mucho que te haya sido útil e interesante.
  
  Responder

Viendo la calima de Canarias con satpy en Python

Meteosat y Pillow

Sentinel-3 y satpy

Referencias

2 comentarios en «Viendo la calima de Canarias con satpy en Python»

Deja una respuesta Cancelar la respuesta

Meteosat y `Pillow`

Sentinel-3 y `satpy`