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:

Query de imágenes de satélite usando aws
Fig. 1: Realizar una consulta con aws

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:

Visualización de las imágenes de satélite en color natural
Fig. 2: Visualización de la imagen en color natural

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:

Visualización de las imágenes de satélite en falso color
Fig. 3: Visualización de la imagen en falso color

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.

Visualización del NDVI de las imágenes de satélite
Fig. 3: Resultado de calcular el índice NDVI

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.

Josep Sitjar
Geógrafo y máster en medio ambiente, análisis y gestión del territorio. Actualmente estoy cursando un grado superior en desarrollo de aplicaciones web. Trabajo en el Servicio de SIG y Teledetección (SIGTE) de la Universitat de Girona participando en numerosos proyectos técnicos vinculados al mundo de los SIG y al desarrollo de aplicaciones web map.


Suscríbete a nuestra newsletter