Inmersión profunda en la API del Sentinel de la ESA

Profundizando en la API del Sentinel de la ESA

Cómo adquirir, analizar y visualizar imágenes satelitales usando Python

Todas las imágenes de este artículo fueron creadas por el autor.

La Agencia Espacial Europea ha estado llevando a cabo sus misiones de Sentinel en apoyo a Copernicus, el componente de Observación de la Tierra del programa espacial de la Unión Europea, con diversos tipos de datos de teledetección, como radar e instrumentos de imágenes multiespectrales para el monitoreo de la tierra, el océano y la atmósfera.

Existen seis misiones de Sentinel en curso, tres de las cuales se pueden acceder fácilmente a través de su API de Python. Estas son, citando la fuente oficial:

Sentinel-1 es una misión de imágenes de radar polar, todo tiempo, día y noche, para servicios de tierra y océano. Sentinel-1A fue lanzado el 3 de abril de 2014, y Sentinel-1B el 25 de abril de 2016. Ambos fueron llevados a órbita en un cohete Soyuz desde el Puerto Espacial de Europa en la Guayana Francesa. La misión de Sentinel-1B finalizó en 2022 y hay planes en marcha para lanzar Sentinel-1C lo antes posible.

• Análisis de proximidad para encontrar el bar más cercano utilizando Python
• Análisis estratégico de datos (Parte 3) Preguntas de diagnóstico
• Simplificando la integración y la automatización el CTO de Boomi comparte su visión para los desarrolladores

Sentinel-2 es una misión de imágenes multiespectrales de alta resolución, en órbita polar, para el monitoreo de la tierra, que proporciona, por ejemplo, imágenes de vegetación, coberturas de suelo y agua, vías fluviales interiores y áreas costeras. Sentinel-2 también puede proporcionar información para servicios de emergencia. Sentinel-2A fue lanzado el 23 de junio de 2015 y Sentinel-2B le siguió el 7 de marzo de 2017.

Sentinel-3 es una misión multiinstrumental para medir la topografía de la superficie del mar, la temperatura de la superficie del mar y la tierra, el color del océano y el color de la tierra, con alta precisión y confiabilidad. La misión respalda a los sistemas de pronóstico oceánico, así como al monitoreo ambiental y climático. Sentinel-3A fue lanzado el 16 de febrero de 2016 y Sentinel-3B se unió a su gemelo en órbita el 25 de abril de 2018.

Después de buscar un poco más, podemos aprender que los datos de Sentinel-1 llegan a unos pocos metros en cuanto a resolución espacial. Luego, la resolución más alta para datos visuales de Sentinel-2 es de 10 metros, mientras que Sentinel-3 opera a una escala mucho más grande, dependiendo del tipo de sensor, a escala de 100 km.

De acuerdo, ahora sabemos dónde obtener los datos satelitales, y también parece que hay una gran selección de fuentes (sensores) y resoluciones espaciales. Uno puede señalar que esto todavía es la punta del iceberg, como también se describe en esta lista de fuentes de datos satelitales. Y luego, ¿para qué usamos estos diferentes tipos de datos satelitales? Bueno, para empezar, aquí hay una selección de más de 50 casos de uso.

Como regla general, yo diría que el caso de uso, los detalles del problema y las características geoespaciales y del terreno del área objetivo son factores importantes para determinar la fuente de datos adecuada para ti. Sin embargo, en la práctica diaria, en mi experiencia, estos son los principales factores determinantes:

• precio (preferiblemente para explorar, gratis, lo que se aplica a Sentinel)
• tiene una resolución espacial de unos pocos metros, e incluso puede capturar estructuras urbanas más pequeñas
• con al menos algunas bandas, como visible e infrarrojo cercano
• frecuencia temporal

Estos aspectos hacen que Sentinel-2 sea probablemente la fuente de datos satelitales más utilizada en la comunidad de datos geoespaciales. Basándome en estos componentes, en este artículo te mostraré cómo obtener datos de Sentinel y qué puedes esperar al descargarlos. También exploraré a fondo las diferentes posibilidades y la evolución temporal de los registros de imágenes y la información almacenada en ellos.

En este artículo se utilizaron datos de Copernicus Sentinel 2023 para datos de Sentinel, ya que la legislación de la UE otorga acceso gratuito a los Datos y Servicios de Copernicus Sentinel.

1. Adquisición de datos

Primero, configuraré la conexión API siguiendo la documentación oficial y el código de muestra. Además, necesitaré un área objetivo desde la que descargar imágenes. Para facilitar la depuración, elegí mi ciudad natal, Budapest. Usaré OSMNx para descargar sus límites administrativos.

import osmnx as ox # versión: 1.0.1
import matplotlib.pyplot as plt # versión: 3.7.1
ciudad = 'Budapest'
admin = ox.geocode_to_gdf(ciudad)
admin.plot()

Ahora toca el turno de la API de Sentinel:

from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt # versión 0.14
# para obtener una cuenta, regístrese aquí: https://apihub.copernicus.eu/apihubuser = <agregue su nombre de usuario
password = <agregue su contraseña>
api = SentinelAPI(user, password, 'https://apihub.copernicus.eu/apihub')

Para hacer una consulta, es mejor tener un polígono suave que especifique la ubicación. Para este propósito, creé la envolvente convexa del área administrativa de Budapest:

# para simplificar la consulta, extraigo la envolvente convexa del polígono de entrada
admin_polygon = admin.convex_hull.geometry.to_list()[0]
admin_polygon

Salida:

Busque imágenes de satélite en nuestra ubicación específica, dentro de un marco de tiempo determinado y una plataforma de nuestra elección. Esta última será Sentinel-A. Además, también podemos filtrar según la cobertura de nubes, lo que significa que si una imagen está demasiado nublada, la descartamos de inmediato.

# aquí podemos especificar la ubicación (basada en un polígono)
# el marco de tiempo
# la sonda espacial
# y el nivel de cobertura de nubes aceptada
products = api.query(admin_polygon,
                     date=('20150623', '20231006'),
                     platformname='Sentinel-2',
                     cloudcoverpercentage=(0, 100))
len(products)

Como muestra la salida de estas celdas, siguiendo la documentación de Sentinel, resulta que entre el 23-06-2015 (comienzo de la misión) y el 06-10-2023 (cuando escribí este artículo), hubo un total de 3876 imágenes satelitales registradas que se superponen con los límites administrativos de Budapest. Establecí el porcentaje de cobertura de nubes en 100, lo que significa que no se aplicó ningún filtro en función de la cobertura de nubes. Por lo tanto, deberíamos tener todos los identificadores de imagen de los últimos ocho años.

También señalo aquí que la lista resultante de productos contiene el identificador y metadatos de todas las imágenes de satélite, sin embargo, no las propias imágenes. Además, si repitiera lo mismo con Sentinel-3, obtendría casi 20 mil registros de imágenes, aunque a una resolución mucho menor.

2. Explorar los metadatos

¡Transformemos la lista de productos en un DataFrame de Pandas y comencemos a analizarlo!

import pandas as pd # versión: 1.4.2
products_gdf = api.to_geodataframe(products)
products_gdf = products_gdf.sort_values(['beginposition'], ascending=[True])
print(products_gdf.keys())
print(len(products_gdf.keys()))
products_gdf.head(3)

Resultado de este bloque:

Después de contar el número de claves en la tabla, indexada por los identificadores de imágenes satelitales, uno puede tener una idea de lo rico que es este dato con las 41 columnas de características que contiene.

Aunque hay mucha información técnica en estos campos, me gustaría echar un vistazo más de cerca a solo algunas características. Por un lado, las dimensiones espaciales y temporales se codifican en la fecha de generación y la posición de inicio (como información de fecha y hora) y la geometría (como un polígono, tipo de dato GIS). Por otro lado, hay un par de métricas interesantes que caracterizan el tipo de cobertura del suelo basado en las imágenes: porcentaje de cobertura de nubes que ya hemos visto en las consultas, porcentaje de vegetación, porcentaje de agua y porcentaje de hielo/nieve. Estos índices ambientales se derivan de las diferentes propiedades espectrales de los diferentes materiales. Nota: estos valores son puntajes agregados que capturan el promedio general de la superficie completa. Más sobre este tema aquí.

3. La dimensión espacial

Como tenemos una dimensión de geometría, ¡echemos un vistazo a cómo se ve esto en un mapa! Lo haré visualizando un conjunto de mosaicos aleatorios, que, después de algunas ejecuciones, son totalmente representativos. Para la visualización, he utilizado Folium con el mapa base CartoDB Dark_Matter.

import foliumimport geopandas as gpdx, y = admin_polygon.centroid.xym = folium.Map(location=[y[0], x[0]], zoom_start=8, tiles='CartoDB Dark_Matter')# visualizar un conjunto de mosaicos aleatoriospolygon_style = { 'fillColor': '#39FF14', 'color': 'black',  'weight': 3, 'opacity': 0}geojson_data = products_gdf[['geometry']].sample(10).to_json()folium.GeoJson(    geojson_data,    style_function=lambda feature: polygon_style).add_to(m)# agregar los límites administrativosencimaadmin_style = {'fillColor': '#00FFFF',  'color': 'black','weight': 3, 'opacity': 100.0  }admin_geojson_data = admin[['geometry']].to_json()folium.GeoJson(    admin_geojson_data,    style_function=lambda feature: admin_style).add_to(m)# mostrar el mapam

La salida de este bloque de código:

Al observar esta visualización, queda claro que algunas secciones de mosaicos se repiten. También es evidente que, por alguna razón, algunos de estos mosaicos dividen los límites administrativos de la ciudad en dos. Esto puede llevar a una situación indeseada en la que desees analizar datos que cubren completamente tu área objetivo deseada, solo para darte cuenta de que se dividió por la mitad. Una posible solución alternativa es filtrar los mosaicos que no se superponen completamente con el área administrativa deseada:

def calcular_area_superpuesta(mosaico, admin):    return mosaico.intersection(admin_polygon).area / admin_polygon.areaproducts_gdf['fraccion_area_superpuesta'] = products_gdf.geometry.apply(lambda x: calcular_area_superpuesta(x, admin_polygon))products_gdf_f = products_gdf[products_gdf.fraccion_area_superpuesta==1]print(len(products_gdf))print(len(products_gdf_f))products_gdf_f.head(3)

El resultado de esta celda:

Al aplicar este filtro, se eliminaron aproximadamente la mitad de los mosaicos. Ahora echemos un vistazo al mapa y veamos cómo se superponen todos estos mosaicos con los límites de la ciudad sin que ninguno divida la ciudad en dos:

import foliumimport geopandas as gpdx, y = admin_polygon.centroid.xym = folium.Map(location=[y[0], x[0]], zoom_start=8, tiles='CartoDB Dark_Matter')# visualizar un conjunto de mosaicos aleatoriospolygon_style = { 'fillColor': '#39FF14', 'color': 'black',  'weight': 3, 'opacity': 0}geojson_data = products_gdf_f[['geometry']].sample(10).to_json()folium.GeoJson(    geojson_data,    style_function=lambda feature: polygon_style).add_to(m)# agregar los límites administrativosencimaadmin_style = {'fillColor': '#00FFFF',  'color': 'black','weight': 3, 'opacity': 100.0  }admin_geojson_data = admin[['geometry']].to_json()folium.GeoJson(    admin_geojson_data,    style_function=lambda feature: admin_style).add_to(m)# mostrar el mapam

4. La dimensión temporal del conjunto de datos

Primero, veamos el número de imágenes que tenemos por día, semana y mes, que cubren Budapest. Para medir el tiempo, me basaré en el campo beginposition.

# Suponiendo que 'beginposition' es una columna de tiempo en su GeoDataFrame# Puede convertirla en un índice de fecha y hora.products_gdf_f_cntr = products_gdf_f.copy()products_gdf_f_cntr['beginposition'] = pd.to_datetime(products_gdf_f_cntr['beginposition'])products_gdf_f_cntr.set_index('beginposition', inplace=True)# Reorganizar los datos para contar las filas por día, semana y mesdaily_counts = products_gdf_f_cntr.resample('D').count()weekly_counts = products_gdf_f_cntr.resample('W').count()monthly_counts = products_gdf_f_cntr.resample('M').count()fig, ax = plt.subplots(1, 3, figsize=(15, 5))for idx, (count_name, count_val) in enumerate([('Recuento diario', daily_counts), ('Recuento semanal', weekly_counts), ('Recuento mensual', monthly_counts), ]):     ax[idx].plot(count_val.index[0:250], count_val['geometry'].to_list()[0:250])    ax[idx].set_xlabel('Fecha')    ax[idx].set_ylabel('Cantidad')    ax[idx].set_title(count_name)plt.tight_layout()plt.suptitle('Número de imágenes de satélite tomadas en varios marcos de tiempo', fontsize = 20, y = 1.15)plt.show()

Estas figuras muestran los primeros 250 días, las primeras 250 semanas y los primeros 250 meses (todo el período) de la sonda Sentinel-2. La primera figura muestra que se toma una instantánea cada dos días. El segundo gráfico muestra, tomando el promedio semanal del anterior, que durante los dos primeros años, los satélites capturaron Budapest una o dos veces por semana, y luego, de 2017 a 2018, aumentó a 5-6 imágenes por semana. El último gráfico, que muestra todo el período de tiempo, muestra las mismas tendencias y cómo, después de 3 años de trabajo, estas 25 imágenes por mes se convirtieron en el nivel estándar.

5. La dimensión temporal de las variables de cobertura terrestre

Ahora, veamos la evolución temporal de la vegetación, el agua, la nieve/hielo y la cobertura de nubes. Como muestra la imagen anterior, los primeros años pueden mostrar diferentes resultados, probablemente ruido, así que seamos cautelosos. Aquí, no eliminaría esos años de datos, ya que tenemos un total de ocho años, y eliminar 3 de ellos podría significar perder demasiada información. Primero, veamos los valores brutos a lo largo del tiempo con una agregación semanal:

import pandas as pdimport matplotlib.pyplot as plt# Suponiendo que 'beginposition' es una columna de tiempo en su GeoDataFrame# Puede convertirla en un índice de fecha y horaproducts_gdf_f_cntr = products_gdf_f.copy()products_gdf_f_cntr['beginposition'] = pd.to_datetime(products_gdf_f_cntr['beginposition'])products_gdf_f_cntr.set_index('beginposition', inplace=True)# Reorganizar los datos para calcular promedios semanalesweekly_averages = products_gdf_f_cntr[['vegetationpercentage', 'waterpercentage', 'snowicepercentage', 'cloudcoverpercentage']].resample('W').mean()# Crear una figura de múltiples subgráficos con cuatro subgráficos(fig, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1, figsize=(10, 15)))# Graficar 'vegetationpercentage' con una línea verdeax1.plot(weekly_averages.index, weekly_averages['vegetationpercentage'], color='green', label='Promedio semanal del porcentaje de vegetación')ax1.set_xlabel('Fecha')ax1.set_ylabel('Porcentaje')ax1.set_title('Promedio semanal del porcentaje de vegetación')ax1.legend()# Graficar 'waterpercentage' con una línea azulax2.plot(weekly_averages.index, weekly_averages['waterpercentage'], color='blue', label='Promedio semanal del porcentaje de agua')ax2.set_xlabel('Fecha')ax2.set_ylabel('Porcentaje')ax2.set_title('Promedio semanal del porcentaje de agua')ax2.legend()# Graficar 'snowicepercentage' con una línea cianax3.plot(weekly_averages.index, weekly_averages['snowicepercentage'], color='cyan', label='Promedio semanal del porcentaje de nieve/hielo')ax3.set_xlabel('Fecha')ax3.set_ylabel('Porcentaje')ax3.set_title('Promedio semanal del porcentaje de nieve/hielo')ax3.legend()# Graficar 'cloudcoverpercentage' con una línea grisax4.plot(weekly_averages.index, weekly_averages['cloudcoverpercentage'], color='gray', label='Promedio semanal del porcentaje de cobertura de nubes')ax4.set_xlabel('Fecha')ax4.set_ylabel('Porcentaje')ax4.set_title('Promedio semanal del porcentaje de cobertura de nubes')ax4.legend()plt.tight_layout()plt.show()

Y lo mismo con la agregación mensual:

Estas series temporales nos están mostrando un par de cosas interesantes:

• El porcentaje de vegetación muestra claramente cómo cambia la estacionalidad, cómo en cada primavera todo se vuelve verde y luego esta verdor se desvanece durante el otoño, disminuyendo del 50 al 60% casi a cero.
• En comparación con eso, el porcentaje de agua fluctúa alrededor del 0.8% durante todo el año y todo el período de observación. Esto se debe a que la cantidad de agua superficial es muy pequeña en el área estudiada. Aún así, las caídas parecen ser más frecuentes durante el invierno, lo que implica que alguna de las entidades de agua dulce se congela.
• En cuanto a la nieve, los picos más prominentes, alrededor del 4 al 8%, ocurren durante el invierno. Sin embargo, ya que he estado presente durante la mayoría de estos inviernos, por experiencia personal, puedo decir que no tuvimos mucha nieve. En este sentido, los valores que miden solo el 1 al 2%, especialmente fuera del invierno, pueden generar algo de ruido o incluso una clasificación errónea de las nubes.
• En cuanto a las nubes, vemos que generalmente van de la mano con la vegetación, siguiendo patrones estacionales.

Algunas de estas observaciones también son visibles en la correlación de estas medidas:

products_gdf_f_cntr[['vegetationpercentage', 'waterpercentage', 'snowicepercentage', 'cloudcoverpercentage']].corr()

6. Descargar imágenes de satélite

Primero, haz una consulta tanto para Sentinel-2 como para Sentinel 3, eligiendo la misma semana de agosto de este año y limitando la cobertura de nubes tanto como sea posible. Observa la cantidad de imágenes disponibles:

# ids de productos de la consulta products_sent = api.query(admin_polygon, date=('20230806', '20230813'), platformname='Sentinel-2', cloudcoverpercentage=(0, 1))products_sent = api.to_geodataframe(products_sent)f, ax = plt.subplots(1,1,figsize=(6,4))admin.plot(ax=ax, color = 'none', edgecolor = 'k')ax.set_title('Sentinel-2, número de imágenes = ' + str(len(products_sent)))products_sent.plot(ax=ax, alpha = 0.3)# filtrar las imágenes que no se superponen completamente con Budapestproducts_sent['overlapping_area_fraction'] = products_sent.geometry.apply(lambda x: compute_overlapping_area(x, admin_polygon))products_sent = products_sent[products_sent.overlapping_area_fraction==1]f, ax = plt.subplots(1,1,figsize=(6,4))admin.plot(ax=ax, color = 'none', edgecolor = 'k')ax.set_title('Sentinel-2, número de imágenes = ' + str(len(products_sent)))products_sent.plot(ax=ax, alpha = 0.3)len(products_sent)

El resultado de este bloque de código:

Ahora descarga los datos de las imágenes basándote en

# descargar las primeras imágenes como imágenes satelitalesproduct_ids = products_sent.index.to_list()for prod in product_ids:    api.download(prod)

Nota — aquí es posible que obtenga esta notificación, en este caso simplemente espere algunas horas y vuelva a ejecutar el programa de descarga.

El producto a3c61497-d77d-48da-9a4d-394986d2fe1d no está en línea. Desencadenando la recuperación del archivo de archivo a largo plazo.

7. Abrir y visualizar la imagen descargada

Aquí encontrará una descripción detallada del formato de datos con visuals muy agradables sobre la estructura de carpetas. Una vez que abra el directorio de la imagen, podrá encontrar las diferentes bandas. El significado de cada banda, así como su resolución espacial, está bien resumido en este artículo, ya que la resolución espacial de las 13 bandas varía de 10 a 60 metros. Algunos puntos destacados:

• Los canales azul (B2), verde (B3), rojo (B4) e infrarrojo cercano o NIR (B8) tienen una resolución de 10 metros.
• Después, luego, su borde rojo de vegetación (B5), infrarrojo cercano NIR (B6, B7 y B8A) e infrarrojos de onda corta SWIR (B11 y B12) tienen una resolución de 10 metros.
• Finalmente, su aerosol costero (B1) y banda de cirrus SWIR (B10) tienen un tamaño de píxel de 60 metros.

Así es como se vería después de descomprimir la carpeta descargada:

# importar ospath_imagen = 'S2B_MSIL1C_20230810T094549_N0509_R079_T34TCT_20230810T124346.SAFE/GRANULE/L1C_T34TCT_A033567_20230810T095651/IMG_DATA'sorted(os.listdir(path_imagen))

Así es como se ve una mosaico completo al visualizar B4, la banda roja, usando rasterio:

import rasteriofrom rasterio.plot import showfile = 'T34TCT_20230810T094549_B04.jp2' with rasterio.open(path_imagen + '/' + file) as src:        imagen = src.read(1)  # Cambie el índice de la banda según sea necesario    plt.figure(figsize=(10, 10))    plt.imshow(imagen, cmap='Reds')  # Puede cambiar la escala de colores    plt.title(file)    plt.colorbar()    plt.show()

Salida:

Ahora enfoquémonos en Budapest y enmascaremos el mosaico completo de ráster con los límites administrativos de la ciudad, las bandas R, G y B por separado:

from rasterio import maskf, ax = plt.subplots(1,3,figsize=(15,5))for idx, (nombre_banda, num_banda, mapa_color) in enumerate([('Azul', 'B02', 'Blues'), ('Verde', 'B03', 'Greens'), ('Rojo', 'B04', 'Reds')]):       ruta_raster = path_imagen + '/T34TCT_20230810T094549_' + num_banda + '.jp2'    with rasterio.open(ruta_raster) as src:        poligonos = admin.copy().to_crs(src.crs)        geom = poligonos.geometry.iloc[0]        imagen_enmascarada, _ = mask.mask(src, [geom], crop=True)    ax[idx].imshow(imagen_enmascarada[0], cmap=mapa_color)    ax[idx].set_title('Sentinel 2 de Budapest - banda ' + nombre_banda)

Finalmente, vamos a ensamblar estos para obtener una imagen RGB de Budapest. Primero, ensamblamos la imagen completa en una imagen RGB, luego la leemos, realizamos una ecualización del histograma siguiendo las instrucciones oficiales, y luego llegamos al resultado final.

# Obtén las ubicaciones de las bandasbanda_azul = '/T34TCT_20230810T094549_B02.jp2'banda_verde = '/T34TCT_20230810T094549_B03.jp2'banda_roja = '/T34TCT_20230810T094549_B04.jp2'# Lee las bandas y crea la imagen completa RGBazul   = rasterio.open(ruta_imagen + '/' + banda_azul)verde   = rasterio.open(ruta_imagen + '/' + banda_verde)roja   = rasterio.open(ruta_imagen + '/' + banda_roja)# Exporta la imagen completa como un archivo tifmeta = roja.metameta.update({"count": 3})ruta_rgb_final = 'budapest_rgb.tif'with rasterio.open(ruta_rgb_final, 'w', **meta) as dest:    dest.write(azul.read(1),1)    dest.write(verde.read(1),2)    dest.write(roja.read(1),3)        # Lee y muestra la versión recortadafrom skimage import exposureimg = rasterio.open('budapest_rgb_cropped.tif')imagen = np.array([img.read(3), img.read(2), img.read(1)])imagen = imagen.transpose(1,2,0)# Realiza la ecualización del histogramap2, p98 = np.percentile(imagen, (2,98))imagen = exposure.rescale_intensity(imagen, in_range=(p2, p98)) / 100000f, ax = plt.subplots(1,1,figsize=(15,15))rasterio.plot.show(imagen.transpose(2,0,1), transform=img.transform, ax = ax)ax.axis('off')plt.savefig('budapest_rgb_cropped.png', dpi = 700, bbox_inches = 'tight')

Resultado:

Conclusiones

Para resumir rápidamente, veamos lo que sucedió en este artículo:

• Una visión general rápida de las plataformas de satélite Sentinel
• Un ejemplo de cómo consultar muchos identificadores de imagen, incluidos sus metadatos
• Cómo realizar análisis temporales basados ​​simplemente en la información agregada de las imágenes, almacenada y dirigida en los metadatos
• Cómo descargar, almacenar y visualizar imágenes individuales

Todos estos pasos apuntan a agregar el procesamiento y análisis de imágenes satelitales a tu conjunto diario de herramientas de ciencia de datos geoespaciales, lo que podría cubrir muchos casos de uso, desde la planificación urbana hasta el monitoreo ambiental y agrícola.

We will continue to update Zepes; if you have any questions or suggestions, please contact us!