Creando un mosaico usando `s2froms3` en Python a partir de imágenes de Sentinel-2

He creado un pequeño paquete para poder descargar Cloud Optimized Geotiffs (COGs) del satélite Sentinel-2 desde S3. El satélite Sentinel-2 es un satélite operado por la ESA (European Space Agency) dentro del programa Copernicus de la UE. El programa Copernicus es un ambicioso proyecto para monitorizar la Tierra desde el espacio. Las imágenes de Sentinel-2 son la opción más sencilla para poder tener imágenes de alta resolución del planeta obtenida a partir de distintos canales. Con un poco de código puedes descargar Gigabytes de datos de los distintos canales del satélite Sentinel-2.

Voy a empezar importando algunas cosas que necesitaré. Si quieres seguir el tutorial tendrás que instalar s2froms3, xarray, rasterio, pyepsg, cartopy, matplotlib. Si tienes algún problema coméntalo e intentaremos solucionarlo.

Como es habitual, me voy a centrar en una localización conocida. La defino usando su latitud y su longitud:

Como las imágenes (COGs) vienen en teselas (divididas siguiendo MGRS) puede que mi localización quede en una esquina de la tesela que voy a descargar y, al final, la zona de interés no esté centrada. Para saber donde quedaría localizada la localización dentro de la tesela a descargar he incluido una función de ayuda dentro de s2froms3 llamada point_in_tile. Se puede usar así:

El resultado de lo anterior mostraría:

Se ve que mi punto de interés está muy escorado hacia el este. Si mi localización indica el centro de lo que quiero mostrar sé que parte hacia el este seguramente quedará cortado. Luego lo tendremos en cuenta.

Como he comentado más arriba, podemos obtener diferentes productos del satélite Sentinel-2. Para saber la información que tenemos disponible para descargar podemos usar la información que he incluido en s2froms3.products.Properties:

Lo anterior mostraría:

En lo anterior veis que disponemos de varias bandas o productos derivados como el TCI (True Color). Se puede ver donde se centra la longitud de onda de la banda, la resolución de la imagen (en metros) y algo más de información. Dependiendo de lo que queramos analizar podemos usar unas bandas u otras o combinarlas de forma inteligente para obtener productos derivados. Las imágenes TCI se obtienen a partir de las bandas 4, 3 y 2, que serán las bandas espectrales del rojo, del verde y del azul, respectivamente, y se usan para obtener una imagen que vendría a representar como vemos los humanos (visible).

Tenemos, más o menos, todos los ingredientes que necesitamos. ¿Vamos a por los COGs? Para descargar los COGs podemos usar la función download_S2. El satélite Sentinel-2 es un satélite de órbita polar que va barriendo la tierra en pasadas de Sur a Norte y de Norte a Sur. Por tanto, no está mirando siempre sobre el mismo punto y puede que pasen varios días hasta que vuelva a pasar sobre la misma localización. Hemos de definir las fechas para descargar las imágenes, si el satélite no pasó sobre la localización durante esas fechas pues no habrá imágenes. Hemos de definir qué queremos descargar, qué bandas o productos. Hemos de definir la cubierta nubosa máxima admisible que queremos en la imagen ya que si hay nubes alguna información igual nos resulta menos útil. Etcétera.

Vamos con un ejemplo durante el verano boreal (estamos en el hemisferio Norte), periodo en el que suele haber menos nubes, y durante unas fechas que sé que ha pasado Sentinel-2 sobre mi área de interés. Por mantener esto sencillo voy a descargar TCI (imágenes True Color).

Definimos lo que necesitamos:

Con la anterior información y con la localización definida más arriba puedo proceder a descargar las imágenes:

Con la opción folder defino donde quiero que se descarguen las imágenes. Con la opción also indico si, además, quiero descargar imágenes adyacentes a la tesela que estoy descargando. Le indico que sí, que quiero, además de la tesela que contiene a lon y lat también quiero las imágenes que se localizan en las teselas al Este, al Noreste y al Norte.

La función download_S2 habrá descargado las imágenes y devuelve las rutas de esas imágenes en S3:

Lo anterior mostrará:

Vamos a leer las imágenes que hemos descargado usando xarray (y rasterio entre bambalinas):

Lo anterior mostrará en pantalla lo siguiente:

Tenemos 4 datasets con la misma forma (3 bandas y el mismo tamaño en x e y. Vamos a visualizar los datasets usando cartopy y matplotlib:

El resultado mostrará:

Tenemos cuatro piezas que unir para resolver nuestro pequeño rompecabezas. Vamos a seguir la misma técnica que antes usando cartopy y matplotlib pero modificando algunas pequeñas cosas.

El crs (el sistema de referencia) se puede ver que es el mismo en todos los casos. Como quiero unir las 4 piezas necesito una imagen que incluya a las 4. Para ello necesito saber las localizaciones mínimas y máximas del conjunto de las 4 imágenes (x0, x1, y0, y1 en el código de más abajo). Y necesito saber la referencia de cada una de las imágenes individuales en la imagen común que estamos creando (x_min, x_max, y_min, y_max en el código de más abajo).

Con lo anterior podemos obtener la imagen final:

¡¡¡Qué isla tan bonita!!!

Pero, ¿qué son esas bandas que se ven por ahí? Puedes leer un poco sobre ello en la sección ‘A few words on the post-processing‘ en este artículo de Pierre Markuse. Corregir eso con Python daría para otro artículo. ¿Te animas a escribirlo?

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

twenty seven − twenty six =