05 Sep 2022 Webmapping y Teledetección
En este post vamos a mostrar cómo visualizar imágenes de satélite y realizar algunas operaciones entre bandas utilizando la librería OpenLayers. Para ello utilizaremos el módulo WebGLTile, que añade soporte de renderizado WebGL a imágenes GeoTIFF. En cuanto a las imágenes, trabajaremos con las que ofrece el satélite Sentinel-2 del programa Copernicus.
Obtención de las imágenes
En otros posts de este blog UNIGIS Cómo obtener imágenes de satélite Sentinel ya hemos comentado cómo acceder a imágenes del programa Copernicus sin embargo, en esta ocasión, usaremos el cloud de Amazon Web Services en el que se encuentra todo el catálogo de imágenes disponibles para la descarga.
Para encontrar la ruta a una determinada imagen, se puede utilizar AWS Command Line Interface. Como parámetros, habrá que especificar la fecha y el identificador de celda de la imagen. Este último parámetro se puede encontrar muy rápidamente utilizando, por ejemplo, el portal EO Browser.
La instalación de AWS CLI se puede llevar a cabo tanto en sistemas operativos Windows, Linux como macOS. En la página web de AWS se detalla el proceso: instalación de AWS CLI.
Una vez instalado AWS CLI en el sistema y conociendo tanto la fecha de la imagen a buscar como su identificador de celda, podemos ejecutar la siguiente consulta a través de la terminal de comandos:
aws s3 ls s3://sentinel-cogs/sentinel-s2-l2a-cogs/31/T/DG/2022/7/ --no-sign-request
En esta ocasión, realizamos una consulta para obtener la ruta a las escenas de la zona de Catalunya, en julio de 2022. Este será el resultado:
La ruta para acceder a la banda 4, por ejemplo, de una de estas escenas, sería:
https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/31/T/DG/2021/2/S2A_31TDG_20210201_0_L2A/B04.tif
Mientras que, por ejemplo, para acceder a la composición en color natural, utilizaríamos:
https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/31/T/DG/2021/2/S2A_31TDG_20210201_0_L2A/TCI.tif
OpenLayers
Tal y como hemos comentado, utilizaremos la librería OpenLayers para visualizar las imágenes de forma interactiva en un entorno web.
En primer lugar, crearemos un proyecto de node e instalaremos la librería OpenLayers y el parcel-bundler.
$npm install ol $npm install parcel-bundler
A continuación crearemos un fichero index.html en el que añadiremos el código html de la aplicación, y otro fichero main.js en el que habrá el código javascript.
También configuraremos el documento package.json, con las funciones start y build.
{
"name": "post_ol_satellite",
"version": "1.0.0",
"description": "",
"main": "index.js",
"scripts": {
"test": "echo \"Error: no test specified\" && exit 1",
"start": "parcel index.html",
"build": "parcel build --public-url . index.html"
},
"author": "josep sitjar",
"license": "ISC",
"dependencies": {
"ol": "^6.14.1",
"parcel-bundler": "^1.12.5"
}
}El código html, con la estructura de la página web, puede ser el siguiente:
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<title>OpenStreetMap Reprojection</title>
<!-- Pointer events polyfill for old browsers, see https://caniuse.com/#feat=pointer -->
<script src="https://unpkg.com/elm-pep@1.0.6/dist/elm-pep.js"></script>
<!-- The lines below are only needed for old environments like Internet Explorer and Android 4.x -->
<script src="https://cdn.polyfill.io/v3/polyfill.min.js?features=fetch,requestAnimationFrame,Element.prototype.classList,TextDecoder"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/core-js/3.18.3/minified.js"></script>
<style>
.map {
width: 100%;
height:100vh;
}
</style>
</head>
<body>
<div id="map" class="map"></div>
<script src="main.js"></script>
</body>
</html>
Se trata de una página muy simple, dado que por ahora únicamente la queremos para mostrar las imágenes, sin ninguna otra funcionalidad adicional.
Y finalmente, en el documento main.js añadiremos el código javascript necesario para crear el mapa con OpenLayers y visualizar la imagen en color natural de una de las escenas Sentinel-2 que hemos identificado en el apartado anterior:
import 'ol/ol.css';
import Map from 'ol/Map';
import OSM from 'ol/source/OSM';
import TileLayer from 'ol/layer/WebGLTile';
import View from 'ol/View';
import GeoTIFF from 'ol/source/GeoTIFF';
// Almacenamos la fuente de datos en formato GeoTIFF
// Een este caso, una imagen Sentinel-2 en color natural
const source = new GeoTIFF({
sources: [
{
url: 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/31/T/DG/2021/2/S2A_31TDG_20210201_0_L2A/TCI.tif',
},
],
});
// Creamos la capa TileLayer a partir de la fuente de datos creada anteriormente
const layer = new TileLayer({
source: source,
});
const map = new Map({
layers: [
// Añadimos la capa (layer) al mapa
layer
],
target: 'map',
view: source.getView(),
});El resultado será el siguiente:
Composición en falso color
A partir del proyecto de OpenLayers con el que trabajamos, resulta relativamente fácil visualizar la misma escena en falso color, utilizando por ejemplo las bandas correspondientes al NIR (banda 8), al Rojo (banda 4) y al Verde (banda 3).
Simplemente tendremos que añadir cada una de estas bandas en el objeto GeoTIFF:
const source = new GeoTIFF({
sources: [
{
// banda correspondiente al NIR
url: 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/31/T/DG/2021/2/S2A_31TDG_20210201_0_L2A/B08.tif',
max: 5000,
},
{
// banda correspondiente al Rojo
url: 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/31/T/DG/2021/2/S2A_31TDG_20210201_0_L2A/B04.tif',
max: 5000,
},
{
// banda correspondiente al Verde
url: 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/31/T/DG/2021/2/S2A_31TDG_20210201_0_L2A/B03.tif',
max: 5000,
},
],
});Este será el resultado:
Cálculo de índices de vegetación
Las capas WebGLTile acceptan la propiedad style, que puede ser utilizada para controlar el renderizado de la fuente de datos, así como ejecutar expresiones matemáticas. Gracias a ello, podremos generar y representar al vuelo una imagen con los valores relativos al índice NDVI.
Será necesario modificar las capas que forman parte de la fuente de datos, añadiendo aquellas necesarias para el cálculo del NDVI (Roja y NIR):
const source = new GeoTIFF({
sources: [
{
url: 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/31/T/DG/2021/2/S2A_31TDG_20210201_0_L2A/B04.tif',
max: 10000,
},
{
url: 'https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/31/T/DG/2021/2/S2A_31TDG_20210201_0_L2A/B08.tif',
max: 10000,
},
],
});Y a continuación, añadir la propiedad style de la capa WebGLTile, quedando del siguiente modo:
const layer = new TileLayer({
// indicamos la fuente de datos
source: source,
// añadimos la propiedad 'style'
style: {
color: [
// aplicamos una interpolación
'interpolate',
['linear'],
// expresión para calcular el NDVI a partir de las bandas de la fuente de datos
['/', ['-', ['band', 2], ['band', 1]], ['+', ['band', 2], ['band', 1]]],
// asignación de colores a cada valor de la capa NDVI
-0.2, // a los valores ndvi <= -0.2 se les asignará el color RGB 191,191,191
[191, 191, 191],
0, // a los valores ndvi entre -0.2 y 0, se les asignará un color resultado de la interpolación entre el color anterior y posterior
[255, 255, 224],
0.2,
[145, 191, 82],
0.4,
[79, 138, 46],
0.6,
[15, 84, 10],
],
},
});En la documentación de OpenLayers relativa a la definición de color de la propiedad style, se encuentra toda la documentación relativa al uso de expresiones, asignación de color y operadores de transformación: color.
Conclusión
De una forma relativamente rápida y fácil, podemos visualizar imágenes de satélite a través de un mapa web desarrollado con la librería OpenLayers. Asimismo, agregando diferentes capas a una fuente de datos GeoTIFF, y utilizando las propiedades style que ofrece la clase WebGLTile, conseguimos crear diferentes combinaciones de bandas y llevar a cabo operaciones entre ellas.
Se trata, en definitiva, de una metodología a tener en consideración para visualizar imágenes de satélite en entornos web.
