Taller de especialización
Herramientas SIG:
Google Earth Engine y R.
Clase 5.
Métricas espaciales.
Autor: Nelson Mattie
Fecha: 17 de febrero de 2021
1 Preparar sesión de trabajo.
1.1 Limpiando las consolas de trabajo
Limpiar la consola de variables
rm(list=ls())Limpiar la consola de trabajo
cat("\014")
1.2 Instalar paquetes.
OBJETIVO: Cargar los paquetes necesarios para desarrollar este práctico Si no tienen los paquetes, deben buscarlos en la consola de abajo a la derecha (opción packages) y bajarlos
if(!require(rgdal)) suppressMessages( suppressWarnings( install.packages("rgdal") )) # para manipulación y visualización de datos geoespaciales"
if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") )) # para manipulación y visualización de datos raster
if(!require(sp)) suppressMessages( suppressWarnings( install.packages("sp") )) # para trabajar con datos vectoriales, debiera cargar junto al paquete raster
if(!require(landscapemetrics)) suppressMessages( suppressWarnings( install.packages("landscapemetrics") )) # Para cálculo de métricas espaciales
if(!require(landscapetools)) suppressMessages( suppressWarnings( install.packages("landscapetools") )) # Herramientas para manipulación y visualización de rasters (esto es opcional)
if(!require(rlang)) suppressMessages( suppressWarnings( install.packages("rlang") )) # Paquetes complementarios para compatibilidad entre códigos y datos en R
if(!require(glue)) suppressMessages( suppressWarnings( install.packages("glue") )) # Paquetes complementarios para compatibilidad entre códigos y datos en RPueden ver el pdf del paquete "landscapemetrics" en: https://cran.r-project.org/web/packages/landscapemetrics/landscapemetrics.pdf
1.3 Definir directorios.
Definir el directorio donde obtendremos y dejaremos los archivos generados
ruta_proyecto <- dirname(getwd())
ruta_clase_01 <- file.path(ruta_proyecto, "Clase_01_GEE")
ruta_clase_02 <- file.path(ruta_proyecto, "Clase_02_R")
ruta_clase_03 <- file.path(ruta_proyecto, "Clase_03_Modelamiento")
ruta_clase_04 <- file.path(ruta_proyecto, "Clase_04_Maxent")
ruta_clase_05 <- file.path(ruta_proyecto, "Clase_05_Metricas")
ruta_clase_06 <- file.path(ruta_proyecto, "Clase_06_amc")
ruta_datos <- file.path(ruta_proyecto, "datos")2 Manipulación de datos.
2.1 Cargar rasters.
Cargamos el raster "Clasificación Final" que generamos la clase pasada con MaxEnt (ojo que la ruta en este caso es ditinta a la que tenemos como directorio de trabajo)
Abierto_MaxEnt <- raster(file.path(ruta_clase_04,"Práctico_Maxent/Resultados_Clasificacion/Abierto_MaxEnt.tif"))
check_landscape(Abierto_MaxEnt) #Para evaluar si el raster cumple con los criterios para poder ser analizado. OJO: el CRS del raster debe estar en metros!## Warning: Caution: Land-cover classes must be decoded as integer values.## Warning: Caution: Coordinate reference system not metric - Units of results based on cellsizes and/or distances may be incorrect.## Warning: Caution: More than 30 land cover-classes - Please check if discrete land-cover classes are present.## layer crs units class n_classes OK
## 1 1 geographic degrees non-integer 991975 ✖### str(Abierto_MaxEnt)Si el raster no cumple con los criterios necesarios, debemos reproyectarlos a coordenadas planas, en este caso EPSG:32719 o UTM_19S
Paisaje_reproject <- projectRaster(Abierto_MaxEnt, crs = "+init=epsg:32719", method="ngb") #Reproyectamos con la opción "nearestneighbor"
check_landscape(Paisaje_reproject) #Con este cambio ahora el resultados debiese ser OK="v"## Warning: Caution: Land-cover classes must be decoded as integer values.## Warning: Caution: More than 30 land cover-classes - Please check if discrete land-cover classes are present.## layer crs units class n_classes OK
## 1 1 projected m non-integer 991498 ✖class(Paisaje_reproject) #Con este cambio ahora el resultados debiese ser OK="v"## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"print(check_landscape(Paisaje_reproject)) #Esto es sólo para ver el resultado en formato tabla## Warning: Caution: Land-cover classes must be decoded as integer values.
## Warning: Caution: More than 30 land cover-classes - Please check if discrete land-cover classes are present.## layer crs units class n_classes OK
## 1 1 projected m non-integer 991498 ✖Observamos dos cosas. Primero es que la reproyección cambia la unidad de grados a metros que es lo esperado. En segundo lugar tenemos que la columna OK nos arroja una "x", esto debido a que el campo clase es "non-integer". Esto es debido al resultado de la clasificación usando MaxEnt puesto que la clase resultante no es un valor entero sino un valor probabilístico entre 0 y 1.
Por esto al final del laboratorio anterior generamos dos imágenes binarias, una usando el estadístico MaxKappa y la otra usando el MaxTSS como puntos de corte para la binarización. Probaremos con alguna de estas imágenes ver si el resultado usando check_landscape se modifica.
Abierto_reclass_maxkappa <- raster(file.path(ruta_clase_04,"Práctico_Maxent/Resultados_Clasificacion/Abierto_reclass_maxkappa.tif"))
check_landscape(Abierto_reclass_maxkappa)## Warning: Caution: Coordinate reference system not metric - Units of results based on cellsizes and/or distances may be incorrect.## layer crs units class n_classes OK
## 1 1 geographic degrees integer 2 ✖El mensaje mejoró porque ahora solamente tenemos la advertencia de que las unidades no son metros, pero no tenemos más advertencias sobre el tipo de datos que es lo que andamos buscando.
Paisaje_reproject_maxkappa <- projectRaster(Abierto_reclass_maxkappa, crs = "+init=epsg:32719", method="ngb")
check_landscape(Paisaje_reproject_maxkappa)## layer crs units class n_classes OK
## 1 1 projected m integer 2 ✔Por fin, logrado el propósito de tener un visto bueno en la columna OK de la función check_landscape, por lo que continuaremos usando esta imagen como insumo. El ejercicio puede repetirse usando la imagen binaria generada con el estadístico MaxTSS.
Solamente mencionamos que intentamos usar la función reclassify para binarizar la imagen, pero no tuvimos éxito, sino que usamos la función calc al final del laboratorio anterior. De todas formas dejamos el código de reclassify con comentarios como evidencia para evitar usarlo en casos futuros.
### matriz_reclass <- matrix(1:9,ncol=3,byrow=TRUE)
### print(matriz_reclass)Nos habría interesado tener una imagen de más categorías que solamente dos que fue parte de lo que intentamos realizar, pero lo descartamos debido a que reclassify no respondió adecuadamente.
### paisaje_reclass <- reclassify( Paisaje_reproject, matriz_reclass , include.lowest=TRUE )
### print(paisaje_reclass)
### plot(paisaje_reclass)
### freq(paisaje_reclass)
### unique(paisaje_reclass)2.2 Calcular las métricas de paisaje.
Primero revisaremos cuáles son las métricas disponibles en el paquete "landcapemetrics" de R
print(head(list_lsm(),10)) #Con esto revisamos todas las métricas que se pueden calcular## # A tibble: 10 x 5
## metric name type level function_name
## <chr> <chr> <chr> <chr> <chr>
## 1 area patch area area and edge metric patch lsm_p_area
## 2 cai core area index core area metric patch lsm_p_cai
## 3 circle related circumscribing circle shape metric patch lsm_p_circle
## 4 contig contiguity index shape metric patch lsm_p_contig
## 5 core core area core area metric patch lsm_p_core
## 6 enn euclidean nearest neighbor distance aggregation metric patch lsm_p_enn
## 7 frac fractal dimension index shape metric patch lsm_p_frac
## 8 gyrate radius of gyration area and edge metric patch lsm_p_gyrate
## 9 ncore number of core areas core area metric patch lsm_p_ncore
## 10 para perimeter-area ratio shape metric patch lsm_p_paraprint(head(list_lsm(level='landscape'),10)) #Si quisiéramos ver sólamente las métricas a nivel de paisaje, podríamos hacer esto## # A tibble: 10 x 5
## metric name type level function_name
## <chr> <chr> <chr> <chr> <chr>
## 1 ai aggregation index aggregation metric landscape lsm_l_ai
## 2 area_cv patch area area and edge metric landscape lsm_l_area_cv
## 3 area_mn patch area area and edge metric landscape lsm_l_area_mn
## 4 area_sd patch area area and edge metric landscape lsm_l_area_sd
## 5 cai_cv core area index core area metric landscape lsm_l_cai_cv
## 6 cai_mn core area index core area metric landscape lsm_l_cai_mn
## 7 cai_sd core area index core area metric landscape lsm_l_cai_sd
## 8 circle_cv related circumscribing circle shape metric landscape lsm_l_circle_cv
## 9 circle_mn related circumscribing circle shape metric landscape lsm_l_circle_mn
## 10 circle_sd related circumscribing circle shape metric landscape lsm_l_circle_sdAhora calcularemos algunas métricas a escala de paisaje
Area_Total_Paisaje<- lsm_l_ta(Paisaje_reproject_maxkappa)
print(Area_Total_Paisaje) #El resultado de esta métrica es en hectáreas## # A tibble: 1 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA ta 79244.Borde_Total_Paisaje<- lsm_l_te(Paisaje_reproject_maxkappa)
print(Borde_Total_Paisaje) #El resultado de esta métrica es en metros## # A tibble: 1 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA te 399987.Número_Parches_Paisaje<- lsm_l_np(Paisaje_reproject_maxkappa)
print(Número_Parches_Paisaje) #El resultado es en número total de parches## # A tibble: 1 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA np 1070Forma_Paisaje<- lsm_l_shape_mn(Paisaje_reproject_maxkappa)
print(Forma_Paisaje) #En este índice 0 es la forma más simple, mientras más alto, la forma es más compleja## # A tibble: 1 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA shape_mn 1.26Diversidad_Shannon_Paisaje<- lsm_l_shdi(Paisaje_reproject_maxkappa)
print(Diversidad_Shannon_Paisaje) #En este índice 0 es diversidad mínima, mientras más alto, más diverso## # A tibble: 1 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA shdi 0.0575Para hacer más simple la visualización, podemos generar una tabla que integre todos estos resultados
Nombres_Metricas_Paisaje<- c("Area_Total_Paisaje", "Borde_Total_Paisaje", "Número_Parches_Paisaje", "Forma_Paisaje", "Diversidad_Shannon_Paisaje")
Tabla_Paisaje<- rbind(Area_Total_Paisaje, Borde_Total_Paisaje, Número_Parches_Paisaje, Forma_Paisaje, Diversidad_Shannon_Paisaje)
print(Tabla_Paisaje)## # A tibble: 5 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA ta 79244.
## 2 1 landscape NA NA te 399987.
## 3 1 landscape NA NA np 1070
## 4 1 landscape NA NA shape_mn 1.26
## 5 1 landscape NA NA shdi 0.0575Métricas_Paisaje<- cbind(Nombres_Metricas_Paisaje, cbind(Tabla_Paisaje[5:6])) #Con esto "adecuamos" nuestra tabla para mejorar su visualización
print(Métricas_Paisaje)## Nombres_Metricas_Paisaje metric value
## 1 Area_Total_Paisaje ta 7.924371e+04
## 2 Borde_Total_Paisaje te 3.999865e+05
## 3 Número_Parches_Paisaje np 1.070000e+03
## 4 Forma_Paisaje shape_mn 1.256767e+00
## 5 Diversidad_Shannon_Paisaje shdi 5.745469e-02Ahora calcularemos algunas métricas a escala de clase
lsm_c_pland(Paisaje_reproject)## # A tibble: 1 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class NA NA pland NAPorcentaje_Paisaje<- lsm_c_pland(Paisaje_reproject_maxkappa)
print(Porcentaje_Paisaje)## # A tibble: 2 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 0 NA pland 99.0
## 2 1 class 1 NA pland 1.03Nuestras clases en esta tabla son 1, 2, 3, 4, 5 y 6. Mejor ponerles nombre no?"
Clases <- c('Abierto', 'Agrícola', 'Agua', 'Construido', 'Nativo', 'Plantación') #Esto es para ponerle nombre a la tablaGeneramos la misma tabla, pero ahora con los nombres que corresponden
Porcentaje_Paisaje_Clase<- cbind(Clases, Porcentaje_Paisaje[,5:6]) #Esto es para filtrar los campos que nos interesean
print(Porcentaje_Paisaje_Clase) #Porcentaje del paisaje por cada clase## Clases metric value
## 1 Abierto pland 98.968267
## 2 Agrícola pland 1.031733
## 3 Agua pland 98.968267
## 4 Construido pland 1.031733
## 5 Nativo pland 98.968267
## 6 Plantación pland 1.031733Calculamos las otras métricas y generamos las tablas con los resultados
Número_Parches<- lsm_c_np(Paisaje_reproject_maxkappa)
Número_Parches_Clase<- cbind(Clases, Número_Parches[,5:6])
print(Número_Parches_Clase) #Número de parches por cada clase## Clases metric value
## 1 Abierto np 32
## 2 Agrícola np 1038
## 3 Agua np 32
## 4 Construido np 1038
## 5 Nativo np 32
## 6 Plantación np 1038Índice_Forma_Paisaje<- lsm_c_lsi(Paisaje_reproject_maxkappa)
Índice_Forma_Paisaje_Clase<- cbind(Clases, Índice_Forma_Paisaje[,5:6])
print(Índice_Forma_Paisaje_Clase) #Índice de forma y agregación de parches. Valores más altos indican formas más complejas y desagregadas## Clases metric value
## 1 Abierto lsi 4.561858
## 2 Agrícola lsi 34.645933
## 3 Agua lsi 4.561858
## 4 Construido lsi 34.645933
## 5 Nativo lsi 4.561858
## 6 Plantación lsi 34.645933Índice_Parche_Mayor<- lsm_c_lpi(Paisaje_reproject_maxkappa)
Índice_Parche_Mayor_Clase<- cbind(Clases, Índice_Parche_Mayor[,5:6])
print(Índice_Parche_Mayor_Clase) #Porcentaje del paisaje cubierto por el parche de mayor tamaño por cada clase## Clases metric value
## 1 Abierto lpi 98.95576540
## 2 Agrícola lpi 0.08031117
## 3 Agua lpi 98.95576540
## 4 Construido lpi 0.08031117
## 5 Nativo lpi 98.95576540
## 6 Plantación lpi 0.08031117Distancia_Euclideana_Promedio<- lsm_c_enn_mn(Paisaje_reproject_maxkappa) #Esta métrica se demora en computar. ¿Por qué será?
Distancia_Euclideana_Promedio_Clase<- cbind(Clases, Distancia_Euclideana_Promedio[,5:6])
print(Distancia_Euclideana_Promedio_Clase) #Distancia promedio (metros) entre parches más cercanos por cada clase## Clases metric value
## 1 Abierto enn_mn 67.83454
## 2 Agrícola enn_mn 194.72602
## 3 Agua enn_mn 67.83454
## 4 Construido enn_mn 194.72602
## 5 Nativo enn_mn 67.83454
## 6 Plantación enn_mn 194.72602Para hacer más simple la visualización, podemos generar una tabla que integre todos estos resultados
Nombres_Métricas_Clases<- c('Porcentaje Paisaje', 'Número Parches', 'Índice Forma', 'Índice Parche Mayor', 'Distancia Euclidiana')
Tabla_Clases<- cbind(Porcentaje_Paisaje_Clase[3], Número_Parches_Clase[3], Índice_Forma_Paisaje_Clase[3], Índice_Parche_Mayor_Clase[3], Distancia_Euclideana_Promedio_Clase[3])
colnames(Tabla_Clases) <- Nombres_Métricas_Clases #Para ponerle nombre a cada columna de la Tabla_Clases
Métricas_Clase<- cbind(Clases, Tabla_Clases)
print(Métricas_Clase)## Clases Porcentaje Paisaje Número Parches Índice Forma Índice Parche Mayor Distancia Euclidiana
## 1 Abierto 98.968267 32 4.561858 98.95576540 67.83454
## 2 Agrícola 1.031733 1038 34.645933 0.08031117 194.72602
## 3 Agua 98.968267 32 4.561858 98.95576540 67.83454
## 4 Construido 1.031733 1038 34.645933 0.08031117 194.72602
## 5 Nativo 98.968267 32 4.561858 98.95576540 67.83454
## 6 Plantación 1.031733 1038 34.645933 0.08031117 194.72602#Si quisieramos guardar la tabla, este es el comando.Ojo "sep=" es para señalar el separador de filas qe utiliza el equipo
#write.table(Métricas_Clase, "./Métricas_Clase.csv", sep=',')2.3 Calcular métricas para distintas partes del paisaje.
Primero debemos "separar" nuestro paisaje original en varios "trozos", para luego analizar cada uno de los trozos
Grilla <- shapefile(file.path(ruta_clase_05,"Cuadrícula_36_UTM19S.shp"), verbose=TRUE) #Cargamos una cuadrícula de muestreo previamente generada en QGIS## OGR data source with driver: ESRI Shapefile
## Source: "/home/matbox/Documents/TrabajosExtra/nelson.mattie/Taller_GEE_R/Clase_05_Metricas", layer: "Cuadrícula_36_UTM19S"
## with 30 features
## It has 5 fields
## Integer64 fields read as strings: idplot(Paisaje_reproject) #para visualizar el mapa cargado
plot(Grilla, border="red", add=TRUE) #añadimos la grilla de muestreo sobre el mapa para visualizarla
Paisaje_1 <- crop(Paisaje_reproject_maxkappa, Grilla[1,]) #Ahora cortamos nuestro raster al tamaño del primero polígono
plot(Paisaje_1, main = "Paisaje 1") #Con este comando vemos el resultado de nuestro nuevo "pedazo" de paisaje
Paisaje_2 <- crop(Paisaje_reproject_maxkappa, Grilla[2,]) #Ahora cortamos nuestro raster al tamaño del segundo polígono
plot(Paisaje_2, main = "Paisaje 2") #Con este comando vemos el resultado de nuestro nuevo "pedazo" de paisaje
Paisaje_3 <- crop(Paisaje_reproject_maxkappa, Grilla[3,]) #Hagamos el mismo proceso con el resto de los polígonos
Paisaje_4 <- crop(Paisaje_reproject_maxkappa, Grilla[4,])
Paisaje_5 <- crop(Paisaje_reproject_maxkappa, Grilla[5,])
Paisaje_6 <- crop(Paisaje_reproject_maxkappa, Grilla[6,])
Paisaje_7 <- crop(Paisaje_reproject_maxkappa, Grilla[7,])
Paisaje_8 <- crop(Paisaje_reproject_maxkappa, Grilla[8,])
Paisaje_9 <- crop(Paisaje_reproject_maxkappa, Grilla[9,])
Paisaje_10 <- crop(Paisaje_reproject_maxkappa, Grilla[10,])
Paisaje_11 <- crop(Paisaje_reproject_maxkappa, Grilla[11,])
Paisaje_12 <- crop(Paisaje_reproject_maxkappa, Grilla[12,])
Paisaje_13 <- crop(Paisaje_reproject_maxkappa, Grilla[13,])
Paisaje_14 <- crop(Paisaje_reproject_maxkappa, Grilla[14,])
Paisaje_15 <- crop(Paisaje_reproject_maxkappa, Grilla[15,])
Paisaje_16 <- crop(Paisaje_reproject_maxkappa, Grilla[16,])
Paisaje_17 <- crop(Paisaje_reproject_maxkappa, Grilla[17,])
Paisaje_18 <- crop(Paisaje_reproject_maxkappa, Grilla[18,])
Paisaje_19 <- crop(Paisaje_reproject_maxkappa, Grilla[19,])
Paisaje_20 <- crop(Paisaje_reproject_maxkappa, Grilla[20,])
Paisaje_21 <- crop(Paisaje_reproject_maxkappa, Grilla[21,])
Paisaje_22 <- crop(Paisaje_reproject_maxkappa, Grilla[22,])
Paisaje_23 <- crop(Paisaje_reproject_maxkappa, Grilla[23,])
Paisaje_24 <- crop(Paisaje_reproject_maxkappa, Grilla[24,])
Paisaje_25 <- crop(Paisaje_reproject_maxkappa, Grilla[25,])
Paisaje_26 <- crop(Paisaje_reproject_maxkappa, Grilla[26,])
Paisaje_27 <- crop(Paisaje_reproject_maxkappa, Grilla[27,])
Paisaje_28 <- crop(Paisaje_reproject_maxkappa, Grilla[28,])
Paisaje_29 <- crop(Paisaje_reproject_maxkappa, Grilla[29,])
Paisaje_30 <- crop(Paisaje_reproject_maxkappa, Grilla[30,])Ahora calculemos algunas métricas espaciales a NIVEL DE PAISAJE para los distintos paisajes creados Diversidad de Shannon por Paisaje (DSP)"
DSP_Paisaje_1 <- lsm_l_shdi(Paisaje_1)
DSP_Paisaje_2 <- lsm_l_shdi(Paisaje_2)
DSP_Paisaje_3 <- lsm_l_shdi(Paisaje_3)
DSP_Paisaje_4 <- lsm_l_shdi(Paisaje_4)
DSP_Paisaje_5 <- lsm_l_shdi(Paisaje_5)
DSP_Paisaje_6 <- lsm_l_shdi(Paisaje_6)
DSP_Paisaje_7 <- lsm_l_shdi(Paisaje_7)
DSP_Paisaje_8 <- lsm_l_shdi(Paisaje_8)
DSP_Paisaje_9 <- lsm_l_shdi(Paisaje_9)
DSP_Paisaje_10 <- lsm_l_shdi(Paisaje_10)
DSP_Paisaje_11 <- lsm_l_shdi(Paisaje_11)
DSP_Paisaje_12 <- lsm_l_shdi(Paisaje_12)
DSP_Paisaje_13 <- lsm_l_shdi(Paisaje_13)
DSP_Paisaje_14 <- lsm_l_shdi(Paisaje_14)
DSP_Paisaje_15 <- lsm_l_shdi(Paisaje_15)
DSP_Paisaje_16 <- lsm_l_shdi(Paisaje_16)
DSP_Paisaje_17 <- lsm_l_shdi(Paisaje_17)
DSP_Paisaje_18 <- lsm_l_shdi(Paisaje_18)
DSP_Paisaje_19 <- lsm_l_shdi(Paisaje_19)
DSP_Paisaje_20 <- lsm_l_shdi(Paisaje_20)
DSP_Paisaje_21 <- lsm_l_shdi(Paisaje_21)
DSP_Paisaje_22 <- lsm_l_shdi(Paisaje_22)
DSP_Paisaje_23 <- lsm_l_shdi(Paisaje_23)
DSP_Paisaje_24 <- lsm_l_shdi(Paisaje_24)
DSP_Paisaje_25 <- lsm_l_shdi(Paisaje_25)
DSP_Paisaje_26 <- lsm_l_shdi(Paisaje_26)
DSP_Paisaje_27 <- lsm_l_shdi(Paisaje_27)
DSP_Paisaje_28 <- lsm_l_shdi(Paisaje_28)
DSP_Paisaje_29 <- lsm_l_shdi(Paisaje_29)
DSP_Paisaje_30 <- lsm_l_shdi(Paisaje_30)Veamos un ejemplo de los resultados para el Paisaje_6
plot(Paisaje_6, main = "Paisaje 6")
print(DSP_Paisaje_6)## # A tibble: 1 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA shdi 0.0118Ahora generamos una tabla para ver los resultados integrados, pero sólo de la fila de datos que nos interesa
DSP_Integrado<- rbind(DSP_Paisaje_1[,6], DSP_Paisaje_2[,6], DSP_Paisaje_3[,6], DSP_Paisaje_4[,6], DSP_Paisaje_5[,6],
DSP_Paisaje_6[,6], DSP_Paisaje_7[,6], DSP_Paisaje_8[,6], DSP_Paisaje_9[,6], DSP_Paisaje_10[,6],
DSP_Paisaje_11[,6], DSP_Paisaje_12[,6], DSP_Paisaje_13[,6], DSP_Paisaje_14[,6], DSP_Paisaje_15[,6],
DSP_Paisaje_16[,6], DSP_Paisaje_17[,6], DSP_Paisaje_18[,6], DSP_Paisaje_19[,6], DSP_Paisaje_20[,6],
DSP_Paisaje_21[,6], DSP_Paisaje_22[,6], DSP_Paisaje_23[,6], DSP_Paisaje_24[,6], DSP_Paisaje_25[,6],
DSP_Paisaje_26[,6], DSP_Paisaje_27[,6], DSP_Paisaje_28[,6], DSP_Paisaje_29[,6], DSP_Paisaje_30[,6])Veamos la tabla generada
colnames(DSP_Integrado) <- "DSP" #Esto nos sirve para ponerle un nombre ala columna
print(DSP_Integrado)## # A tibble: 30 x 1
## DSP
## <dbl>
## 1 0.00351
## 2 0.0308
## 3 0.0315
## 4 0.0212
## 5 0.00759
## 6 0.0118
## 7 0.101
## 8 0.0760
## 9 0.0685
## 10 0.00673
## 11 0.00389
## 12 0.0629
## 13 0.0691
## 14 0.121
## 15 0.0245
## 16 0.0466
## 17 0.0140
## 18 0.0928
## 19 0.0675
## 20 0.0412
## 21 0.0855
## 22 0.0260
## 23 0.108
## 24 0.201
## 25 0.0169
## 26 0.0131
## 27 0.0964
## 28 0.0744
## 29 0.103
## 30 0.0177Ahora realizaremos un simple gráfico para ver cómo nuestra variable varía entre los diferentes paisajes
Eje_X<- c(1:30) #Con este comando generamos una nueva columna de 25 números que usaremos como eje X en nuestro gráfico
Tabla <- cbind(Eje_X,DSP_Integrado) #Integramos los datos y graficamos
plot(Tabla$Eje_X,Tabla$DSP, type = "b", xlab = "Paisaje", ylab = "N° Parches", col = 'blue')
Y cómo se ve nuestra variable DSP si la mapeamos? Para ello, primero vamos a hacer un "merge" de los datos entre nuestra Tabla y la Grilla que usamos
names(Tabla) #Con esto podemos ver el nombre de las variables de nuestra tabla## [1] "Eje_X" "DSP"names(Grilla) #Lo mismo para ver el nombre de las columnas de la tabla## [1] "left" "top" "right" "bottom" "id"Claramente no tenemos ninguna columna que se repita en ambas bases de datos...¿qué podemos hacer?
colnames(Tabla) <- c("id","DSP") #La variable Eje_x es equivalente al id del polígono, por lo que podemos cambiarle el nombre a "id"
Grilla_DSP<- merge(Grilla,Tabla, by="id") #Ahora podemos unir nuestro shape (Grilla) con nuestra tabla, usando "id" como columna común.
spplot(Grilla_DSP["DSP"], main = "DSP") #Ahora ya podemos maperar nuestro resultado de DSP". Ojo que usamos la función spplot() en vez de plot().
Guardemos nuestro resultado, ya que nos servirá para la próxima clase
if(!file.exists(file.path(ruta_clase_05,"Resultados")))
{ dir.create(file.path(ruta_clase_05,"Resultados")) }
writeOGR(Grilla_DSP, dsn = file.path(ruta_clase_05,"Resultados"),
layer = "Shannon_Paisaje", driver = "ESRI Shapefile",
overwrite=TRUE)Ahora calculemos algunas métricas espaciales a NIVEL DE CLASE para los distintos paisajes creados Porcentaje de Paisaje por clase (PPC)"
PPC_Paisaje_1 <- lsm_c_pland(Paisaje_1)
PPC_Paisaje_2 <- lsm_c_pland(Paisaje_2)
PPC_Paisaje_3 <- lsm_c_pland(Paisaje_3)
PPC_Paisaje_4 <- lsm_c_pland(Paisaje_4)
PPC_Paisaje_5 <- lsm_c_pland(Paisaje_5)
PPC_Paisaje_6 <- lsm_c_pland(Paisaje_6)
PPC_Paisaje_7 <- lsm_c_pland(Paisaje_7)
PPC_Paisaje_8 <- lsm_c_pland(Paisaje_8)
PPC_Paisaje_9 <- lsm_c_pland(Paisaje_9)
PPC_Paisaje_10 <- lsm_c_pland(Paisaje_10)
PPC_Paisaje_11 <- lsm_c_pland(Paisaje_11)
PPC_Paisaje_12 <- lsm_c_pland(Paisaje_12)
PPC_Paisaje_13 <- lsm_c_pland(Paisaje_13)
PPC_Paisaje_14 <- lsm_c_pland(Paisaje_14)
PPC_Paisaje_15 <- lsm_c_pland(Paisaje_15)
PPC_Paisaje_16 <- lsm_c_pland(Paisaje_16)
PPC_Paisaje_17 <- lsm_c_pland(Paisaje_17)
PPC_Paisaje_18 <- lsm_c_pland(Paisaje_18)
PPC_Paisaje_19 <- lsm_c_pland(Paisaje_19)
PPC_Paisaje_20 <- lsm_c_pland(Paisaje_20)
PPC_Paisaje_21 <- lsm_c_pland(Paisaje_21)
PPC_Paisaje_22 <- lsm_c_pland(Paisaje_22)
PPC_Paisaje_23 <- lsm_c_pland(Paisaje_23)
PPC_Paisaje_24 <- lsm_c_pland(Paisaje_24)
PPC_Paisaje_25 <- lsm_c_pland(Paisaje_25)
PPC_Paisaje_26 <- lsm_c_pland(Paisaje_26)
PPC_Paisaje_27 <- lsm_c_pland(Paisaje_27)
PPC_Paisaje_28 <- lsm_c_pland(Paisaje_28)
PPC_Paisaje_29 <- lsm_c_pland(Paisaje_29)
PPC_Paisaje_30 <- lsm_c_pland(Paisaje_30)Veamos un ejemplo de los resultados para el Paisaje_10
print(PPC_Paisaje_10) #¿Cuántas clases hay en este paisaje?## # A tibble: 2 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 0 NA pland 99.9
## 2 1 class 1 NA pland 0.0831Para asegurarnos de no tener problemas con el número de clases, primero usaremos la función "merge Esta función nos permitirá asegurar que todas nuestros resultados tengan el mismo número de filas
Clases_ok<- data.frame(1:6) #Con esto creamos una tabla de una columna que representa las 6 clases
colnames(Clases_ok) <- "class"PPC_Paisaje_1 <- merge(Clases_ok, PPC_Paisaje_1, by="class", all=TRUE)
PPC_Paisaje_2 <- merge(Clases_ok, PPC_Paisaje_2, by="class", all=TRUE)
PPC_Paisaje_3 <- merge(Clases_ok, PPC_Paisaje_3, by="class", all=TRUE)
PPC_Paisaje_4 <- merge(Clases_ok, PPC_Paisaje_4, by="class", all=TRUE)
PPC_Paisaje_5 <- merge(Clases_ok, PPC_Paisaje_5, by="class", all=TRUE)
PPC_Paisaje_6 <- merge(Clases_ok, PPC_Paisaje_6, by="class", all=TRUE)
PPC_Paisaje_7 <- merge(Clases_ok, PPC_Paisaje_7, by="class", all=TRUE)
PPC_Paisaje_8 <- merge(Clases_ok, PPC_Paisaje_8, by="class", all=TRUE)
PPC_Paisaje_9 <- merge(Clases_ok, PPC_Paisaje_9, by="class", all=TRUE)
PPC_Paisaje_10 <- merge(Clases_ok, PPC_Paisaje_10, by="class", all=TRUE)
PPC_Paisaje_11 <- merge(Clases_ok, PPC_Paisaje_11, by="class", all=TRUE)
PPC_Paisaje_12 <- merge(Clases_ok, PPC_Paisaje_12, by="class", all=TRUE)
PPC_Paisaje_13 <- merge(Clases_ok, PPC_Paisaje_13, by="class", all=TRUE)
PPC_Paisaje_14 <- merge(Clases_ok, PPC_Paisaje_14, by="class", all=TRUE)
PPC_Paisaje_15 <- merge(Clases_ok, PPC_Paisaje_15, by="class", all=TRUE)
PPC_Paisaje_16 <- merge(Clases_ok, PPC_Paisaje_16, by="class", all=TRUE)
PPC_Paisaje_17 <- merge(Clases_ok, PPC_Paisaje_17, by="class", all=TRUE)
PPC_Paisaje_18 <- merge(Clases_ok, PPC_Paisaje_18, by="class", all=TRUE)
PPC_Paisaje_19 <- merge(Clases_ok, PPC_Paisaje_19, by="class", all=TRUE)
PPC_Paisaje_20 <- merge(Clases_ok, PPC_Paisaje_20, by="class", all=TRUE)
PPC_Paisaje_21 <- merge(Clases_ok, PPC_Paisaje_21, by="class", all=TRUE)
PPC_Paisaje_22 <- merge(Clases_ok, PPC_Paisaje_22, by="class", all=TRUE)
PPC_Paisaje_23 <- merge(Clases_ok, PPC_Paisaje_23, by="class", all=TRUE)
PPC_Paisaje_24 <- merge(Clases_ok, PPC_Paisaje_24, by="class", all=TRUE)
PPC_Paisaje_25 <- merge(Clases_ok, PPC_Paisaje_25, by="class", all=TRUE)
PPC_Paisaje_26 <- merge(Clases_ok, PPC_Paisaje_26, by="class", all=TRUE)
PPC_Paisaje_27 <- merge(Clases_ok, PPC_Paisaje_27, by="class", all=TRUE)
PPC_Paisaje_28 <- merge(Clases_ok, PPC_Paisaje_28, by="class", all=TRUE)
PPC_Paisaje_29 <- merge(Clases_ok, PPC_Paisaje_29, by="class", all=TRUE)
PPC_Paisaje_30 <- merge(Clases_ok, PPC_Paisaje_30, by="class", all=TRUE)Ahora integraremos nuestros resultados en una sola tabla, recogiendo sólo los datos de nuestro interés
PPC_Integrado <- cbind(PPC_Paisaje_1[,6], PPC_Paisaje_2[,6], PPC_Paisaje_3[,6], PPC_Paisaje_4[,6], PPC_Paisaje_5[,6],
PPC_Paisaje_6[,6], PPC_Paisaje_7[,6], PPC_Paisaje_8[,6], PPC_Paisaje_9[,6], PPC_Paisaje_10[,6],
PPC_Paisaje_11[,6], PPC_Paisaje_12[,6], PPC_Paisaje_13[,6], PPC_Paisaje_14[,6], PPC_Paisaje_15[,6],
PPC_Paisaje_16[,6], PPC_Paisaje_17[,6], PPC_Paisaje_18[,6], PPC_Paisaje_19[,6], PPC_Paisaje_20[,6],
PPC_Paisaje_21[,6], PPC_Paisaje_22[,6], PPC_Paisaje_23[,6], PPC_Paisaje_24[,6], PPC_Paisaje_25[,6],
PPC_Paisaje_26[,6], PPC_Paisaje_27[,6], PPC_Paisaje_28[,6], PPC_Paisaje_29[,6], PPC_Paisaje_30[,6])Estos serán los nombres de nuestras columnas
Columnas <- c("PPC1" , "PPC2" , "PPC3" , "PPC4" , "PPC5" , "PPC6" , "PPC7" , "PPC8" , "PPC9" , "PPC10",
"PPC11", "PPC12", "PPC13", "PPC14", "PPC15", "PPC16", "PPC17", "PPC18", "PPC19", "PPC20",
"PPC21", "PPC22", "PPC23", "PPC24", "PPC25", "PPC26", "PPC27", "PPC28", "PPC29", "PPC30")Ahora le asignamos el nombre a las filas y columnas y vemos el resultado
PPC_Integrado <- PPC_Integrado[1:6,]
rownames(PPC_Integrado) <- Clases
colnames(PPC_Integrado) <- Columnas
print(PPC_Integrado)## PPC1 PPC2 PPC3 PPC4 PPC5 PPC6 PPC7 PPC8 PPC9 PPC10 PPC11 PPC12 PPC13 PPC14 PPC15 PPC16 PPC17 PPC18 PPC19 PPC20 PPC21 PPC22 PPC23 PPC24 PPC25 PPC26
## Abierto 99.96024647 99.5124548 99.5000843 99.6873856 99.9045212 99.8410707 97.92515 98.545054 98.720186 99.9168952 99.95532349 98.847614 98.706335 97.387739 99.6287986 99.1999105 99.8060888 98.132912 98.743455 99.3093846 98.31492 99.601097 97.731239 94.902906 99.7589961 99.8210191
## Agrícola 0.03975353 0.4875452 0.4999157 0.3126144 0.0954788 0.1589293 2.07485 1.454946 1.279814 0.0831048 0.04467651 1.152386 1.293665 2.612261 0.3712014 0.8000895 0.1939112 1.867088 1.256545 0.6906154 1.68508 0.398903 2.268761 5.097094 0.2410039 0.1789809
## Agua NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## Construido NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## Nativo NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## Plantación NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
## PPC27 PPC28 PPC29 PPC30
## Abierto 98.043113 98.582568 97.880813 99.7465887
## Agrícola 1.956887 1.417432 2.119187 0.2534113
## Agua NA NA NA NA
## Construido NA NA NA NA
## Nativo NA NA NA NA
## Plantación NA NA NA NAFinalmente desarrollaremos algunos análisis exploratorios de nuestros resultados
Como pudimos ver con el comando "View(PPC_INtegrado)", nuestros datos están al revés, por lo que tenemos que transponerlos
Resultados_Clase <- as.data.frame(t(PPC_Integrado)) #Con esto transponemos nuestros datos y forzamos a tratarlos como tabla de datos
print(Resultados_Clase)## Abierto Agrícola Agua Construido Nativo Plantación
## PPC1 99.96025 0.03975353 NA NA NA NA
## PPC2 99.51245 0.48754518 NA NA NA NA
## PPC3 99.50008 0.49991574 NA NA NA NA
## PPC4 99.68739 0.31261441 NA NA NA NA
## PPC5 99.90452 0.09547880 NA NA NA NA
## PPC6 99.84107 0.15892932 NA NA NA NA
## PPC7 97.92515 2.07484972 NA NA NA NA
## PPC8 98.54505 1.45494627 NA NA NA NA
## PPC9 98.72019 1.27981385 NA NA NA NA
## PPC10 99.91690 0.08310480 NA NA NA NA
## PPC11 99.95532 0.04467651 NA NA NA NA
## PPC12 98.84761 1.15238649 NA NA NA NA
## PPC13 98.70634 1.29366464 NA NA NA NA
## PPC14 97.38774 2.61226073 NA NA NA NA
## PPC15 99.62880 0.37120142 NA NA NA NA
## PPC16 99.19991 0.80008952 NA NA NA NA
## PPC17 99.80609 0.19391119 NA NA NA NA
## PPC18 98.13291 1.86708773 NA NA NA NA
## PPC19 98.74346 1.25654450 NA NA NA NA
## PPC20 99.30938 0.69061543 NA NA NA NA
## PPC21 98.31492 1.68507979 NA NA NA NA
## PPC22 99.60110 0.39890302 NA NA NA NA
## PPC23 97.73124 2.26876091 NA NA NA NA
## PPC24 94.90291 5.09709410 NA NA NA NA
## PPC25 99.75900 0.24100391 NA NA NA NA
## PPC26 99.82102 0.17898093 NA NA NA NA
## PPC27 98.04311 1.95688690 NA NA NA NA
## PPC28 98.58257 1.41743247 NA NA NA NA
## PPC29 97.88081 2.11918686 NA NA NA NA
## PPC30 99.74659 0.25341131 NA NA NA NAcor(Resultados_Clase) #Esto genera una matriz de correlación entre nuestras variables## Abierto Agrícola Agua Construido Nativo Plantación
## Abierto 1 -1 NA NA NA NA
## Agrícola -1 1 NA NA NA NA
## Agua NA NA 1 NA NA NA
## Construido NA NA NA 1 NA NA
## Nativo NA NA NA NA 1 NA
## Plantación NA NA NA NA NA 1### cor(Resultados_Clase, use="complete.obs") #Esto genera una matriz de correlación entre nuestras variablesVeamos si estas relaciones son estadísticamente significativas
### cor.test(Resultados_Clase$Abierto, Resultados_Clase$Nativo) #Test paramétrico
### cor.test(Resultados_Clase$Abierto, Resultados_Clase$Nativo, method = "spearman") #Test no paramétrico
### cor.test(Resultados_Clase$Abierto, Resultados_Clase$Plantación)
### cor.test(Resultados_Clase$Abierto, Resultados_Clase$Plantación, method = "spearman")Finalmente, traslademos algunos de los resultados a un mapa para ver su distribución
Tabla_Final <- cbind(Tabla, Resultados_Clase)
Grilla_Final <- merge(Grilla,Tabla_Final, by="id")
spplot(Grilla_Final["Abierto"], main = "% Cobertura Abierto")
spplot(Grilla_Final["Agrícola"], main = "% Cobertura Agrícola")
### spplot(Grilla_Final["Agua"], main = "% Cobertura Agua")
### spplot(Grilla_Final["Construido"], main = "% Cobertura Construido")
### spplot(Grilla_Final["Nativo"], main = "% Cobertura Nativo")
### spplot(Grilla_Final["Plantación"], main = "% Cobertura Plantación")Guardemos nuestro resultado, ya que nos servirá para la próxima clase
writeOGR(Grilla_Final, dsn = file.path(ruta_clase_05,"Resultados"),
layer = "Proporción_Paisaje_Clase", driver = "ESRI Shapefile",
overwrite=TRUE)