Autor
Una guía paso a paso sobre cómo recoger y preprocesar imágenes de satélite data para algoritmos de aprendizaje automático
Este artículo es la primera parte de una serie de 2 partes que explora los usos de los algoritmos de aprendizaje automático en las imágenes de satélite. Aquí nos centraremos en los pasos de recogida y preprocesamiento de data, mientras que la segunda parte mostrará el uso de diversos algoritmos de aprendizaje automático.
TL;DR
Este artículo lo hará:
Utilizaremos como ejemplo la detección y clasificación de campos agrícolas. Para los pasos de preprocesamiento se utilizan Python y jupyter notebook. El cuaderno jupyter para reproducir los pasos está disponible en github.
Este artículo asume fundamentos básicos en ciencia data y python.
La visión por ordenador es una fascinante rama del aprendizaje automático que se ha desarrollado rápidamente en los últimos años. Con el acceso público a muchas bibliotecas de algoritmos off the shelf (Tensorflow, OpenCv, Pytorch, Fastai ...) así como datasets de código abierto (CIFAR -10, Imagenet, IMDB ...) es muy fácil para un científico data adquirir experiencia práctica en este tema.
Pero ¿sabía que también puede acceder a imágenes de satélite de código abierto? Estos tipos de imágenes pueden utilizarse en multitud de casos de uso (agricultura, logística, energía...). Sin embargo, suponen retos adicionales para un científico data debido a su tamaño y complejidad. He escrito este artículo para compartir algunos de mis principales aprendizajes trabajando con ellas.
Paso 1 - Seleccione la fuente data abierta por satélite
En primer lugar, debe elegir qué satélite desea utilizar. Las características clave a tener en cuenta aquí son :
Estos son los satélites cuyos data están disponibles de forma gratuita :
Referencia de los sistemas de satélites
Dado que nuestro objetivo es detectar e identificar cultivos agrícolas, el espectro visible (que se utiliza para proporcionar imágenes en color estándar) debería ser bastante útil. Además, existen correlaciones entre el nivel de nitrógeno de los cultivos y las bandas espectrales NIR (infrarrojo cercano). En cuanto a las opciones gratuitas, Sentinel 2 es actualmente la mejor en términos de resolución con 10 m para esas bandas espectrales, por lo que vamos a utilizar este satélite.
Por supuesto, la resolución de 10 m es muy inferior en comparación con el estado del arte en imágenes por satélite, con satélites comerciales que alcanzan una resolución de 30 cm. El acceso a los satélites comerciales data puede desbloquear muchos casos de uso, pero a pesar de la creciente competencia con la entrada de nuevos actores en este mercado, siguen teniendo un precio elevado y, por tanto, son difíciles de integrar en un caso de uso a escala.
Paso 2 - Recoger las imágenes de satélite data de Sentinel 2
Ahora que hemos seleccionado la fuente de nuestro satélite, el nuevo paso es descargar las imágenes. Sentinel 2 forma parte del Programa Copérnico de la Agencia Espacial Europea. Podemos acceder a su data utilizando la función API python de SentinelSat.
Primero tenemos que crear una cuenta en Centro de acceso abierto Copernicus. Como alternativa, tenga en cuenta que este sitio web también ofrece una forma de consultar imágenes con una interfaz gráfica de usuario.
Una forma de especificar las zonas que desea explorar es proporcionar a la API un archivo geojson, que es un json que contiene las coordenadas GPS de posición (latitud y longitud) de una zona. Los geojson pueden crearse rápidamente en página web.
De este modo, podrá consultar todas las imágenes de satélite que se crucen con la zona proporcionada.
Aquí se enumeran todas las imágenes disponibles en a para una zona y un intervalo de tiempo determinados :
api = SentinelAPI(
credenciales[“nombre de usuario”],
credenciales[“contraseña”],
“https://scihub.copernicus.eu/dhus”
)
shape = geojson_to_wkt(read_geojson(geojson_path))
imágenes = api.query(
forma,
date=(date(2020, 5, 1), date(2020, 5, 10)),
platformname=”Sentinel-2″,
processinglevel = “Nivel-2A”,
1TP37Porcentaje de cobertura=(0, 30)
)
imágenes_df = api.to_dataframe(imágenes)
Observe el uso de argumentos:
imágenes_df
Imágenes_df es un pandas dataframe que contiene todas las imágenes que coinciden con nuestra consulta. Seleccionamos una con una cobertura cloud baja y la descargamos utilizando la api :
api.download(uuid)
Paso 3 - Comprender qué bandas espectrales utilizar y gestionar la cuantización
Una vez descargado el archivo zip podemos ver que en realidad contiene varias imágenes.
Esto se debe a que, como dijimos en el paso 1, Sentinel 2 capta muchas bandas espectrales y no sólo el rojo, el verde y el azul, y cada banda se capta en una sola imagen.
Además sus instrumentos no tienen todos la misma resolución espacial, por lo que hay una subcarpeta para cada resolución diferente de 10, 20 y 60 metros.
Si comprobamos el GRANULE//IMG_DATA carpeta vemos que hay una subcarpeta para cada resolución diferente: 10, 20 y 60 metros como las bandas espectrales. Vamos a utilizar las que tienen una resolución máxima de 10 metros ya que proporcionan las bandas 2, 3, 4 y 8 que son azul, verde, rojo y NIR. Las bandas 7, 8, 9 que son el borde rojo de la vegetación (bandas espectrales entre el rojo y el NIR que marcan una transición en el valor de reflectancia para la vegetación) también podrían ser interesantes para detectar cultivos pero como son de menor resolución vamos a dejarlas de lado ahora
Aquí cargamos cada banda como una matriz numpy 2D:
def obtener_banda(carpeta_imagen, banda, resolución=10):
subcarpeta = [f for f in os.listdir(carpeta_imagen + "/GRANULE") if f[0] == "L"][0]
ruta_carpeta_imagen = f"/GRANULE//IMG_DATA/Rm"
archivos_de_imagen = [im para im en os.listdir(ruta_carpeta_de_imagen) si im[-4:] == ".jp2"].
archivo_seleccionado = [im for im in archivos_de_imagen if im.split("_")[2] == banda][0]
con rasterio.open(f"/") como infile:
img = infile.read(1)
devolver img
band_dict =
para banda en ["B02", "B03", "B04", "B08"]:
diccionario_banda[banda] = get_banda(carpeta_imagen, banda, 10)
Una prueba rápida para asegurarse de que todo funciona es intentar recrear la imagen en color verdadero (es decir, RGB) y visualizarla utilizando opencv y matplotlib:
img = cv2.merge((diccionario_banda[“B04”], diccionario_banda[“B03”], diccionario_banda[“B02”]))
plt.imshow(img)
Imagen de satélite - Problema de intensidad de color (Copernicus Sentinel data 2020)
Sin embargo, existe un problema con la intensidad del color. Esto se debe a que la cuantización (número de valores posibles) difiere del estándar de matplotlib, que asume un valor int 0-255 o un valor float 0-1. Las imágenes estándar suelen utilizar una cuantización de 8 bits (256 valores posibles), pero las imágenes de Sentinel 2 utilizan una cuantización de 12 bits y un reprocesado convierte los valores a un entero de 16 bits (65536 valores).
Si queremos escalar un entero de 16 bits a un entero de 8 bits, en teoría deberíamos dividir por 256, pero en realidad esto da como resultado una imagen muy oscura, ya que no se utiliza toda la gama de valores posibles. El valor máximo de nuestra imagen es en realidad 16752 de 65536 y pocos píxeles alcanzan realmente valores superiores a 4000 por lo que dividir por 8 da en realidad una imagen con un contraste decente
img_procesado = img / 8
img_processed = img_processed.astype(int)
plt.imshow(img_procesada)
Imagen de satélite - Intensidad de color derecha (Copernicus Sentinel data 2020)
Ahora tenemos un tensor de 3 dimensiones de forma (10980, 10980, 3) de enteros de 8 bits en forma de matriz numpy
Tenga en cuenta que :
Paso 4 - Dividir las imágenes en tamaños listos para el aprendizaje automático
Podríamos intentar aplicar un algoritmo a nuestra imagen actual, pero en la práctica eso no funcionaría con la mayoría de las técnicas debido a su tamaño. Incluso con una gran potencia de cálculo, la mayoría de los algoritmos (y especialmente el aprendizaje profundo) quedarían fuera de juego.
Así que debemos dividir la imagen en fragmentos. Una cuestión es cómo establecer el tamaño del fragmento.
Decidimos optar por dividir cada imagen en una cuadrícula de 45 * 45 (siendo 45 convenientemente un divisor de 10980)
Pasada esta fijación de parámetros, este paso es sencillo utilizando la manipulación de matrices numpy . Almacenamos nuestros valores en un dict utilizando tuplas (x, y) como claves.
frag_count = 45
tamaño_frag = int(img_processed.shape[0] / frag_count)
frag_dict =
para y, x en itertools.product(range(frag_count), range(frag_count)):
frag_dict[(x, y)] = img_processed[y*tamaño_frag: (y+1)*tamaño_frag,
x*tamaño_frag: (x+1)*tamaño_frag, :]
plt.imshow(frag_dict[(10, 10)])
Fragmento de la imagen del satélite (Copernicus Sentinel data 2020)
Paso 5 - Vincule su data a las imágenes de satélite utilizando la conversión de GPS a UTM
Por último, necesitamos vincular nuestras imágenes con el data que tenemos sobre los campos agrícolas para permitir un enfoque de aprendizaje automático supervisado. En nuestro caso disponemos de data sobre campos agrícolas en forma de una lista de coordenadas GPS (latitud, longitud) , similar al geojson que utilizamos en el paso 2 y queremos ser capaces de encontrar la posición de sus píxeles. Sin embargo, en aras de la claridad, empezaré mostrando una conversión de píxeles a coordenadas GPS y luego mostraré la inversa
El metadata recogido al solicitar la API SentinelSat proporciona las coordenadas GPS de las esquinas de la imagen (columna de la huella) y nosotros queremos las coordenadas de un píxel concreto. La latitud y la longitud, que son valores angulares, no evolucionan linealmente dentro de una imagen, por lo que no podemos hacer una simple relación utilizando la posición del píxel. Hay formas de resolver esto utilizando ecuaciones matemáticas pero ya existen soluciones estándar en python.
Una solución rápida es convertir la posición del píxel a UTM (Universal Transverse Mercator) en el que una posición está definida por una zona así como una posición (x, y) (en unidad metro) que sí evoluciona linealmente con la imagen. Podemos entonces utilizar una conversión de UTM a Latitud / Longitud proporcionada en la biblioteca utm. Para ello, primero necesitamos obtener la zona UTM y la posición de la esquina superior izquierda de la imagen, que pueden encontrarse en el metadata de la imagen en color verdadero.
#Obteniendo meatadata
transform = rasterio.open(tci_ruta_archivo, driver=’JP2OpenJPEG’).transform
número_zona = int(tci_ruta_archivo.split(“/”)[-1][1:3])
zona_letra = tci_ruta_archivo.split(“/”)[-1][0]
utm_x, utm_y = transform[2], transform[5]
# Conversión de la posición del píxel a utm
este = utm_x + pixel_columna * 10
norte = utm_y + fila_pixel * - 10
# Conversión de UTM a latitud y longitud
latitud, longitud = utm.to_latlon(este, norte, número_zona, letra_zona)
La conversión de la posición GPS en píxeles se realiza utilizando las fórmulas inversas. En ese caso tenemos que asegurarnos de que la zona obtenida coincide con nuestra imagen. Esto puede hacerse especificando la zona. Si nuestro objeto no está presente en la imagen, el resultado será un valor de píxel fuera de límites.
# Conversión de latitud y longitud a UTM
este, norte, número_zona, letra_zona = utm.from_latlon(
latitud, longitud, número_zona_fuerza=número_zona
)
# Conversión de UTM a columna y fila
pixe_column = round((este - utm_x) / 10)
fila_pixel = round((norte - utm_y) / -10)
ver en bruto
Conclusión
¡Ya estamos listos para utilizar las imágenes de satélite para el aprendizaje automático !
En el próximo artículo de esta serie de 2 partes veremos cómo podemos utilizar esas imágenes para detectar y clasificar superficies de campos utilizando el aprendizaje automático supervisado y no supervisado.
Gracias por leer y no dude en siga el blog técnico de Artefact ¡si desea recibir una notificación cuando se publique este próximo artículo !
