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.

import datetime as dt
from pathlib import Path

import xarray as xr
import s2froms3
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

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

lat, lon = 39.14151, 2.94708

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í:

print(s2froms3.point_in_tile(lon, lat))

El resultado de lo anterior mostraría:

+ + + + + + + + + + + 
+                   + 
+                   + 
+                   + 
+                   + 
+                   + 
+                   + 
+                 O + 
+                   + 
+                   + 
+ + + + + + + + + + + 

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:

for e in s2froms3.products.Properties:
    print(e.name)
    for k, v in e.describe().items():
        print(f"\t{k}: {v}")

Lo anterior mostraría:

TCI
	resolution: 10
	title: True color image
	center wavelength: None
B01
	resolution: 60
	title: Band 1 (coastal)
	center wavelength: 0.4439
B02
	resolution: 10
	title: Band 2 (blue)
	center wavelength: 0.4966
B03
	resolution: 10
	title: Band 3 (green)
	center wavelength: 0.56
B04
	resolution: 10
	title: Band 4 (red)
	center wavelength: 0.6645
B05
	resolution: 20
	title: Band 5
	center wavelength: 0.7039
B06
	resolution: 20
	title: Band 6
	center wavelength: 0.7402
B07
	resolution: 20
	title: Band 7
	center wavelength: 0.7825
B08
	resolution: 10
	title: Band 8 (nir)
	center wavelength: 0.8351
B8A
	resolution: 20
	title: Band 8A
	center wavelength: 0.8648
B09
	resolution: 60
	title: Band 9
	center wavelength: 0.945
B11
	resolution: 20
	title: Band 11 (swir16)
	center wavelength: 1.6137
B12
	resolution: 20
	title: Band 12 (swir22)
	center wavelength: 2.22024
AOT
	resolution: 60
	title: Aerosol Optical Thickness (AOT)
	center wavelength: None
WVP
	resolution: 10
	title: Water Vapour (WVP)
	center wavelength: None
SCL
	resolution: 20
	title: Scene Classification Map (SCL)
	center wavelength: None

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:

start_date = dt.date(2020, 7, 1) # Start date to search images
end_date = dt.date(2020, 7, 1) # End date to search images
what = ['TCI'] # What we want to download
cc = 25 # Minimum cloud cover on each image, 25 is 25%
folder = Path('.') # Where the files will be downloaded

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

downloaded = s2froms3.download_S2(
    lon=lon,
    lat=lat,
    start_date=start_date,
    end_date=end_date,
    what=what,
    cloud_cover_le=cc,
    folder=folder,
    also=['E', 'NE', 'N'],
)

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:

print(downloaded)

Lo anterior mostrará:

['sentinel-cogs/sentinel-s2-l2a-cogs/31/S/DD/2020/7/S2B_31SDD_20200701_0_L2A/TCI.tif', 'sentinel-cogs/sentinel-s2-l2a-cogs/31/S/ED/2020/7/S2B_31SED_20200701_0_L2A/TCI.tif', 'sentinel-cogs/sentinel-s2-l2a-cogs/31/T/DE/2020/7/S2B_31TDE_20200701_0_L2A/TCI.tif', 'sentinel-cogs/sentinel-s2-l2a-cogs/31/T/EE/2020/7/S2B_31TEE_20200701_0_L2A/TCI.tif']

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

files = folder.glob("S2*.tif")

datasets = []

for f in files:
    datasets.append(xr.open_rasterio(f))

print(datasets)

Lo anterior mostrará en pantalla lo siguiente:

[<xarray.DataArray (band: 3, y: 10980, x: 10980)>
[361681200 values with dtype=uint8]
Coordinates:
  * band     (band) int64 1 2 3
  * y        (y) float64 4.5e+06 4.5e+06 4.5e+06 ... 4.39e+06 4.39e+06 4.39e+06
  * x        (x) float64 5e+05 5e+05 5e+05 ... 6.098e+05 6.098e+05 6.098e+05
Attributes:
    transform:           (10.0, 0.0, 499980.0, 0.0, -10.0, 4500000.0)
    crs:                 +init=epsg:32631
    res:                 (10.0, 10.0)
    is_tiled:            1
    nodatavals:          (0.0, 0.0, 0.0)
    scales:              (1.0, 1.0, 1.0)
    offsets:             (0.0, 0.0, 0.0)
    AREA_OR_POINT:       Area
    OVR_RESAMPLING_ALG:  AVERAGE, <xarray.DataArray (band: 3, y: 10980, x: 10980)>
[361681200 values with dtype=uint8]
Coordinates:
  * band     (band) int64 1 2 3
  * y        (y) float64 4.4e+06 4.4e+06 4.4e+06 ... 4.29e+06 4.29e+06 4.29e+06
  * x        (x) float64 4e+05 4e+05 4e+05 ... 5.097e+05 5.097e+05 5.098e+05
Attributes:
    transform:           (10.0, 0.0, 399960.0, 0.0, -10.0, 4400040.0)
    crs:                 +init=epsg:32631
    res:                 (10.0, 10.0)
    is_tiled:            1
    nodatavals:          (0.0, 0.0, 0.0)
    scales:              (1.0, 1.0, 1.0)
    offsets:             (0.0, 0.0, 0.0)
    AREA_OR_POINT:       Area
    OVR_RESAMPLING_ALG:  AVERAGE, <xarray.DataArray (band: 3, y: 10980, x: 10980)>
[361681200 values with dtype=uint8]
Coordinates:
  * band     (band) int64 1 2 3
  * y        (y) float64 4.4e+06 4.4e+06 4.4e+06 ... 4.29e+06 4.29e+06 4.29e+06
  * x        (x) float64 5e+05 5e+05 5e+05 ... 6.098e+05 6.098e+05 6.098e+05
Attributes:
    transform:           (10.0, 0.0, 499980.0, 0.0, -10.0, 4400040.0)
    crs:                 +init=epsg:32631
    res:                 (10.0, 10.0)
    is_tiled:            1
    nodatavals:          (0.0, 0.0, 0.0)
    scales:              (1.0, 1.0, 1.0)
    offsets:             (0.0, 0.0, 0.0)
    AREA_OR_POINT:       Area
    OVR_RESAMPLING_ALG:  AVERAGE, <xarray.DataArray (band: 3, y: 10980, x: 10980)>
[361681200 values with dtype=uint8]
Coordinates:
  * band     (band) int64 1 2 3
  * y        (y) float64 4.5e+06 4.5e+06 4.5e+06 ... 4.39e+06 4.39e+06 4.39e+06
  * x        (x) float64 4e+05 4e+05 4e+05 ... 5.097e+05 5.097e+05 5.098e+05
Attributes:
    transform:           (10.0, 0.0, 399960.0, 0.0, -10.0, 4500000.0)
    crs:                 +init=epsg:32631
    res:                 (10.0, 10.0)
    is_tiled:            1
    nodatavals:          (0.0, 0.0, 0.0)
    scales:              (1.0, 1.0, 1.0)
    offsets:             (0.0, 0.0, 0.0)
    AREA_OR_POINT:       Area
    OVR_RESAMPLING_ALG:  AVERAGE]

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:

for ds in datasets:
    crs = ds.crs
    proj = ccrs.epsg(crs.split(':')[-1])

    fig, ax = plt.subplots(subplot_kw={'projection': proj})
    ax.imshow(ds.data.transpose(1, 2, 0), extent=extent)
    fig.show()

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:

crs = datasets[0].crs
proj = ccrs.epsg(crs.split(':')[-1])

x0 = min([ds.x.values.min() for ds in datasets])
x1 = max([ds.x.values.max() for ds in datasets])
y0 = min([ds.y.values.min() for ds in datasets])
y1 = max([ds.y.values.max() for ds in datasets])

fig, ax = plt.subplots(figsize=(20, 20), subplot_kw={'projection': proj})
ax.set_xlim(x0, x1)
ax.set_ylim(y0, y1)

for ds in datasets:
    x_min = ds.x.values.min()
    x_max = ds.x.values.max()
    y_min = ds.y.values.min()
    y_max = ds.y.values.max()
    extent = (x_min, x_max, y_min, y_max)
    ax.imshow(ds.data.transpose(1, 2, 0), extent=extent)
fig.show()

¡¡¡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?