Laboratorio 1 - Diplomado Geómatica Aplicada U. de Chile Realizado por Rodrigo Aros Bustamante, Geógrafo-Especialista GIS, Universidad de Concepción, email: , Concepción, 2023.

1. Introducción

En este informe se presentará los resultados de ánalisis para caso de estudio proyecto de inversión en energía, donde se pretende construir en el área de la cuenca de abastecimiento (cuenca H). En donde se debe evaluar la instalación de torre de vigilancia y su mejor ubicación, basado en visibilidad de bosques y matorrales. Posterior a esto un ánalisis de los índice NDVI y Shannon, para subusos de suelo de la comuna de Constitución, comuna y ciudad de la provincia de Talca en la VII región del Maule de Chile.

2. Metodología

Se cuenta con información geográfica la cual se encuentra georreferenciada con EPSG: 32718, correspondiente a SRC: WGS84 UTM Zona 18 S. Para el caso se utiliza información vectorial y raster. Se construye MDE a partir de interpolación de puntos que contiene información de altura(h). Posterioremente se aplica eliminación de sumideros(Wang & Liu). Utilzando software Saga GIS se aplica Upslope Area (Interactivo) con método “Deterministic 8” O’Callaghan, J.F. / Mark, D.M. (1984), para marcar punto de cuenca de abastecimiento denominada “H”. En una etapa posterior se aplica Viewshed con altura obervador = 65 mts, utilizando MDE (No sinks). Luego se normaliza salida con calculadora raster.Finalmente Se corrigen Geometrías con algoritmo de QGIS. Ademas de aplicar correciones autómaticas en Propiedades-> Digitalización -> Correciones Automáticas.

En relación a la parte 2 del trabajo se crea un srcript en Rstudio y se aplica de la siguiente forma:

##### Cargo librerias ,instaldas por tanto no es necesario usar install.packages####

#library(tidyverse)
#library(raster)
#library(sf)
#library(exactextractr)
#library(ggplot2)


##### Seteo mi WD #####
#setwd("C:/Users/Rodrigo Aros/Desktop/Laboratorio/Nueva carpeta/Nueva carpeta/R")

#getwd() # devuelve mi WD

#### procedo a calcula NDVI ####

#imagen <-stack("Landsat_8_2016.tif")
#class(imagen) # consulto que tipo de objeto es img


#RED<-imagen$Landsat_8_2016_3
#NIR<-imagen$Landsat_8_2016_4

#NDVI<-(NIR-RED)/(NIR+RED)
#plot(NDVI, main=" Indice NDVI") # ploteo objeto NDVI con titulo "Indice NDVI"


#######################################################################
#### cargo Tabla Indice de Shanon #####

#IndiceShanon <-  "biodiversidad_subusos.csv" #traigo path desde win

#bio_div <- read.csv2(IndiceShanon)  #lectura csv
#class(bio_div)
#colnames(bio_div) # consulta nombre de columnas




#### cargo shape constitucion ####

#constitucion <- st_read("usos_constitucion.shp")

#### join con bio_div para sumar indice de shanon a Shape constitucion

#consti2<- left_join(constitucion, bio_div, by = "ID_SUB_USO")


#### Paso a calcular promedio ####

#PromNDVI <- exactextractr::exact_extract(NDVI,consti2, 'mean')


#### aplico mutate para calcular promNDVI ya que exactetractr me devuelve vector
#### con la funciòn mutate calculo el campo en el propio objeto SF llamado consti2 #####
#constitucionNDVI <- consti2 %>%
  #mutate(PRomNDVI = exact_extract(NDVI,consti2,'mean'))

#st_write(constitucionNDVI, "constiNDVI.shp") #escribo mi archivo de salida

#a<- st_read("ConstiNDVI.shp") # creo objeto


#plot(a['SHANNON']) # ploteo de acuerdo a INDICE SHANON


#a$SHANNON<- as.numeric(as.character(a$SHANNON)) #coerciono Shannon ya que se encontraba como character lo paso a numeric
################################################################################
############### Resultados graficos y correlacion entre variables  #############

#plot(a['SHANNON']) # ploteo nuevamente con indice shannon en numeric

#cor(a$SHANNON,a$PRomNDVI) # consulto correlacion en las variables

##### aplico boxplot a Prom NDVI #####

#boxplot(a$PRomNDVI~a$ID_SUB_USO, main = "Boxplot de PRomNDVI por ID_SUB_USO",
        #xlab = "ID_SUB_USO", ylab = "PRomNDVI") 

#### aplico boxplot a shanon #####
#boxplot(a$SHANNON~a$ID_SUB_USO, main = "Boxplot de Indice Shanon por ID_SUB_USO",
        #xlab = "ID_SUB_USO", ylab = "Indice Shanon")

3. Resultados

Resultados Parte 1 En el proceso la evaluación resultante se pudieron constatar que el algoritmo kriging de QGIS conlleva una serie de paràmetros que modifican el resultado esperado. En Capa resultante se obtuvo resultado óptimo con Cell size = 50, maximum search distance = 2000 m y aplicando “all points within search distance”. Por otro lado se pudo observar que el modelo digital de elevación puede variar bastante aplicando la eliminación de sumideros, esto observado de acuerdo a sus histogramas. Una vez aplicado kriging se pudo observar que con filtro Fill sinks de Wang Liu, los valores de pixel previos a su aplicaciòn (curva gris) tiende a aumentar el valor de pixel máximo (curva roja), pero a esos valores altos previos los distribuye en otras zonas ” aplanando” de alguna forma la distribución.

a. Histograma MDE Con Sumideros

test

b. Histograma post Fill Sinks

test
Otro ejercicio que se realizò fue comparar el comportamiento del algoritomo viewshed con otro SIG, a modo de ejercicio adicional y como comprobación de resultado. Para evaluaciòn de Punto B se obtuvo similares resultados. Tanto desde QGIS como de ArcMap. utilizando párametros similares en la aplicación de la herramienta. Los parámetros para esta comparativa fueron los soguientes:

Viewshed Punto B

Input = MDE FInal (No sinks)

nbanda = banda 1

Localizacion Observador = Ubicación pto B

Altura Observador = 65 m

Altura objetivo= 1

Distancia maxima = 19000 m

Para el resultado zona en color rosa es la visibilidad aplicada a todo el extent de MDE Final (No sinks). Para QGIS se utiliza la poligonización de el output de viewshed y se deja sin relleno mostrando el resultado obtenido y asé poder comparar. En la generalidad es bastante similar, sin embargo, realizando exploración visual con zoom existen distinciones en los resultados. No se realza comparaciòn estadistica mas profunda, debido a que no es el objetivo del laboratorio, pero si como una inquietud futura a investigar.

Comparativa ArcMap v/s QGIS 3.28.6

test

En relación a los resultados de proporción de bosques y matorrales para cada punto A y B. Se puede concluir que la “Torre de vigilancia” debiera situarse en el Punto A, ya que contiene mayor Ha de bosque + Matorral. Esto se obtiene aplicando Group Stats primero a la delimitación de la cuenca H con los usos de constitución. Esto para poder conocer primero cual es la totalidad de bosque y matorral dentro de este polígono y asi poder cuantificar la relación de proporcian para punto A y B.

1. Tabla Total Há “SUBUSO” en Cuenca H

SUBUSO/ID_SUB_USO 102 103 108 109
Bosque Mixto 419,08
Bosque Nativo 1123,5
Matorral 74,49
Matorral Arborescente 123,64
TOTAL Bosque+Matorral 1740,7

2. Tabla Total Há “Visibilidad Punto A

SUBUSO/ID_SUB_USO 102 103 108 109
Bosque Mixto 315,91
Bosque Nativo 542,87
Matorral 57,32
Matorral Arborescente 115,91
TOTAL Bosque+Matorral 1032,01

3. Tabla Total Há “Visibilidad Punto B”

SUBUSO/ID_SUB_USO 102 103 108 109
Bosque Mixto 330
Bosque Nativo 954,32
Matorral 53,35
Matorral Arborescente 92,08
TOTAL Bosque+Matorral 422,08

4. Mapa Punto A y su relación con “usos_constitución”

test ### 3.2 Resultados Parte 2 Aplicado el script R se puede observar que la correlación entre variables Indice de Shannon y NDVI promedio es de un valor de correlación Pearson R = 0,29, por tanto no es determinante la relación entre ambas. Esto se calculó por medio de la función “Cor” ( por defecto calcula method = “Pearson”. Mapa indice NDVI se promedio por código de SUBSUSO, de acuerdo a obtención de NDVI de imagen satelital LANDSAT 8 y calculando NDVI de acuerdo a NDVI =(NIR-RED)/(NIR+RED).Con la combinación de bandas respectiva esto es: banda 3 para obtener RED y banda 4 para obtener el NIR. (Infrarojo cercano).

test Se puedo obtener boxplot mediante R, en donde se puede apreciar que para subuso ‘Afloramiento Rocoso’, ‘Minería Industrial’ y ‘Ciudades, Pueblos y Zonas Industriales’ el NDVI promedio es cerccano a 0,4 por tanto se puede inferir que se presenta vegetación moderada en general en estado de salud moderado. Por otra parte la presencia de ‘Plantación’ (114) presenta valores átipicos en gran cantidad. Esto quiere decir que si bien el sensor detecta vegetación robusta para este subuso, también existe presencia de sectores con vegatación escasa, probablemente por la degradación de suelos, producto de la alta explotación por el tipo de plantación. Para ‘Rotación, cultivo y pradera’ altos indices de NDVI pero con outliers incluso negativos. Para el caso de Bosque Nativo y Bosque Mixto se detecta presencia de vegetación robusta y saludable pero con valores átipicos entre 0,4 y 0,6 por tanto se puede inferir que existe presencia de este tipo de vegetación escasa o poco saludable.

Tabla SUBUSO / ID SUBUSO

SUBUSO ID_SUB_USO
Afloramientos Rocosos 101
Bosque Mixto 102
Bosque Nativo 103
Cajas de Ríos 105
Ciudades, Pueblos, Zonas Industriales 106
Lagos, Lagunas, Embalses, Tranques 107
Matorral 108
Matorral Arborescente 109
Matorral-Pradera 110
Minería Industrial 111
Otros Terrenos Húmedos 112
Otros sin Vegetación 113
Plantación 114
Playas y Dunas 115
Praderas 116
Rotación Cultivo-Pradera 117
Ríos 118
Terreno de Uso Agrícola 119
Vegas 120

test

Conclusión Se concluye finalmete que para parte 1 la ubicación de la torre de vigilancia debe estar en Punto A, ya que posee mayor visibilidad dentro de la cuenca H. Dado que tiene una mayor proporcion, esta cercana a 60%. Mientras que el punto B para ubicación de esta torre presenta una proporción de visibilidad cercana al 25 %. Esto tomado del total de Há de Bosque + Matorral, cuyo valor es de 1740,7 Há. Para la parte 2 del laboratorio, se puede concluir que los promedios NDVI y Índice de Shannon si bien tienen una correlación baja de Pearson, en torno a 0,29. Posiblemente tiene una relación que podría intensificarse en ciertos subusos como los Bosques, Bosques Mixtos y sectores arbóreos como Matorral Arborescente. Por tanto aplicando otros métodos estadísticos y de analísis espacial, queda abierta una linea de investigación al respecto. Se observó que los indices responden de forma equivalente para las zonas de baja vegetación. Los valores átipicos dejan en claro ver la actividad forestal en Constitución, ya que la presencia de Plantación y rotación de cultivo probablemente sea relacionado a monocultivo de espcies forestales exógenas.

Referencias

Anexos

Parámetros aplicados en Parte 1

Paso Proceso Parámetros
1 Ordinary Kriging (SAGA GIS) entrada = Elevaciones
atributo = Elevación (metros)
Output Extent = “Marco”
cell size = 50 mts
Error Measure = SD
Logaritmic transformation
maximum search distance = 2000
numbers of point = all points within search distance
minimum = 16
maximum = 25
2 Fill Sinks (Wang Liu) Entrada = MdeFinal
Slope = 0.1
Salida = MDE_Final (No sinks)
3 Channel Network Entrada = MDE_Final (no sinks)
Threshold = 4
Salida = Chanels
4 UpslopeArea (interactivo) Elevation MDE_Final (no sinks)
method = Deterministic 8
Se aplica con cursor punto cercano a coordenada informada. x=742400 y=6082331
Salida = Upslope Area
5 Vectorising Grid classes Vectoriza Chanels y Upslope Area
6 Paso 2 ViewShed pto Torres Viewshed pto A
entrada = MDE FINAL (NO sinks)
nbanda = banda1
Localizacion Observador = ubicación pto A
Altura Obs = 65 m
altura objetivo = 1
Distancia max = 19000 m
Viewshed pto B
entrada = MDE FInal (No sinks)
nbanda = banda 1
Localizacion Observador = Ubicación pto B
Altura Observador = 65 m
altura objetivo = 1
Distancia maxima = 19000 m
7 Aplicar mascara entrada = CuencaVisual
mascara = “CuencaH”
Salida = Cuenca vis
8 Normalización Se normalizan salidas con calculadora raster en valores 0 y 1 con la siguiente expresión:
(“cortadoA” > 0.5) * 1
(“cortadoB” > 0.5) * 1
9 Vectorización Se vectorizan salidas raster intermedias
10 Clip Se aplica proceso de Clip de “usos_constitucion” con feature Cuenca H; con esto se obtiene el área que vamos a trabajar.
11 Corregir topología Se procede a corregir topología con la herramienta “Corregir Geometrías” y se aplican ajustes de Propiedades->Digitalización->Correcciones Automáticas->Allowed Gaps (Vertices Faltantes, Solapa) en propiedades de la capa “GeometríasCorregidas”.
12 Clip Entrada = Geometría Corregida (usos_Constitución)
Salida = HaVisible_PtoA
13 Cálculo de Bosque + Matorral Se realiza el cálculo de Bosque + Matorral con GroupStats, llevando a excel valores de acuerdo a “SUBUSO” y “SUPERF_HA”