Omitir los comandos de cinta
Saltar al contenido principal
Inicio de sesión

Saber más

Administrar permisosAdministrar permisos
Caso 1: Imágenes Sentinel2 - Detección de cambios en zonas rurales (post 1/2)
01/04/2020 - Producto, Formación, Noticia

Primer post en el que compararemos imágenes NDVI (y CN de consulta) de junio de 2019 respecto al mismo mes de 2017; para conocer los cambios en la Demarcación Forestal de Aezkoa-Quinto Real.

Para tener una referencia geográfica, incorporaremos el Mapa Base de IDENA por servicio teselado (WMTS) dando de alta el servicio:

 

http://idena.navarra.es/ogc/wmts/1.0.0/WMTSCapabilities.xml

Así el proyecto queda fijado al SRC EPSG:25830, que es el oficial. QGIS hará una transformación al vuelo excelente con el resto de información, aunque tenga otro SRC.

Importaremos las imágenes a seleccionar desde:

https://filescartografia.navarra.es/3_ORTOFOTOGRAFIA/3_00_SATELITES/3_00_1_SENTINEL2/

Deben de ser de épocas del año similares (en el ejemplo de junio) y con una distancia temporal de 1 o más años. En caso contrario, los falsos positivos y negativos en la detección aumentan de forma importante. Luego lo veremos, pero los criterios de aquello a lo que se llama "cambio" también deberían ser ajustados a otras circunstancias.

Para este ejercicio, descargaremos e incorporaremos a QGIS:

Ahora descargamos e incorporamos las Demarcaciones Forestales desde IDENA, la capa FOREST_Pol_DemarGuarde.shp. Es un shapefile, también en EPSG:25830, pero no importa que este en ZIP, porque QGIS lo incorpora igual al proyecto.

A través del método gráfico (para quién sepa cuál es) o de la tabla de atributos, seleccionamos la Demarcación de Aezkoa_Quinto-Real y exportamos esa feature a un shapefile en local bajo el nombre "Demarcación_Aezkoa_Quinto-Real" y la añadimos al proyecto. Podemos eliminar del proyecto (y del disco) la capa de Demarcaciones completa.

De momento la queremos solo de marco de referencia (su bounding box o envolvente). Sobre la capa, con botón derecho, solicitaremos "Zoom a la capa". Incluso, modificando el lienzo del mapa para que se ajuste lo máximo posible, ya que tener imágenes más extensas de lo necesario penalizará el rendimiento y no aportará nada al análisis.

El truco es: Como la Demarcación es muy "cuadrada" comparada con nuestros monitores panorámicos, se puede hacer el Panel de capas más ancho, lo que ayuda a enmarcar mejor el área; luego ya volveremos a darle su tamaño habitual. El método formal es crear una máscara que ajuste a la feature con el algoritmo «Geometría mínima delimitadora», con el método Envolver y creando una capa temporal, pero se ajusta en exceso al área. Hablaremos después de los dos métodos.

Ahora se trata de "recortar" las 4 imágenes a la zona a la vista en el lienzo (método 1) o a la envolvente (método 2). Sobre cada imagen, con botón derecho del ratón, Exportar – Guardar como GeoTIFF (mantener nombre del archivo original con _R [recortado]), igual SRC. En el método 1, con "Extensión de la vista del mapa", en el 2, con "Calcular a partir de la capa": Geometría delimitadora. Se mantiene también la resolución original de 10 x 10m. Dejar el check para que añada el archivo al mapa. Así, con las 4 imágenes.

Truco: Para no tener que estar escribiendo el nombre del archivo, antes de Exportar, en Propiedades – Fuente, podemos copiar al portapapeles el nombre y así solo habrá que poner el añadido "_R".

Al terminar, se pueden eliminar las imágenes fuente del proyecto, y, sobre todo, del disco, que son voluminosas y siempre podríamos volver a descargarlas si hacen falta en el futuro ¡Están más seguras en filescartografia que en vuestros equipos!

Los que siguieron el método 1, ya pueden restaurar la interface habitual de QGIS. Los del método 2, eliminar el archivo temporal "Geometría delimitadora". En ambos casos, modificar la simbología de la Demarcación de relleno sencillo a línea sencilla, para hacer transparente el polígono. Un color llamativo que no esté entre los de las imágenes y un grosor de 1 mm estará bien.

Un aspecto así en el Panel de Capas es el adecuado:

 

Conectar /desconectar las imágenes del NDVI de ambos años ya dan algunas pistas de los cambios. Las grandes manchas rojas son agua (embalse de Eugi y de Itoiz) y suelo desnudo (las explotaciones de la mina de Magnesitas de Navarra). Pero, mejor si hacemos el análisis con las propias imágenes.

Recordar la leyenda NDVI:

 

Y la naturaleza de las imágenes: 3 canales (B1: rojo, B2: verde y B3: azul), con valores de un byte, comprendidos entre 0 y 255, aunque la B3 es despreciable. Luego trabajaremos con la intensidad de B1 y B2, que producen una imagen que transita entre rojo intenso: B1>245 y B2<10 y verde intenso: B1<25 y B2>240, más o menos.

Interesante, analizar lugares destacados con Identificar:  , la visibilidad de ambos NDVI activada y el método de consulta "de arriba abajo". Conocidos los valores, pasamos a:

LA ECUACIÓN, paso a paso:

Distinguiremos 3 situaciones y les asignaremos clases (valores discretos):

  1. Sin cambios relevantes
  2. Pérdida de vegetación: Pasa de banda 2 (verde) alto en 2017 a bajo en 2019 y a la inversa en banda 1 (roja) de bajo a alto.
  3. Revegetación: De B2 bajo a alto y de B1 alto a bajo.

Construimos una ecuación que ejecutaremos en Ráster – Calculadora ráster (copiar más adelante ya completa)

Valor por defecto (sin cambios): 0 +         En realidad no es necesario porque es el valor cuando no se cumplen las otras condiciones, pero tiene intención didáctica para los que pregunten ¿y dónde está el cero?

Pérdida de vegetación:

("NDVI_20190628_R@1" - "NDVI_20170618_R@1" > 25 AND "NDVI_20190628_R@2" - "NDVI_20170618_R@2" < 25) * 1

25 es el parámetro de "sensibilidad" del cambio. En el dominio de valores (0 a 255) es el 10%, pero cada cuál puede poner otro valor y comparar resultados.

La primera parte: "NDVI_20190628_R@1" - "NDVI_20170618_R@1" > 25 recoge un aumento de la Banda 1 (rojo).

AND: Y a la vez…

La segunda parte: "NDVI_20190628_R@2" - "NDVI_20170618_R@2" < 25 recoge un descenso equivalente de la Banda 2 (verde).

* 1: Si el resultado anterior es "verdadero" devolverá un 1. Como en el caso del "0" anterior no es necesario, pero así se comprende mejor la asignación de las clases que tendrá el resultado.

+: Imprescindible para encadenar con la siguiente parte.

Revegetación: realiza la operación inversa (es excluyente) de la anterior.

("NDVI_20190628_R@1" - "NDVI_20170618_R@1" < 25 AND "NDVI_20190628_R@2" - "NDVI_20170618_R@2" > 25) * 2

En este caso, si es cierta, multiplica el valor para asignarle la clase=2.

Así, la ecuación completa es:

0 +

("NDVI_20190628_R@1" - "NDVI_20170618_R@1" > 25 AND "NDVI_20190628_R@2" - "NDVI_20170618_R@2" < 25) * 1 +

("NDVI_20190628_R@1" - "NDVI_20170618_R@1" < 25 AND "NDVI_20190628_R@2" - "NDVI_20170618_R@2" > 25) * 2

Que puede copiarse "tal cual" en Ráster – Calculadora ráster:

 

Capa de salida: Cambios-20170618-20190628.tif

Extensión de la capa seleccionada (seleccionar cualquier NDVI, por ejemplo)

SRC: el del proyecto, EPSG:3345 (pero podría transformarse a otro, como 25830)

Check de Añadir al proyecto, activado.

Abajo a la izquierda, debe figurar "Expresión válida". En caso contrario, no se activa el botón [Aceptar].

 

RESULTADO: Cambios-20170618-20190628

Situar la capa en la parte superior del Panel de capas, bajo la Demarcación.

Desactivar la visualización de todos los demás ráster excepto el Mapa Base y esta nueva capa.

La simbología por defecto es "Gris monobanda", porque el resultado solo contiene una única banda (1 Gray - gris), con los valores 0 a 2.

Cambiar la simbología a "Pseudocolor monobanda", Mín=0, Máx=2, Interpolación Lineal o Exacto. Sirven muchas rampas de color, pero se recomienda la socorrida "Spectral", aunque modificaremos los colores. Modo "Intervalo igual" de 3 clases. Clic sobre [Clasificar]

A continuación, se cambian colores y asignan valores a las etiquetas como en la figura:

 

En la pestaña "Transparencia", Valor sin datos, check activado y valor adicional 'sin dato'=0

[Aceptar]

Los valores 0 (sin cambios) se han hecho transparentes porque interesan en la leyenda, pero no queremos que impidan la visibilidad del resto del mapa.

Activando/desactivando la visualización de esta capa, las NDVI (incluso las CN) se puede apreciar la capacidad de detección de los cambios con esa "sensibilidad" del 10% (el valor 25 de la ecuación).

Sigue en el post 2/2.

SITNA

Datos adjuntos