Etiqueta: satellite imagery

Tras varios artículos sobre el uso de información satélite vamos a ver cómo unirlo (casi) todo en un ejemplo práctico e impactante. Y es que el tamaño de los últimos incendios ocurridos en España han llamado mucho mi atención y no era capaz de hacerme una idea de lo brutales que han sido (aunque nada que ver con otros sitios como Chile, Australia, o EEUU). Pero hagamos el ejercicio sin gastar demasiados GB de información geográfica.

Lo que quiero es mostrar la extensión del incendio ocurrido en Asturias en marzo, pero quiero también intentar mostrar el impacto retirando los árboles afectados por el fuego. ¡Vamos allá!

Descarga de datos

Usaré un Modelo Digital de Superficies (que incluye árboles y estructuras), una ortofoto tomada durante el incendio, y un Modelo Digital del Terreno (al que han eliminado los árboles y las estructuras) para reemplazarlo por las zonas afectadas por el incendio.

1. Modelos del terreno

Usaré los modelos del genial IGN, descargando los productos MDS5 y MDT5 de la zona.

http://centrodedescargas.cnig.es/CentroDescargas/index.jsp

2. Ortofotos

Para la imagen satélite, finalmente me decanto por Landsat ya que contaba con una imagen despejada durante los últimos días del incendio.

Usaré las imágenes tomadas el día 17 de febrero de 2023 (antes del incendio) y del 6 de abril de 2023 (ya en sus últimos días).

https://search.earthdata.nasa.gov

Procesar imágenes satélite

Usaremos el proceso i.group de GRASS en QGIS para agrupar las distintas bandas capturadas por el satélite en un único ráster RGB, tal y como vimos en este post:

https://theroamingworkshop.cloud/b/1725/procesar-imagenes-satelite-de-landsat-o-sentinel-2-en-qgis/

Tendremos que hacerlo para cada una de las regiones descargadas, cuatro en mi caso, que luego volveremos a unir usando el proceso Construir ráster virtual.

1. Imagen de color verdadero (TCI)

Combinamos las bandas 4, 3, 2.

2. Imagen de falso color

Combinamos las bandas 5, 4, 3.

3. Ajuste de la tonalidad

Para obtener un mejor resultado, puedes regular los valores mínimos y máximos que se consideran en cada banda que compone la imagen. Estos valores se encuentran en el Histograma de las propiedades de la capa.

Aquí te dejo los valores que yo he usado para obtener el resultado de arriba:

Banda	TCI min	TCI max	FC min	FC max
1 Rojo	-100	1500	-50	4000
2 Verde	0	1500	-100	2000
3 Azul	-10	1200	0	1200

Extensión del incendio

Como ves, la imagen en falso color nos muestra claramente la extensión del incendio. Con ella, vamos a generar un polígono que delimite el alcance del incendio.

Primero consultaremos los valores de la banda 1 (rojo) que ofrece mayor contraste para la zona del incendio. Más o menos están en el rango 300-1300.

Usando el proceso Reclasificar por tabla, asignaremos el valor 1 a las celdas dentro del rango, y el valor 0 al resto.

Vectorizamos el resultado con el proceso Poligonizar y, contrastando con la imagen satélite, seleccionamos aquellos polígonos que correspondan con el incendio.

Usaremos la herramienta Disolver para unir todos los polígonos en un elemento, y Suavizar para redondear ligeramente los contornos.

Ahora obtenemos su inverso. Extraemos la extensión de la capa Landsat y, posteriormente, hacemos la Diferencia con el polígono del incendio.

Procesar terreno

1. Combinar los datos del terreno

Lo primero que haremos es combinar los distintos archivos que conforman los modelos en un archivo único (un único fichero MDS y un único fichero MDT).

Usamos el proceso GDAL - Miscelánea ráster - Construir ráster virtual

2. Extraer datos del terreno

Extraemos los datos que nos interesan de cada modelo:

• Del MDS extraemos la superficie afectada por el incendio, de modo que quitaremos los árboles que hayan en él.
• Con el MDT hacemos lo inverso: dejamos el terreno (sin árboles) de la zona del incendio, para sustituir los huecos generados en el otro modelo.

Usaremos el proceso Cortar ráster por capa de máscara empleando las capas generadas en el apartado anterior.

Finalmente unimos ambas capas ráster, para que rellenen una a la otra, usando Construir ráster virtual.

Dale vida con Cesium JS

Ya deberíamos tener un modelo de superficies sin árboles en la zona del incendio, pero vamos a intentar verlo de forma interactiva.

Ya mostré un ejemplo parecido, usando un Modelo Digital del Terreno personalizado, así como una imagen satélite reciente, del volcán Tajogaite de La Palma:

https://theroamingworkshop.cloud/b/1319/cesiumjs-el-visor-gratuito-de-mapas-en-3d-para-tu-web/

En este caso volveré a usar Cesium JS para poder interactuar fácilmente con el mapa (sigue el post anterior para ver cómo subir tus ficheros personalizados al visor Cesium JS).

Para esta ocasión he creado una pantalla dividida (usando dos instancias de CesiumJS) para poder comparar el antes y el después del incendio. Aquí tienes una vista previa:

https://theroamingworkshop.cloud/demos/cesiumJSmirror/

Espero que te guste! Aquí tienes el código completo y el enlace a github para que puedas descargarlo. Y recuerda, comparte tus dudas o comentarios en twitter!

🐦 @RoamingWorkshop

<html lang="en">
<head>
<meta charset="utf-8">
<title>Cesium JS mirror v1.0</title>
<script src="https://cesium.com/downloads/cesiumjs/releases/1.96/Build/Cesium/Cesium.js"></script>
<link href="https://cesium.com/downloads/cesiumjs/releases/1.96/Build/Cesium/Widgets/widgets.css" rel="stylesheet">
</head>
<body style="margin:0;width:100vw;height:100vh;display:flex;flex-direction:row;font-family:Arial">
<div style="height:100%;width:50%;" id="cesiumContainer1">
	<span style="display:block;position:absolute;z-Index:1001;top:0;background-color: rgba(0, 0, 0, 0.5);color:darkorange;padding:13px">06/04/2023</span>
</div>
<div style="height:100%;width:50%;background-color:black;" id="cesiumContainer2">
	<span style="display:block;position:absolute;z-Index:1001;top:0;background-color: rgba(0, 0, 0, 0.5);color:darkorange;padding:13px">17/02/2023</span>
	<span style="display:block;position:absolute;z-Index:1001;bottom:10%;right:0;background-color: rgba(0, 0, 0, 0.5);color:white;padding:13px;font-size:14px;user-select:none;">
		<b><u>Cesium JS mirror v1.0</u></b><br>
		· Use the <b>left panel</b> to control the camera<br>
		· <b>Click+drag</b> to move the position<br>
		· <b>Control+drag</b> to rotate camera<br>
		· <b>Scroll</b> to zoom in/out<br>
		<span><a style="color:darkorange" href="https://theroamingworkshop.cloud" target="_blank">© The Roaming Workshop <span id="y"></span></a></span>
    </span>
</div>
<script>
	// INSERT ACCESS TOKEN
    // Your access token can be found at: https://cesium.com/ion/tokens.
    // Replace `your_access_token` with your Cesium ion access token.

    Cesium.Ion.defaultAccessToken = 'your_access_token';

	// Invoke LEFT view
    // Initialize the Cesium Viewer in the HTML element with the `cesiumContainer` ID.

    const viewerL = new Cesium.Viewer('cesiumContainer1', {
		terrainProvider: new Cesium.CesiumTerrainProvider({
			url: Cesium.IonResource.fromAssetId(1640615),//get your asset ID from "My Assets" menu
		}),	  
		baseLayerPicker: false,
		infoBox: false,
    });    

	// Add Landsat imagery
	const layerL = viewerL.imageryLayers.addImageryProvider(
	  new Cesium.IonImageryProvider({ assetId: 1640455 })//get your asset ID from "My Assets" menu
	);
	
	// Hide bottom widgets
	viewerL.timeline.container.style.visibility = "hidden";
	viewerL.animation.container.style.visibility = "hidden";

    // Fly the camera at the given longitude, latitude, and height.
    viewerL.camera.flyTo({
      destination : Cesium.Cartesian3.fromDegrees(-6.7200, 43.175, 6000),
      orientation : {
        heading : Cesium.Math.toRadians(15.0),
        pitch : Cesium.Math.toRadians(-20.0),
      }
    });
    
    // Invoke RIGHT view

    const viewerR = new Cesium.Viewer('cesiumContainer2', {
		terrainProvider: new Cesium.CesiumTerrainProvider({
			url: Cesium.IonResource.fromAssetId(1640502),//get your asset ID from "My Assets" menu
		}),	  
		baseLayerPicker: false,
		infoBox: false,
    });    

	// Add Landsat imagery
	const layerR = viewerR.imageryLayers.addImageryProvider(
	  new Cesium.IonImageryProvider({ assetId: 1640977 })
	);
	
	// Hide bottom widgets
	viewerR.timeline.container.style.visibility = "hidden";
	viewerR.animation.container.style.visibility = "hidden";

    // Fly the camera at the given longitude, latitude, and height.
    viewerR.camera.flyTo({
      destination : Cesium.Cartesian3.fromDegrees(-6.7200, 43.175, 6000),
      orientation : {
        heading : Cesium.Math.toRadians(15.0),
        pitch : Cesium.Math.toRadians(-20.0),
      }
    });
    
    // Invoke camera tracker
    //define a loop
    var camInterval=setInterval(function(){

	},200);
    clearInterval(camInterval);
    
    document.onmousedown=trackCam();
    document.ondragstart=trackCam();
    
    //define loop function (read properties from left camera and copy to right camera)
    function trackCam(){
		camInterval=setInterval(function(){
			viewerR.camera.setView({
				destination: Cesium.Cartesian3.fromElements(
					  viewerL.camera.position.x,
					  viewerL.camera.position.y,
					  viewerL.camera.position.z
					),
				orientation: {
					direction : new Cesium.Cartesian3(
						viewerL.camera.direction.x,
						viewerL.camera.direction.y,
						viewerL.camera.direction.z),
					up : new Cesium.Cartesian3(
						viewerL.camera.up.x,
						viewerL.camera.up.y,
						viewerL.camera.up.z)
				},
			});
		},50);
	};
	//stop loop listeners (release mouse or stop scroll)
	document.onmouseup=function(){
		clearInterval(camInterval);
	};
	document.ondragend=function(){
		clearInterval(camInterval);
	};
	
	//keep the copyright date updated
	var y=new Date(Date.now());
	document.getElementById("y").innerHTML=y.getFullYear();
  </script>
</div>
</html>