GIS-curriculum

Módulo 9 - Procesamiento y análisis de ráster

Autor: Codrina

Introducción pedagógica

Este módulo se centra en un tipo específico de modelo de datos geográficos: geodatos ráster.

Al final de este módulo, los alumnos tendrán la comprensión básica de los siguientes conceptos:

• modelo de datos ráster
• punto de la cuadrícula vs celda de la cuadrícula
• bandas de un conjunto de datos ráster
• álgebra de mapas
• las cuatro resoluciones de un ráster (espacial, temporal, espectral y radiométrica)
• remuestreo
• procesamiento por lotes
• detección de cambios

Herramientas y recursos necesarios

• Este módulo se ha preparado utilizando QGIS versión 3.16.1 - Hannover
• module9.gpkg
• High Resolution Settlement Layer
• SRTM Digital Elevation Model
• Global Land Cover Map 2015-2019
• The coordinate reference system used is the SLD 99 / Sri Lanka Grid 1999, EPSG 5235.

Requisitos previos

• Conocimientos básicos de funcionamiento de una computadora
• Una sólida comprensión de los módulos 0, 1, 2, 6 y 8 de este plan de estudios.

Introducción temática

Desglose de los conceptos

El modelo de datos ráster

Un ráster es una matriz regular de valores, como en la figura 9.1.

Figura 9.1 - Una matriz de valores

Los valores se pueden asignar a puntos de la cuadrícula, principalmente a los centroides, y en este caso, el ráster puede denominarse enrejado o grilla. La segunda opción es que los valores se asignen a la totalidad de la celda de la cuadrícula, denominada píxeles (ver figura 9.2). Para el primer caso, el ráster generalmente representa un campo continuo, como elevación, temperatura, precipitaciones, concentraciones químicas, etc. Para el segundo caso, cuando los valores se asignan a toda el área del píxel, el ráster generalmente representa una imagen: una imagen de satélite, un mapa escaneado, mapas vectoriales convertidos (ver Fase 2). Cabe señalar que este modelo de datos no es particularmente eficiente para redes y otros tipos de datos que dependen en gran medida de las líneas, como los límites de las propiedades.

Figura 9.2 - En el lado izquierdo, los valores se asignan a centroides.

A la derecha, los valores se asignan al área de la celda de la cuadrícula: el píxel.

Racionalidad del procesamiento ráster

Como en el caso del procesamiento de datos vectoriales, la razón fundamental detrás del procesamiento de datos ráster se basa en la misma capacidad de los sistemas de información geográfica para almacenar, procesar y representar información de los fenómenos del mundo real, solo que la forma en que esto se logra difiere. En lugar de tener distintos puntos, líneas y polígonos almacenados como colecciones de coordenadas x e y, tenemos una matriz de valores que cubre un área específica como una malla. Para tener una imagen más clara en mente, imagine los mapas de temperatura que se muestran en la televisión. La temperatura es un campo continuo, no hay lugares en la superficie de la Tierra sin temperatura, ya sea positiva, negativa o 0.

Se pueden realizar muchas operaciones en conjuntos de datos ráster, el concepto de geoprocesamiento detallado en el módulo 8 también se aplica aquí. El término utilizado para englobar las operaciones que se pueden realizar en rásteres es álgebra de mapas.

El álgebra de mapas representa un conjunto de operaciones primitivas [^ 1] en un SIG que permite que dos o más capas ráster de dimensiones similares produzcan una nueva capa ráster (mapa) usando varias operaciones como suma, resta, comparación, etc.

Hay cuatro categorías de operaciones que se pueden realizar en rásteres, como sigue:

• Operaciones aritméticas: suma, resta, multiplicación, división.
• Operaciones estadísticas: mínimo, máximo, promedio, mediana.
• Las operaciones relacionales proporcionan comparaciones entre celdas que utilizan funciones como mayor que, menor que o igual a.
• Operaciones trigonométricas: seno, coseno, tangente, arcoseno entre dos o más conjuntos de datos ráster.
• Las operaciones exponenciales y logarítmicas usan funciones de exponentes y logaritmos.

Para cada una de estas operaciones, existen algoritmos implementados en la mayoría de los entornos SIG que permiten al usuario aplicarlos en sus datos. A continuación, también implementaremos algunas de las operaciones más comunes para tener una idea de cómo trabajar y qué se debe esperar de este tipo de procesamiento de datos.

El concepto de dimensiones similares se refiere a las características de los conjuntos de datos ráster. Es decir, las operaciones detalladas anteriormente no se pueden realizar con resultados significativos en 2 conjuntos de datos ráster con diferentes resoluciones espaciales, temporales o espectrales. A continuación, presentaremos muy brevemente las cuatro resoluciones que son relevantes para las imágenes rasterizadas.

Recordando cuáles son los 4 tipos de resoluciones de una imagen satelital (raster con valores asignados al área de celda - píxeles):

1. La resolución espacial corresponde al tamaño elemental de la superficie del suelo medida, se expresa en unidades de longitud y representa la longitud de un lado del píxel (ver figura 9.3). Por ejemplo, como ha visto en el módulo 3, los datos de la capa de asentamiento de alta resolución tienen una resolución espacial de 30 metros, es decir, cada píxel del conjunto de datos estima el número de personas que viven dentro de un metro de 30.

2. La resolución temporal está asociada con imágenes aéreas (imágenes adquiridas por satélites, aviones, helicópteros, drones) y corresponde al período entre 2 imágenes consecutivas del mismo punto exacto de la Tierra, tomadas en las mismas condiciones (tanto como sea posible), como como el mismo avión, la misma altitud, etc. Por ejemplo, Landsat 8 tiene una resolución temporal de 16 días, es decir, cada punto de la Tierra es revisado por el satélite Landsat 8 cada 16 días [^ 2].

3. Resolución espectral: los sensores a bordo de satélites o aviones capturan la radiación electromagnética proveniente de todos los objetos de la Tierra (agua, asentamientos, bosques, carreteras, edificios, tierra desnuda, etc.) y los sensores están diseñados específicamente para capturarla en una banda espectral conocida determinada. (o longitud de onda). El ojo humano puede ver una parte muy pequeña del espectro electromagnético: la luz visible (roja, verde y azul), ¡pero los sensores pueden ‘ver’ mucho más! (ver figura 9.4)

Figura 9.4 - Espectro electromagnético (crédito de la foto NASA Science - https://science.nasa.gov/ems/01_intro)

1. La resolución radiométrica está determinada por el número de bits en los que se divide la radiación registrada. En datos de 8 bits, los números digitales (DN) pueden oscilar entre 0 y 255 para cada píxel (28 ¼ 256 números posibles en total). Claramente, más bits significa que el sensor puede detectar cambios más sutiles en la energía que captura, lo que conduce a una imagen ‘más clara’, una mayor precisión radiométrica del sensor, pero también requiere más espacio para almacenar los datos.

Bandas espectrales

Un conjunto de datos ráster contiene una o más capas llamadas bandas. Cada banda almacena otro conjunto de información sobre el área que cubre el ráster. Cada banda tiene exactamente la misma extensión y coordenadas, pero no necesariamente la misma resolución espacial. Además, además de los valores almacenados, hay otras propiedades clave contenidas, como: valores de celda máximo, mínimo y medio e histograma de valores de celda.

Un histograma es una representación aproximada de la distribución de datos numéricos (ver figura 9.5); de otras formas, es una forma en la que uno puede tener una idea de los datos disponibles.

Figura 9.5 - Ejemplo de un histograma, donde x es una capa ráster

¿Por qué es esto importante en nuestro contexto? Porque, como se mencionó al principio, un ráster es una matriz de valores numéricos continuos (recuerde el ejemplo de temperatura) y un histograma ayuda al usuario a comprender cómo se distribuyen los valores de sus rásteres. Cada barra agrupa los valores de celda que se encuentran en un rango específico, mientras más alta sea la barra, más celdas tienen valores en ese rango específico. En caso de que el ráster tenga más de una banda, se calculará el histograma para cada una. No hay espacios entre los rangos representados en un histograma, el histograma se usa solo para datos continuos.

Contenido principal

Fase 1: Comprensión de sus datos ráster.

En las últimas dos décadas, la cantidad de satélites que captan datos de la Tierra ha crecido exponencialmente. Además, una política de acceso a datos abiertos que ha sido adoptada por diferentes agencias espaciales, como la NASA para el programa Landsat o la Agencia Espacial Europea para el programa Copernicus, lo cual ha abierto la puerta a un flujo abrumador de datos de observación de la Tierra. Como consecuencia natural, el progreso científico de los algoritmos, las metodologías y el desarrollo de herramientas más poderosas para procesar datos ráster, especialmente imágenes satelitales, ha sido impresionante y extenso en campos como la agricultura, la silvicultura, el desarrollo urbano, las actividades humanitarias, el océano y el mar. monitoreo de agua, seguridad y muchos más. En las siguientes 3 fases, presentaremos algunas de las técnicas de procesamiento más comunes y qué resultados esperar de ellas.

Paso 1. Prepare su entorno de trabajo.

Comenzaremos agregando a su proyecto QGIS todos los conjuntos de datos ráster con los que trabajaremos, de la siguiente manera:

• Datos de capa de asentamiento de alta resolución
• Modelo de superficie digital global ALOS “ALOS World 3D - 30m (AW3D30)”¹
• Cobertura terrestre dinámica moderada: las 5 épocas disponibles 2015, 2016, 2017, 2018 y 2019.²

Como los datos de la capa de asentamiento de alta resolución se presentaron en un módulo anterior, también detallaremos información sobre los otros 2 conjuntos de datos que emplearemos en este módulo.

1. El modelo de superficie digital global de ALOS (AW3D30) es un conjunto de datos de modelo de superficie digital global (DSM) con una resolución horizontal de aproximadamente 30 metros (malla de 1 segundo de arco). Un DSM incluye la superficie del suelo, la vegetación y los objetos artificiales, como edificios, puentes, etc., a diferencia del modelo de elevación digital que considera estrictamente el terreno.

2. Dynamic Land Cover Map (Mapa dinámico de cobertura terrestre) es un producto de Copernicus Global Land Service (CGLS) derivado de la serie de tiempo PROBA-V 100 my varios otros conjuntos de datos de cobertura terrestre. El producto proporciona clases discretas de cobertura terrestre primaria, así como capas de campo continuas para todas las clases de cobertura terrestre básica que proporcionan estimaciones proporcionales de vegetación / cobertura terrestre para los tipos de cobertura terrestre. El producto tiene 3 bandas: discrete_classification, (clasificación discreta), forest_type (tipo de bosque) y urban_coverfraction (fraccion de cobertura urbana). Las 2 tablas siguientes presentan los valores para cada clase discreta:

Tabla 1 Tabla de valores de la banda discrete_classification (clasificación discreta)

Valor	Color	Descripción
0	282828	Desconocido. No hay suficientes datos satelitales disponibles o no hay datos.
20	FFBB22	Arbustos. Plantas leñosas perennes con tallos persistentes y leñosos y sin tallo principal definido que mida menos de 5 m de altura. El follaje de los arbustos puede ser perenne o caducifolio.
30	FFFF4C	Vegetación herbácea. Plantas sin tallo persistente o brotes por encima del suelo y sin estructura firme definida. La cobertura de árboles y arbustos es inferior al 10%.
40	F096FF	Vegetación / agricultura cultivada y gestionada. Tierras cubiertas con cultivos temporales seguidos de cosecha y un período de suelo desnudo (por ejemplo, sistemas de cultivo único y múltiple). Tenga en cuenta que los cultivos leñosos perennes se clasificarán como el tipo apropiado de cobertura forestal o arbustiva.
50	FA0000	Urbano / edificado. Terreno cubierto por edificios y otras estructuras artificiales.
60	B4B4B4	Descubierto / escasa vegetación. Tierras con suelo, arena o rocas expuestas y nunca tienen más del 10% de cobertura vegetal durante cualquier época del año.
70	F0F0F0	Hielo y nieve. Aterriza bajo cubierta de nieve o hielo durante todo el año.
80	0032C8	Cuerpos de agua permanentes. Lagos, embalses y ríos. Pueden ser masas de agua dulce o salada.
90	0096A0	Humedal herbáceo. Tierras con una mezcla permanente de agua y vegetación herbácea o leñosa. La vegetación puede estar presente en agua salada, salobre o dulce.
100	FAE6A0	Musgos y líquenes.
111	58481F	Bosque cerrado, hoja de aguja siempre verde. Copas de los árboles> 70%, casi todos los árboles con hojas de aguja permanecen verdes todo el año. El dosel nunca está exento de follaje verde.
112	9900	Bosque cerrado, hoja ancha siempre verde. Copas de los árboles> 70%, casi todos los árboles de hoja ancha permanecen verdes durante todo el año. El dosel nunca está exento de follaje verde.
113	70663E	Bosque cerrado, hoja de aguja decidua. La copa de los árboles> 70%, consiste en comunidades de árboles de hojas de aguja estacionales con un ciclo anual de períodos de hojas abiertas y sin hojas.
114	00CC00	Bosque cerrado, hoja ancha caducifolia. La copa de los árboles> 70%, consiste en comunidades de árboles de hoja ancha estacionales con un ciclo anual de periodos de hojas abiertas y sin hojas.
115	4E751F	Bosque cerrado, mixto.
116	7800	Bosque cerrado, que no coincide con ninguna de las otras definiciones.
121	666000	Bosque abierto, hoja de aguja siempre verde. La capa superior - árboles 15-70% y la segunda capa - mezcla de arbustos y pastizales, casi todos los árboles de hojas de aguja permanecen verdes todo el año. El dosel nunca está exento de follaje verde.
122	8DB400	Bosque abierto, hoja ancha siempre verde. La capa superior - árboles 15-70% y la segunda capa - mezcla de arbustos y pastizales, casi todos los árboles de hoja ancha permanecen verdes durante todo el año. El dosel nunca está exento de follaje verde.
123	8D7400	Bosque abierto, hoja de aguja caducifolia. La capa superior - árboles 15-70% y la segunda capa - mezcla de arbustos y pastizales, consiste en comunidades de árboles de hojas de aguja estacionales con un ciclo anual de períodos de hojas abiertas y sin hojas.
124	A0DC00	Bosque abierto, hoja ancha caducifolia. La capa superior - árboles 15-70% y la segunda capa - mezcla de arbustos y pastizales, consiste en comunidades de árboles latifoliados estacionales con un ciclo anual de periodos con y sin hojas.
125	929900	Bosque abierto, mixto.
126	648C00	Bosque abierto, no coincide con ninguna de las otras definiciones.
200	80	Océanos, mares. Pueden ser masas de agua dulce o salada.

Tabla 2: tabla de valores de la banda forest_type (tipo de bosque)

Valor	Color	Descripción
0	282828	Desconocido
1	666000	Hoja de aguja perenne
2	9900	Hoja ancha de hoja perenne
3	70663E	Hoja de aguja caducifolio
4	A0DC00	Hoja ancha caducifolia
5	929900	Mezcla de tipos de bosque

Para organizar mejor sus capas, agrúpelas por categoría, de la siguiente manera: para las 5 capas ráster de cobertura terrestre, cree un grupo denominado Cobertura terrestre (en la tabla de contenidos, haga clic en el pictograma Agregar grupo). Para el modelo de superficie digital, cree un grupo llamado JAXA DSM.

No olvide agregar también el límite del área de trabajo, provincia de Santa Fe, que es un conjunto de datos vectoriales.

Su ventana de mapa QGIS debería verse como en la figura 9.6, tal vez en colores ligeramente diferentes.

Figura 9.6 - conjuntos de datos ráster cargados

Paso 2. Comprenda lo que está viendo.

A continuación, utilizaremos una serie de herramientas que nos permitirán tener una idea de los datos con los que estamos trabajando.

Después de cargar todos los conjuntos de datos, verificaremos el sistema de referencia de coordenadas en el que se encuentran todos nuestros conjuntos de datos. Como sabemos por módulos anteriores, QGIS ofrece la posibilidad de reproyectar todos los conjuntos de datos cargados en el proyecto sobre la marcha, sin embargo, eso podría conducir a problemas en un geoprocesamiento. Así, aunque todas las capas estén correctamente superpuestas, como se puede ver por inspección visual, procederemos a reproyectarlas todas en el sistema de coordenadas oficial de nuestra región de interés, provincia de Santa Fe - EPSG: 5347.

Hay varias formas de obtener información de las capas cargadas en QGIS, algunas proporcionan al usuario más detalles que otras. Para obtener una descripción general rápida de los metadatos de un conjunto de datos, se debe seleccionar la capa, hacer doble clic y abrir Propiedades ‣ Información.

Para la capa HRSL_SantaFe_poblacion, la ventana de información se vería como en la figura 9.7.

Figura 9.7 - Extracción de metadatos básicos de una capa ráster

Con respecto a nuestra primera pregunta sobre qué CRS se está utilizando para los conjuntos de datos que hemos cargado, podemos observar que incluso si el HRSL está correctamente superpuesto, la proyección del conjunto de datos nativo es EPSG 4326 - WGS 84 - Geográfico, con unidades medidas en grados. También identificamos que esta capa ráster específica tiene solo una banda, sin embargo, el tamaño de píxel es difícil de leer ya que la medición está en grados y no en metros, lo que la haría más fácil de entender.

Por lo tanto, lo primero que debe hacer es reproyectar todos los conjuntos de datos con los que trabajaremos en el mismo sistema de coordenadas - EPSG 5347.

Comenzando con los conjuntos de datos HRSL, vamos a Ráster ‣ Proyecciones ‣ Combar (Reproyectar) (ver figura 9.8).

Figura 9.8 - Funcionalidad de reproyección en QGIS

Después de seleccionar la funcionalidad Combar, aparecerá una nueva ventana que le permitirá al usuario configurar los parámetros correctos (ver figura 9.9).

Figura 9.9 - Ventana de Combar (reproyectar) QGIS

If you selected the output to be [Save to temporary file] then there will be a raster layer named Reprojected in the Layers Panel. This is a memory layer and you can rename this layer to Reprojected_HRSL_Santa Fe_Population and save it to make it persistent.

Figure 9.9b - Reprojected HRSL

Notará que, a diferencia de cuando reproyectó conjuntos de datos vectoriales, hay un nuevo parámetro que puede configurar como “Método de muestreo para usar”.

Muestreo representa la interpolación de los valores de la celda para que transforme el ráster como lo indica el usuario. Hay varios métodos de remuestreo disponibles dentro de la funcionalidad Deformación, cada uno con su propio soporte matemático. Sin embargo, las explicaciones detalladas sobre cada uno no son el alcance de este ejercicio. La lectura adicional está disponible en las referencias.

En este caso particular, queremos reproyectar datos de población - valores numéricos. Y según el método de remuestreo seleccionado (vecino más cercano), la coordenada de cada píxel de salida se utilizará para calcular un nuevo valor a partir de los valores de píxel cercanos en la capa de entrada (ver figura 9.10).

Figura 9.10 - Método de remuestreo - vecino más cercano (crédito de la foto: documentación de ILWIS - (http://spatial-analyst.net/ILWIS/htm/ilwisapp/resample_functionality.htm)

Los píxeles de entrada están representados por líneas negras discontinuas, las coordenadas de los píxeles de entrada por puntos negros; los píxeles de salida están representados por líneas sólidas rojas, las coordenadas de los píxeles de salida por signos más rojos. Las flechas grises indican cómo se determinan los valores de salida. Se puede ver en la figura 9.10, algunos valores del mapa de entrada se pueden usar dos veces en el mapa de salida, mientras que otros valores de entrada pueden no usarse en absoluto. Por eso, aunque es uno de los métodos más rápidos para volver a muestrear, no es apropiado en nuestro caso, ya que estamos trabajando con datos numéricos: datos de población. Este método de remuestreo es adecuado para datos categóricos, como valores de cobertura terrestre.

Para reproyectar nuestro conjunto de datos ráster de población, utilizaremos el método de interpolación bilineal para volver a muestrear los valores de píxeles (ver figura 9.11).

Figura 9.11 Método de remuestreo - bilineal (crédito de la foto: ILWIS documentación - (http://spatial-analyst.net/ILWIS/htm/ilwisapp/resample_functionality.htm)

El método bilineal determina el nuevo valor de una celda basado en un promedio de distancia ponderado de los cuatro centros de celda de entrada más cercanos. Es útil para datos continuos y suavizará los datos.

Procedemos a verificar el CRS de los conjuntos de datos de cobertura terrestre que hemos cargado en nuestro proyecto QGIS. Al acceder a la información de Propiedades de capa ‣, podemos ver que las capas de cobertura terrestre se proyectan en EPSG: 4326 - WGS 84. Una solución sería utilizar la herramienta Reproyectar y configurar para cada capa individualmente. Sin embargo, una forma más rápida es utilizar la función de reproyección que se ejecuta como un procesamiento por lotes.

El procesamiento por lotes es la capacidad de ejecutar procesos repetitivos en datos, con una mínima interacción del usuario. La mayoría de las funcionalidades de QGIS tienen esta opción disponible y se puede activar en la ventana del proceso cambiando a la pestaña Ejecutar procesamiento por lotes (ver figura 9.12).

Figura 9.12 - Pestaña de procesamiento por lotes en una ventana de funcionalidad QGIS

Para las capas ráster de cobertura terrestre, usaremos el procesamiento por lotes y como método de remuestreo del vecino más cercano. Para agregar una nueva capa, haga clic en el pictograma +. Para llenar automáticamente el CRS y los parámetros del método de remuestreo, haga clic en el botón de autocompletar en la parte superior de las columnas correspondientes y seleccione Rellenar. Cambie el nombre de los rásteres reproyectados agregando el código EPSG al final del nombre, por ejemplo, LandCover2015, se convertirá en landCover2015_5347. Configure sus parámetros como en la figura 9.13: CRS de origen: EPSG: 4326, CRS de destino EPSG 5347, método de remuestreo a utilizar: vecino más cercano (explicamos en el párrafo anterior por qué), valor de nodata para las bandas de salida: 255 (de la ventana de información, vemos el tipo de datos - yte - entero sin signo de 8 bits - lo que significa que el valor máximo puede ser 255), resolución de salida: 100 m (como los rásteres de cobertura terrestre iniciales). Después de configurar todos los parámetros, marque la casilla en la esquina izquierda de la ventana - Cargar capas al finalizar. Pulsa Ejecutar.

Figura 9.13a - Procesamiento por lotes para reproyectar los rásteres de cobertura terrestre

Figure 9.13b - Autofill output names

Figure 9.13c - Reprojected land cover rasters

Luego vienen los rásteres del modelo de superficie digital. Como se puede observar, para cubrir nuestra región de interés, necesitábamos varios mosaicos ráster. Cuando los archivos ráster se vuelven demasiado grandes, imagine un solo archivo DSM a 30m para Europa, que tiene más de 10 millones de kilómetros cuadrados, se dividen en mosaicos porque, en áreas más pequeñas, son más fáciles de administrar.

Aunque podríamos usar la herramienta de ajuste en modo por lotes para reproyectar todos los archivos ráster DSM, en una inspección visual, también se puede notar que las delimitaciones entre cada mosaico son bastante visibles, lo que lo hace, al menos, visualmente poco atractivo. Lo que sería útil es tener una visión completa del terreno, como un fenómeno continuo, como es, de hecho, sin interrupciones artificiales. Para obtenerlo, usaremos la herramienta de fusión GDAL, disponible en la barra de herramientas de Procesamiento para fusionar todos los rásteres DSM. Para abrirlo, vaya a Procesamiento ‣ Caja de herramientas y en la barra de búsqueda escriba combinar (ver figura 9.14).

Figura 9.14 - Encontrar la herramienta de combinación GDAL en la caja de herramientas de procesamiento

En la ventana de combinación que se abre, seleccione los archivos ráster DSM que queremos volver mosaico y haga clic en ejecutar. Guarde el resultado como DSM_mosaico. El resultado debería verse como en la figura 9.15.

Figure 9.15a - Selecting the ASTER layers to merge

Figure 9.15b - Parameters of the Merge processing algorithm

Figura 9.15c - Mosaico de todos los archivos DSM correspondientes a nuestra región de trabajo

Ahora, podemos proceder a reproyectar el mosaico - un archivo, en lugar de todos los archivos. Vaya a Ráster ‣ proyección ‣ Envoltura (reproyectar) y configure los parámetros conocidos:

• Fuente CRS EPSG 4326,
• Objetivo CRS: EPSG: 5347,
• Método de remuestreo: Vecino más cercano,
• resolución de archivo de salida - 30 m.
• Guarde el resultado como DSM_mosaico_5347.

En este punto, deberíamos tener todas las capas en el mismo CRS - EPSG 5347.

Figure 9.15d - Reproject Merged raster

Podemos hacer otra verificación para asegurarnos de que todos los rásteres con los que estamos trabajando estén proyectados y obtener información adicional sobre los datos ejecutando un proceso por lotes de Información de trama sobre todos. Para abrir la ventana de funciones, vaya a Ráster ‣ Miscelaneos ‣ Información de ráster. La ventana de información ráster de procesamiento por lotes debe verse como en la figura 9.16.

Figura 9.16 - Proceso por lotes para extraer información en un archivo HTML separado para múltiples capas ráster.

Un archivo HTML de información ráster debería verse como a continuación. Un archivo HTML se puede abrir con cualquier editor de texto o navegador web de su elección.

Driver: GTiff/GeoTIFF
Files: /home/malena/docus/okfn/capacitaciones/modulo9/modulo9_capas_qgis/DSM_mosaico_5347.tif
Size is 16546, 29702
Coordinate System is:
PROJCRS["POSGAR 2007 / Argentina 5",
    BASEGEOGCRS["POSGAR 2007",
        DATUM["Posiciones Geodesicas Argentinas 2007",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",5340]],
    CONVERSION["Argentina zone 5",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",-90,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",-60,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",5500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (X)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (Y)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["Argentina - 61.5°W to 58.5°W onshore"],
        BBOX[-39.06,-61.51,-23.37,-58.5]],
    ID["EPSG",5347]]
Data axis to CRS axis mapping: 2,1
Origin = (5202154.917641102336347,7014336.020335459150374)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 5202154.918, 7014336.020) ( 62d59'56.94"W, 26d58' 5.06"S)
Lower Left  ( 5202154.918, 6123276.020) ( 63d15'42.51"W, 34d59'35.07"S)
Upper Right ( 5698534.918, 7014336.020) ( 58d 0' 0.86"W, 26d59' 8.89"S)
Lower Right ( 5698534.918, 6123276.020) ( 57d49'29.65"W, 35d 1' 2.75"S)
Center      ( 5450344.918, 6568806.020) ( 60d31'12.08"W, 31d 1' 7.19"S)
Band 1 Block=16546x1 Type=Float32, ColorInterp=Gray
    Computed Min/Max=-76.000,251.000
  Minimum=-76.000, Maximum=251.000, Mean=72.929, StdDev=29.775
  NoData Value=255
  Metadata:
    STATISTICS_MAXIMUM=251
    STATISTICS_MEAN=72.929210388133
    STATISTICS_MINIMUM=-76
    STATISTICS_STDDEV=29.774935127844
    STATISTICS_VALID_PERCENT=67.29

Después de preparar los ráster reproyectándolos en los CRS que estamos trabajando y leyendo sus metadatos para comprender mejor los archivos, es hora de profundizar en los datos reales. Para lograrlo, calcularemos e interpretaremos los histogramas (consulte la sección Desglose de conceptos para obtener más detalles) de nuestros rásteres.

Para calcular un histograma, seleccione la capa ráster que le interese, ábrala haciendo doble clic en la ventana de diálogo Propiedades y vaya a Histograma (ver figura 9.17).

Figura 9.17 - Ventana de Histograma

Toca el botón Computar histograma y QIS lo computará automáticamente.

Después de computar el histograma, podemos ver que el mouse se convierte en una lupa. Es una herramienta que permite la inspección del histograma, viendo como varía la frecuencia de diferentes valores. Acercándose puede ver algo como en la figura 9.18.

Para volver a la vista completa, haga clic en la izquierda.

Figura 9.18 - Ampliación del histograma calculado HRSL_SantaFe_poblacion_5347

Más que simplemente ver la distribución de los valores numéricos de los píxeles, el histograma permite al usuario reclasificar los valores para la visualización del ráster. Para ello, utilice las 2 herramientas para señalar en el histograma los nuevos valores mínimo y máximo (consulte la figura 9.19).

Figura 9.19 - Selección de valores mínimos y máximos para reclasificar el ráster

Después de presionar Aplicar, el ráster se representará utilizando el nuevo rango mínimo-máximo seleccionado. Esta funcionalidad permite al usuario ignorar valores extremos que pueden estirar de forma anormal el ráster.

Aunque la siguiente herramienta que usamos es un complemento (consulte el módulo 1 para obtener detalles del complemento), lo consideramos muy útil para comenzar a trabajar con ráster. Nos referimos a la Value Tool (Herramienta de Valores) que permite la identificación inmediata de los valores de las celdas al pasar el cursor sobre capas ráster.

Vaya a Complemento ‣ Administrar e instalar complementos, busque el complemento Value Tool y haga clic en instalar. Luego, haga clic derecho en la barra de la ventana principal de QGIS para abrir todos los Paneles y Barras de Herramientas disponibles en su instalación de QGIS (ver figura 9.20).

Figure 9.20 - The Value Tool Panel

La practicidad de esta herramienta reside en su simplicidad de uso, con solo unos pocos clics, uno puede extraer muy fácilmente celdas de valor en las áreas exactas de interés. Además, permite esto para todas las capas ráster cargadas.

Value Tool tiene 3 pestañas: Table, Graph, Options (consulte la figura 9.21).

Figura 9.21 - Value Tool cargado - resaltar en la primera pestaña - Table

La primera pestaña - Table - presenta una lista de todas las capas ráster cargadas y los valores de las celdas, a medida que el usuario mueve el mouse. También existe la posibilidad de seleccionar con cuántos decimales se deben mostrar los valores. Si el mouse se desplaza fuera de la extensión de una capa ráster, en lugar de un valor, se mostrará un mensaje: “out of extent”.

La segunda pestaña, Graph, muestra en un gráfico unido todos los valores que lee en la posición del mouse. Permite al usuario insertar valores para mínimo y máximo en el eje Y, que es el eje de los valores de celda. El eje X enumera todas las bandas del ráster que muestra en la tabla, con los números de orden correspondientes: la banda que tiene el número 1 en la primera pestaña, también será la número 1 en el gráfico.

La tercera pestaña permite al usuario personalizar lo que Value Tool muestra: qué capas (todas, solo las visibles o las seleccionadas) y qué bandas mostrar.

Preguntas de evaluación

1. ¿Se pueden reproyectar capas ráster en otros sistemas de coordenadas?
   • Sí.
   • No
2. ¿Puede haber espacios entre los rangos de valores de un histograma calculado para una capa ráster?
   • Sí.
   • No.
3. ¿Pueden diferentes bandas del mismo conjunto de datos ráster tener diferentes resoluciones?
   • Sí.
   • No.

Fase 2: Introducción al trabajo con datos ráster

Ahora que hemos aprendido cómo extraer información básica en los conjuntos de datos ráster cargados, continuaremos con un procesamiento de datos ráster más detallado para obtener nuevos rásteres derivados y, en consecuencia, más información.

Como habrá notado, debido a la estructura del modelo de datos ráster, las capas que hemos cargado se están expandiendo sobre nuestra región de interés: la provincia de Pampanga. Eso no es deseable debido a varias razones, pero principalmente porque termina procesando más datos de los que realmente necesita, lo que se traduce en mayores necesidades de almacenamiento y procesamiento informático. Por eso, antes de seguir adelante con cualquier otro paso, nos aseguraremos de procesar exactamente la cantidad de datos que necesitamos. Tenga en cuenta, dado que comenzará a trabajar en sus propios conjuntos de datos, que el tamaño de sus archivos es un factor importante cuando se trata de tiempos de procesamiento. Cuanto más grandes sean los archivos, se requerirá más tiempo. Debido a la estructura de datos del modelo (ráster frente a vector), los archivos ráster suelen ser mucho más grandes.

Como ya habrá notado, los conjuntos de datos que cargamos en un SIG, en nuestro caso en QGIS, se pueden procesar juntos incluso si son de diferente naturaleza, como unir tablas csv a capas vectoriales para agregar información a las geometrías. Lo mismo se aplica a los datos ráster y vectoriales, como veremos.

Para trabajar solo con capas ráster que son relevantes para nuestra provincia de Santa Fe, usaremos la capa de extensión vectorial (limite_provincial_santafe) para cortar / recortar todas las capas ráster relevantes. Vaya a Ráster ‣ Extracción ‣ Recortar ráster por capa de máscara (consulte la figura 9.22).

Figura 9.22 - Usando una máscara vectorial para extraer los datos ráster en una región específica

Dado que trabajaremos con 12 capas ráster - las 10 capas de Cobertura Terrestre, el modelo de superficie digital y el HRSLl, usaremos el procesamiento por lotes para recortar todas las capas en una vez. Tenga en cuenta que si ha omitido el paso de reproyección, tiene capas en diferentes proyecciones y el algoritmo no funcionará o producirá resultados inesperados.

La configuración de la ventana de procesamiento por lotes debería verse como en la figura 9.23.

Figura 9.23 - Procesamiento por lotes recortando todas las capas ráster requeridas por geometría de la provincia de Santa Fe

La configuración de los parámetros es la siguiente:

• capa de máscara: limite_provincial_santafe,
• tanto el CRS de origen como el de destino es EPSG 5347,
• seleccione sí para: hacer coincidir la extensión del ráster recortado con la capa de máscara y mantener la resolución de la capa de entrada.
• Tenga en cuenta que para DSM_mosaico también seleccionaremos sí para crear una banda alfa de salida. Cargue las capas al finalizar.

Si todo salió bien, su ventana principal de QGIS debería verse como en la figura 9.24.

Figura 9.24 - Capas ráster recortadas por el contorno de la provincia de Santa Fe.

Ahora, imagine que tiene que presentar un informe sobre dónde vive la mayoría de las personas, pero teniendo en cuenta la altitud [^ 4]. Debes saber cuántas personas viven entre 0 y 100 m de altitud en la provincia de Santa Fe. Hay algunos elementos a considerar. En primer lugar, cuáles son los datos que utilizaremos y cuáles son sus características. Para la población, tenemos los datos de la capa de asentamiento de alta resolución y para el alivio, tenemos el ALOS World 3D - 30m (AW3D30). Ambas capas ráster tienen una resolución espacial de 30 m, lo que nos permite pasar a otras consideraciones. El alivio es un fenómeno continuo, la expansión de la población no lo es, sin embargo, el informe no tendría sentido que lo hiciera el píxel de 30 metros. Necesitamos identificar todos los píxeles con valores de celda de 0 a 100. Considerando el histograma de DSM_mosaico para nuestra región de interés, hemos visto que la mayoría de los valores de celda están entre 0 y 100 m. Podemos proceder a realizar un mapa de relieve básico en base a las siguientes divisiones:

1. 0 - 25m
2. 26 - 50m
3. 51 - 75m
4. 76 - 100m
5. 101m++.

Usando sus conocimientos adquiridos en el módulo 4, puede diseñar la capa DSM por estas categorías. Su mapa debería verse como en la figura 9.25.

Figura 9.25 - Representación DSM_mosaico_recortado

Para calcular el número de personas en base a los datos raster HRSL que viven hasta 50 metros de altura en la provincia de Santa Fe, debemos ver qué píxeles caen en cada una de esas categorías. Para hacer eso, usaremos Raster Calculator. Esta es una funcionalidad que permite al usuario realizar cálculos sobre la base de los valores de píxeles de la trama existentes. Los resultados se escriben en una nueva capa ráster en un formato compatible con GDAL [^ 5].

There are several ways to open the raster calclulator in QGIS. You can do so from the Menu bar Raster ‣ Raster Calculator or by searching raster calculator on the Processing Toolbox or Locator bar. If we run the Raster Calculator under Raster analysis in the Processing Toolbox, the window in Figure 9.26b should appear.

Figure 9.26a - Opening the Raster calculator

Figura 9.26b - Calculadora ráster

En esta ventana podemos reconocer las operaciones detalladas en la sección Conceptos, sustracciones, adiciones, comparaciones y todas las demás (ver página 3). A continuación, ‘cortaremos’ nuestra capa ráster cortado_DSM_mosaico_5347 para extraer solo los píxeles con valores de hasta 50 metros. Sabemos que las celdas de valor de cortado_DSM_mosaico_5347 representan datos numéricos continuos (no valores discretos, como LandCover). Por lo tanto, la operación que debemos emplear en este caso es una comparación de valores de una celda <50 metros. Para obtenerlo, escribiremos lo siguiente en la Calculadora:

"cortado_DSM_mosaico_5347@1" <= 50

Se puede observar la convención de nomenclatura para los rásteres: lo que viene antes de @ es el nombre de la capa ráster, lo que viene después de @ es el número de la banda. Guarde el resultado de la calculadora como DSM_menos50.tiff. Su calculadora de ráster debería verse como en la figura 9.27.

Figura 9.27 - Insertar una fórmula en la Calculadora ráster.

Su resultado debería verse como en la figura 9.28.

Figura 9.28 - Resultado de identificar todos los valores de píxeles que están por debajo de 50 metros usando la Calculadora Ráster

Como podemos ver en la Tabla de Contenidos, la capa ráster que hemos obtenido tiene solo 2 valores: 0 y 1. Eso es porque hemos usado una operación racional, una comparación, por lo tanto, cada píxel que está por debajo de 50 metros recibió un valor = 1 y todos los anteriores, un valor = 0. Podemos probar esto usando Value Tool. La figura 9.29 presenta solo píxeles de valor 1, en otras palabras, los píxeles que nos interesan para nuestro ejercicio.

Figura 9.29 - Distribución espacial de todos los píxeles de valor 1, es decir, con una altitud inferior a 200 metros

Yendo más allá, podemos mostrar la distribución espacial de la población a la resolución espacial de 30 m solo en esta región geográfica específica que hemos seleccionado - provincia de Santa Fe , por debajo de 50 m. Para hacer eso, volvemos a emplear la Calculadora ráster.

La fórmula es bastante simple, dado que todos los valores de celda DSM que nos interesan tienen el valor 1.

Abra la Calculadora e inserte la siguiente fórmula:

"dsm_menos50@1"*"HRSL_SantaFe_poblacion_5347@1"

Figura 9.30 - Uso de la calculadora ráster para identificar la distribución de la población clases basadas en una altitud de hasta 50 m

A diferencia del uso anterior de la calculadora ráster, hemos utilizado 2 conjuntos de datos ráster diferentes para obtener el resultado deseado, sin embargo, observará que incluso si hay píxeles que caen fuera de DSM_menos50.tiff en el resultado, su valor es 0. Utilice Value Tool para comprobar (consulte la figura 9.31).

Figura 9.31 - Usando Value Tool para verificar los resultados de la Calculadora ráster

Puede ver que incluso si hrsl_menos50 tiene valores en esta ubicación específica del mouse, el ráster obtenido con la Calculadora ráster tiene el valor 0.

A continuación, presentamos la distribución espacial de la población que vive por debajo de 50 m en la provincia de Santa Fe. Para elegir una clasificación adecuada, calculamos el histograma. Podemos notar que la mayoría de los valores se encuentran entre 0,1 y 8,7 personas por 30 m. La clasificación que hemos elegido es visible en la figura 9.32.

Figura 9.32 - Distribución de la población que vive por debajo de los 50 m en la provincia de Santa Fe, representada con una resolución de 30 m.

Si estamos interesados ​​en el número total de personas que viven por debajo de 50 m en la provincia de Santa Fe y no en la distribución geográfica por 30 m de resolución espacial, entonces necesitamos sumar todos los valores de píxeles de la capa ráster hrsl_menos50. Una forma de obtener este número es transformar el dsm_menos50 de ráster a vector y luego [….]

Para hacer eso, vaya a Ráster ‣ Conversión Poligonizar (Ráster a Vector) (ver figura 9.33).

Figura 9.33a - Conversión de ráster a vector

Recuerde que esta capa ráster tenía solo 2 valores - 0 y 1, así que elija como parámetro por el cual construir el vector el DN (número digital).

Figure 9.33b - Raster to vector conversion parameters

Note that it’s possible that you will have multiple features in the vectorized layer instead of just 2. In order to just have 2 features, you can run the Dissolve algorithm and choose DN as the Dissolve Field. This will result in having just 2 features in the layer. Also, don’t forget to fix invalide geometries using the Fix invalid geometries algorithm if you encounter some geometry issues. Your result should look like in figure 9.34.

Figura 9.34 - Resultado de convertir un conjunto de datos ráster en uno vectorial.

Elimina la geometría de valor 0.

Para encontrar nuestra respuesta usaremos Estadísticas de zona. Para encontrar rápidamente esta funcionalidad, abra la Caja de herramientas de procesamiento y escriba en el cuadro de búsqueda “zonal” (consulte la figura 9.35).

Figura 9.35 - Identificación de estadísticas zonales en la caja de herramientas de procesamiento

En la ventana que se ha abierto, seleccione los parámetros como en la figura 9.36. A medida que se calculan las estadísticas, seleccione: recuento, suma, mínimo y máximo.

Figura 9.36 - Configuración de los parámetros para Estadísticas zonales

La capa resultante es una capa vectorial que tiene como atributos las estadísticas que hemos seleccionado (ver figura 9.37).

Figura 9.37 - Capa vectorial resultante de Estadísticas zonales

Y con este paso final, respondimos a nuestro ejercicio, cuántas personas (y dónde) viven por debajo de los 50 m en la provincia de Santa Fe.

Preguntas de evaluación

1. ¿Se pueden recortar los datos ráster?
   • Sí.
   • No.
2. ¿Se pueden usar conjuntos de datos vectoriales cargados en el proyecto QGIS en Raster Calculator?
   • Sí.
   • No.
3. ¿Se pueden convertir conjuntos de datos ráster en vectoriales?
   • Sí.
   • No.

Fase 3: Trabajo con datos raster y vectoriales.

En la fase anterior, hemos visto cómo podemos procesar 2 conjuntos de datos ráster para derivar nueva información. Hemos utilizado el modelo de superficie digital y la capa de asentamiento de alta resolución para averiguar cuántas personas viven por debajo de los 200m en la provincia de Pampanga. Antes de hacer cualquier análisis, nos aseguramos de que los conjuntos de datos estuvieran en la misma proyección y, además, si los rásteres tienen la misma resolución espacial para que los resultados que obtuvimos sean viables. Al referirse al sistema de referencia de coordenadas, el razonamiento es claro, pero ¿por qué la misma resolución espacial?

Recordando la resolución espacial, es el tamaño de la superficie del suelo medido en unidades de longitud, en otras palabras, el tamaño del píxel medido en el suelo. Si un ráster tiene una resolución de 30 m, eso significa que el objeto lineal más pequeño que pudimos detectar en esa imagen es de 30 m, más pequeño y no pudimos detectarlo. Continuando con la analogía, podemos compararla con la escala de un mapa. Si un mapa tiene una escala de 1: 25000, eso significa que 1 unidad de longitud en el mapa representa 25000 unidades en el suelo, es decir, 1 cm son 25000 cm, 1 cm en el mapa equivale a 250 m en el suelo. Por ejemplo, una carretera de 2 km tendría 8 cm en el mapa.

¿Por qué es tan importante cuando se trabaja con conjuntos de datos ráster? La figura 9.38 podría ofrecer una explicación.

Figura 9.38 - Ejemplo de diferentes resoluciones específicas para diferentes imágenes satelitales - Landsat y SPOT - para la misma área

(Crédito de la foto: Congedo, L. y Munafò, M, (2013) Evaluación del cambio de cobertura terrestre mediante sensores remotos: objetivos, métodos y Resultados, Roma: Universidad Sapienza. Disponible en: http://www.planning4adaptation.eu/ La

Figura 9.38 detalla la relación entre la resolución de las imágenes satelitales y la información de cobertura terrestre extraída de estas imágenes capturadas. Recuerde como detallamos al principio que el valor de la celda de píxel se atribuye a toda el área que cubre, pero eso no significa que esa sea la realidad en el terreno. Estas representan decisiones tomadas por los expertos de EO que derivan varios productos basados ​​en imágenes de observación de la Tierra, todos documentados en artículos revisados ​​por pares y descripciones de algoritmos. Más explicaciones están más allá del alcance de este módulo, pero es importante tener en cuenta la relación entre lo que un sensor a bordo de un satélite captura y los productos que utilizamos.

Volviendo a nuestra región de interés, la provincia de Santa Fe, podemos probar estas diferencias con los datos que tenemos a mano. Hemos cargado en nuestro proyecto QGIS las 5 capas ráster de LandCover durante 5 años: 2015, 2016, 2017, 2018 y 2019. A continuación, cargaremos un mosaico de imágenes Sentinel-2 [^ 6]. Cargaremos la capa WMS EOX Sentinel-2 sin nubes, disponible aquí. Recordando el módulo 2, para agregar una capa WMS, vaya a Capa ‣ Agregar capa ‣ Agregar capa WMS / WMTS.

Cuando se abra la ventana de agregar, use los siguientes parámetros:

Name: EOX Sentinel-2
URL: https://tiles.maps.eox.at/wms?service=wms&request=getcapabilities

Figura 9.39 - Adición de una capa WMS a QGIS

Después de conectarnos a la capa WMS recién agregada, cargaremos la capa llamada Sentinel-2 cloudless layer para 2019 por EOX - 4326 en QGIS. Después de hacer zoom a la extensión de la región de interés, la ventana de su mapa debería verse como en la figura 9.40.

Figura 9.40 - Capa sin nubes Sentinel-2 para 2019 por EOX - 4326 para la provincia de Santa Fe

Aunque los productos LandCover se han obtenido utilizando otros datos satelitales (Proba-V), compararemos las 2 capas para que podamos tener una idea de las diferentes resoluciones significar. Recuerde que el producto LandCover está a 100 m y las imágenes de Sentinel 2 a 10m. Para lograrlo, abriremos cortado_LandCover2019_5347 y la capa WMS. Para hacer comparaciones entre 2 capas, usaremos un nuevo plugin que deberás instalar. Por lo tanto, vaya a Complementos ‣ administrar e instalar complementos y escriba en el cuadro de búsqueda MapSwipe Tool. Una vez que lo instale, debería aparecer como un nuevo pictograma en su barra de herramientas ().

Figura 9.41 - Comparación de 2 capas ráster usando el complemento de la herramienta MapSwipe

Para activar la herramienta MapSwipe, haga clic en ella mientras la capa ráster que desea cubrir está seleccionada en el Panel de capas. Las diferencias de resolución son obvias, así como el hecho de que el producto de satélite (Land Cover) se ha desarrollado utilizando una imagen de satélite (PROBA-V) de una resolución más gruesa. Sin embargo, las clases generales más grandes están bien identificadas, como se puede ver en la figura 9.42.

Figura 9.42 - LandCover2019 obtenido de PROBA-V (100 m) en la parte superior del mosaico Sentinel 2 (30 m)

Al agregar el HRSL al mapa, se mostrará una buena coincidencia entre HRSL y LandCover. El espacio urbano está representado en rojo y, como se puede ver en la figura 9.43, está cubierto casi por completo por la capa HRSL.

Figura 9.43 - HRSL agregado en la parte superior de cortado_LandCover2019_5347.

Sin embargo, al hacer zoom, puede ver la diferencia en la resolución entre los 2 productos ráster, como en la figura 9.44.

Figura 9.44 - Diferencia en la resolución espacial entre HRSL (30m - blanco) y cortado_LandCover2019_5347 (100m - colores verde, naranja y violeta)

Ahora, si se hiciera algún análisis usando estos rásteres, estos resultados no serían viables, porque estaríamos comparando valores que se aplican a diferentes resoluciones. Como fase de preprocesamiento, el usuario debe volver a muestrear uno de los 2 para que coincida.

El remuestreo se refiere a cambiar los valores de celda debido a cambios en la cuadrícula de celdas ráster y solo hay 2 opciones: (1) el muestreo superior se refiere a los casos en los que estamos convirtiendo a celdas más pequeñas / de mayor resolución y (2) el muestreo reducido es un remuestreo a una resolución más baja / tamaños de celda más grandes.

Imaginemos el siguiente ejercicio. Necesitamos identificar los números de población para cada categoría de cobertura terrestre que hemos definido en la provincia de Pampanga. Como se explicó anteriormente, necesitamos preprocesar los datos que tenemos para obtener resultados viables de nuestro análisis, es decir, en nuestro caso, necesitamos llevar nuestros dos conjuntos de datos ráster a la misma resolución espacial. Como se detalla anteriormente, podemos aumentar o disminuir las dimensiones de los píxeles. Debe destacarse que volver a muestrear, con escala ascendente o descendente, implicará un proceso de interpolación (consulte la página 12 para obtener más detalles), por lo que el resultado introduce un error estadístico. La práctica habitual es volver a muestrear todos los rásteres para que se correspondan con el ráster con la resolución más baja, pero nuevamente esta decisión debe tomarse teniendo en cuenta todos los factores. Detallar el proceso de toma de decisiones para remuestrear capas ráster excede con mucho el alcance de este plan de estudios.

Debe destacarse una diferencia entre los 2 productos: el producto LandCover cubre toda el área de la extensión considerada, a diferencia del producto HSRL, donde la capa ráster contiene estrictamente la celda donde existen valores superiores a 0. Esta situación plantea problemas cuando se vuelven a muestrear los valores de las celdas de interpolación, porque sin importar el método de interpolación elegido, eso tomaría en consideración los píxeles circundantes, siguiendo un algoritmo específico bien definido y, en esta situación particular, los píxeles de borde no están en el borde de la región de estudio, pero dentro de ella. Por lo tanto, en nuestro caso de demostración, consideraremos aumentar el muestreo del producto Land Cover de una resolución de 100m a 30 m para que coincida con la resolución del producto LandCover. El método de remuestreo que elegimos es crucial, ya que los resultados pueden variar significativamente. Para fines de demostración, volveremos a muestrear el producto LandCover utilizando 2 métodos diferentes: Vecino más cercano y Modo.

Para volver a muestrear, vaya a Ráster ‣ Proyecciones ‣ Ajustar (reproyectar). En la ventana de funcionalidad, configure los siguientes parámetros:

• capa de entrada: cortado_LandCover2019_5347,
• CRS de origen y CRS de destino: EPSG: 5347,
• Método de remuestreo: Vecino más próximo,
• Sin datos: 255, resolución de archivo de salida: 30,
• Tipo de datos de salida: usar datos de capa de entrada tipo,
• extensión de georeferencia: seleccione la capa cortado_LandCover2019_5347.

Guarde la capa de salida como LC2019_NearestNeighbour.

Figura 9.45 - Remuestreo de la capa de Cobertura Terrestre

Siga los pasos exactos, excepto que para el parámetro del método de remuestreo elija Modo. Guarde la capa de salida como LC2019_Mode.

Ahora, comparemos los resultados: consulte las figuras 9.46a y 9.46b.

Figura 9.46a - Producto de remuestreo de cobertura terrestre usando el método de remuestreo del vecino más cercano

Figura 9.46b - Producto de remuestreo de cobertura terrestre usando el método de remuestreo de modo

Ambos rásteres tienen aplicada la misma simbología y podemos observar que en la figura 9.45 hay valores que no caen en ninguna categoría, los píxeles no se muestran. Sin embargo, sabemos que el producto Land Cover es una capa uniforme: no tiene espacios entre las mismas categorías definidas. Profundicemos y observemos los histogramas de las tres capas. Para eso, vaya a Propiedades ‣ Histograma y elija Prefs / Acciones para mostrar solo la Banda seleccionada. Guarde el histograma haciendo clic en el icono de guardar en el lado derecho de la ventana (consulte la figura 9.47).

Figura 9.47 - Mostrar valores solo para una banda seleccionada en el histograma.

La figura 9.48 (a), (b) y (c) presenta los 3 histogramas de interés.

Figura 9.48 - Histogramas de (a) cortado_LandCover2019_5347 (100m), (b) LC2019_NearestNeighbour y (c) LC2019_Mode

Podemos observar las diferencias en la distribución de valores para los 3 conjuntos de datos, énfasis en (b), donde usamos el vecino más cercano método de remuestreo, que resulta, como era de esperar, en valores de píxeles calculados y, por lo tanto, da como resultado valores que no se corresponden con ningún valor en el producto Land Cover (ver tabla 1, página 7). De ahí los espacios en blanco: píxeles sin color asignado en la figura 9.45. Como se infiere, el método de remuestreo de modo, o método de remuestreo mayoritario, selecciona el valor que aparece con más frecuencia.

En este punto, las ventanas de su mapa deberían verse como en la figura 9.49.

Figura 9.49 - Los 2 productos ráster: Land Cover 2019 y HRSL superpuestos.

Volviendo a nuestro ejercicio, los requisitos eran identificar los números de población para cada categoría de cobertura terrestre que hemos definido en la provincia de Pampanga. En este punto, hemos preprocesado nuestros datos ráster, para tenerlos en el mismo sistema de coordenadas, con la misma resolución espacial. Continuaremos con un algoritmo de conversión, transformaremos el dataset ráster Land Cover en un dataset vectorial - tipo polígono. Esto nos permitirá identificar más fácilmente el recuento de población para cada categoría de cobertura terrestre.

Para conversiones, de ráster a polígono, así como de vector a polígono, vaya a Ráster ‣ Conversión y aquí tenemos más opciones. Elegiremos Polygonize (Raster to vector). Convertiremos a vector el último dataset ráster que hayamos obtenido: LC2019_Mode. Su resultado debería verse como en la figura 9.50b.

Figura 9.50a - Poligonización de la capa ráster LC2019_Mode (30m)

Figura 9.50b - Resultado de poligonizar la capa ráster LC2019_Mode (30m)

Como podemos observar, cada píxel - cada grupo de píxeles adyacentes con el mismo código de categoría (ver tabla 1, página 7) se convirtió en un polígono. Sin embargo, no estamos interesados ​​en cada región distinta, sino en toda la categoría. Así, disolveremos la capa vectorial que hemos obtenido usando la identificación de clase como atributo común (Vector ‣ Herramientas de geoprocesamiento ‣ Disolver - para más detalles, ver el módulo 8).

Adjuntando el mismo estilo que para el ráster, notaremos que las mismas clases vectoriales corresponden a las del ráster (ver figura 9.51).

Figura 9.51 - LC2019_Mode poligonos

Como podemos observar, los píxeles de la HRSL, como se esperaba, no están perfectamente contenidos en un polígono del producto de cobertura terrestre (ver imagen 9.52).

Figura 9.52 - Discordancia de píxeles HRSl y el conjunto de datos vectoriales de cobertura terrestre

Para eliminar este inconveniente, también vectorizamos la capa HRSL, solo que esta vez transferiremos el valor de la celda del píxel a un punto: el centro geométrico de cada píxel.

Vaya a Procesamiento ‣ Caja de herramientas. En la barra de búsqueda, escriba las palabras clave punto ráster y elija el Píxeles ráster a puntos algoritmo (consulte la figura 9.53a). Guarde la salida como HRSL_puntos.

Figura 9.53a - Píxeles ráster a puntos

Figura 9.53b - Ejecución del algoritmo Píxeles ráster a puntos

Teniendo en cuenta la extensión de su área de estudio, esta operación puede prolongarse significativamente en el tiempo. La Figura 9.54 muestra cuántos puntos hemos obtenido para la provincia de Santa Fe.

Figura 9.54 - Carga de datos vectoriales puntuales obtenidos

El número de características es considerablemente alto y sin importarlo a una base de datos, cualquier tipo de procesamiento o visualización requeriría demasiado tiempo. En este tipo de situaciones, la solución razonable es dividir los conjuntos de datos que tenemos que procesar en fragmentos manejables. Por lo tanto, consideraremos procesar los cálculos necesarios en áreas más pequeñas y bien definidas. Para dividir la capa HRSL usaremos la opción para crear un VRT. Seleccione la capa HRSL y elija Exportar como .. En la nueva ventana, marque la opción Crear VRT y configure los siguientes parámetros: busque una carpeta donde se exportará el ráster dividido, CRS: EPSG: 5347, mosaicos VRT: máx. columnas 1000, filas máximas: 1000 (consulte la figura 9.55).

Figura 9.55 - Creación de un archivo VRT con mosaicos ráster específicos

Después de exportar, cargue todos los mosaicos ráster en su proyecto QGIS. El resultado debería verse como en la figura 9.56.

Figura 9.56 - Mosaicos ráster para la HRSL del departamento de Rosario

A continuación, volveremos a ejecutar el algoritmo de píxeles ráster a puntos para cada uno de los mosaicos ráster. Debido a que tenemos varios mosaicos, usaremos la función de procesamiento por lotes (ver figura 9.57).³

Figura 9.57 - Ejecución de ráster a puntos para todos los mosaicos ráster

Los puntos vectoriales resultantes deben verse como en la figura 9.58.

Figura 9.58 - Conjunto de datos de vectores de puntos para la capa HRSL con valores de píxeles almacenados en la columna VALOR

Al observar más de cerca los resultados del algoritmo, podemos ver que los puntos vectoriales caen exactamente en el centro del píxel del que extraen el valor ( ver figura 9.59.).

Figura 9.59 - Verificación de los valores de los datos vectoriales de puntos frente a los valores de los píxeles de la trama

Para resolver el ejercicio, se debe calcular la suma de los valores de los puntos extraídos del producto HRSL que se encuentran dentro de cada polígono del producto de cobertura terrestre. Para hacer eso, usaremos una función que está disponible en la Calculadora de campo: aggregate(). Esta función es bastante poderosa ya que hace uniones espaciales sobre la marcha permitiendo varios cálculos. En nuestro caso, la función debe identificar los puntos que caen en cada polígono y luego sumar los valores de esos puntos. Para hacer eso, vaya a la tabla de atributos del conjunto de datos vectoriales LC2019_Mode⁴ y abra la calculadora de campo. Creando un nuevo campo, introduzca la siguiente expresión:

aggregate( layer:= ‘hrsl_rosario1’, aggregate:=’sum’, expression:=VALUE, filter:=intersects($geometry, geometry(@parent))

)

Donde capa es el nombre del conjunto de datos del que queremos extraer información (en nuestro caso, el conjunto de datos puntuales que contiene los valores de píxeles de la capa ráster HRSL), agregado: indica la acción que se realizará una vez que se confirme la unión espacial (suma, recuento, media, mediana, concatenar, etc.), expresión: indica de qué columna queremos extraer datos, filtro: indica la función de geometría (intersección, dentro, etc.) (ver figura 9.60).

Figura 9.60 - Uso de funciones integradas con calculadora de campo

Los resultados, como se esperaba, se guardan en la tabla de atributos del LC2019_Mode. Y así, se ha completado el ejercicio.

¡Por favor, tenga cuidado! ¡Este último paso puede ser excesivamente largo dependiendo del volumen de sus datos! Pruebe para encontrar sus mejores dimensiones de mosaico.

Si usa la Calculadora de campo y la función agregada no funciona, también podemos usar el algoritmo Unir atributos por ubicación (resumen) para calcular la suma de los valores de los puntos dentro de las características LC_Mode_vector.

Figura 9.61 - Parámetros del algoritmo Unir atributos por ubicación (resumen)

Figura 9.62 - Salida del algoritmo Unir atributos por ubicación (resumen)

En esta tercera fase de nuestro módulo 9, hemos trabajado con raster, así como con datos vectoriales. Hemos analizado más de cerca lo que significa el remuestreo y cuáles son las implicaciones de las características de los rásteres, como el sistema de referencia de coordenadas, la resolución espacial, etc. La mayoría de las veces, se debe trabajar con conjunto de datos ráster y vectoriales y, por lo tanto, es importante saber que es posible de varias formas. También probamos varias conversiones, de ráster a vector y viceversa. Sin embargo, en nuestro procesamiento no debemos perder de vista los fenómenos que estamos estudiando, así como la escala inicial a la que se han recolectado los datos, ya sea vectorial o raster. Al realizar conversiones, es importante conocer la correspondencia entre la escala de un mapa y la resolución de la trama. Según el cartógrafo Waldo Tobler “la regla es dividir el denominador de la escala del mapa por 1000 para obtener el tamaño detectable en metros. La resolución es la mitad de esta cantidad “. En otras palabras, la fórmula sería: escala del mapa = resolución de la trama (en m) X 2 X 1000.

Las soluciones dadas a nuestros ejercicios en este módulo, como también en el módulo 8, son solo una forma de resolver los requisitos. A medida que se avanza en el estudio de las herramientas de código abierto para geoespacial, se descubre que hay varias posibilidades de llegar al mismo resultado, algunas mejores que otras.

Preguntas de evaluación

1. ¿Cuál es el nombre del proceso mediante el cual la resolución de un dataset ráster se puede aumentar o disminuir?
   • Remuestreo.
2. ¿Qué nos muestra un histograma?
   • La frecuencia de los valores de los píxeles, organizados en rangos de valores adyacentes.
3. ¿Se puede convertir un dataset ráster en un polígono? ¿Qué acerca del otro camino alrededor?
   • Sí y sí.

Referencias:

Centro para la Red Internacional de Información sobre Ciencias de la Tierra - CIESIN - Universidad de Columbia. 2018. Gridded Population of the World, Versión 4 (GPWv4): Densidad de población, Revisión 11. Palisades, NY: Centro de Aplicaciones y Datos Socioeconómicos de la NASA (SEDAC). https://doi.org/10.7927/H49C6VHW.

Métodos de remuestreo: https://gisgeography.com/raster-resampling/

Notas

Notes

1. https://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm ↩

2. https://lcviewer.vito.be/download ↩

3. Se recomienda hacer el proceso por lotes en sólo una porción del territorio de la provincia para evitar procesar demasiadas capas (Departamento de Rosario) ↩

4. Se puede usar un recorte por departamento al igual que los píxel a raster ↩

This site is open source. Improve this page.