GeografĂa Ambiental : Ejercicio Práctico en R: Paisajes
Dra. Monica Ortiz, monicaortiz@udec.cl
IntroducciĂłn
¡Hola! Este sencillo script ayudará a echar un vistazo al análisis espacial de paisajes con R.
rm(list=ls()) # Borrar todo
## Instalar los paquetes necesarios
# install.packages("terra") #análisis raster
# install.packages("sf") #análisis vector
# install.packages("landscapemetrics") #análisis paisaje
# install.packages("landscapetools")
## Cargarlos
library(terra)## Warning: package 'terra' was built under R version 4.3.3## terra 1.7.78library(landscapemetrics)## Warning: package 'landscapemetrics' was built under R version 4.3.3library(sf)## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE## Especificar el directorio de la carpeta y los archivos
setwd("~/Desktop/practica")
## Crear un raster de 5x5
r <- rast(nrows=5, ncols=5)
## Asignar valores que representen diferentes tipos de cobertura del suelo
values(r) <- c(1, 1, 2, 2, 3,
1, 1, 2, 3, 3,
3, 3, 2, 2, 2,
1, 1, 2, 2, 2,
2, 1, 1, 2, 3)
## Graficar el raster
plot(r, main="Tipos de Cobertura del Suelo, 2000")
cover <- c("plantaciĂłn", "bosque", "agricultura")
plot(r, main="Tipos de Cobertura del Suelo, 2000", type="classes", levels=cover)
## Chequear el paisaje, cuantos tipos de clases hay
check_landscape(r, verbose=FALSE)## layer crs units class n_classes OK
## 1 1 geographic degrees integer 3 ✖## Calcular el área de cada uno de los parches
lsm_p_area(r)## # A tibble: 7 Ă— 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 patch 1 1 area 1.04
## 2 1 patch 1 2 area 1.04
## 3 1 patch 2 3 area 0.259
## 4 1 patch 2 4 area 2.59
## 5 1 patch 3 5 area 0.518
## 6 1 patch 3 6 area 0.778
## 7 1 patch 3 7 area 0.259## Y su area
rc<-lsm_c_ca(r)
## Area total del paisaje
lsm_l_ta(r)## # A tibble: 1 Ă— 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA ta 6.48##### Pausa: Discutir los resultados y su significado#####
## Crear otro raster de 5x5, para el año 2020
r1 <- rast(nrows=5, ncols=5)
values(r1) <- c(1, 1, 1, 1, 3,
1, 1, 2, 3, 3,
3, 3, 2, 2, 2,
1, 1, 1, 1, 2,
1, 1, 1, 2, 2)
par(mfrow=c(2,1))
plot(r, main="Tipos de Cobertura del Suelo, 2000", type="classes", levels=cover)
plot(r1, main="Tipos de Cobertura del Suelo, 2020", type="classes", levels=cover)
## Cuantos tipos de clases hay
check_landscape(r1, verbose=FALSE)## layer crs units class n_classes OK
## 1 1 geographic degrees integer 3 ✖## Calcular el área de cada uno de los parches
lsm_p_area(r1)## # A tibble: 5 Ă— 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 patch 1 1 area 1.56
## 2 1 patch 1 2 area 1.81
## 3 1 patch 2 3 area 1.81
## 4 1 patch 3 4 area 0.518
## 5 1 patch 3 5 area 0.778## Y su area
r1c<-lsm_c_ca(r1)
## Area total del paisaje
lsm_l_ta(r1)## # A tibble: 1 Ă— 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA ta 6.48## Resumir el análisis
resumen<-data.frame(cover, rc$value, r1c$value)
names(resumen)<-c("Land Cover", "Year_2000", "Year_2020")
resumen$porcentaje<- abs(resumen$Year_2020 - resumen$Year_2000)/(resumen$Year_2000)* 100
## Muestra los resultados/ resumen
resumen## Land Cover Year_2000 Year_2020 porcentaje
## 1 plantaciĂłn 2.0736 3.3696 62.50000
## 2 bosque 2.8512 1.8144 36.36364
## 3 agricultura 1.5552 1.2960 16.66667##### Pausa: Discutir los resultados y su significado#####
## Leer ráster de uso del suelo
lc<-rast("reproj_lc_2014.tif")
## Recalibrar el plot y visualizar el mapa
par(mfrow=c(1,1))
## Leer lĂmites administrativos
chile_shp<-vect("chl_adm_shp/chl_admbnda_adm1_bcn_20211008.shp")
head(chile_shp) #ver las primeras lineas del vector## Shape_Leng Shape_Area ADM1_ES ADM1_PCODE
## 1 10.716495 3.657773 Región de Tarapacá CL01
## 2 21.600173 11.156584 RegiĂłn de Antofagasta CL02
## 3 20.211632 6.907941 RegiĂłn de Atacama CL03
## 4 17.168212 3.817481 RegiĂłn de Coquimbo CL04
## 5 13.951085 1.562133 RegiĂłn de ValparaĂso CL05
## 6 8.057616 1.603334 RegiĂłn del Libertador Bernardo O'Higgins CL06
## ADM1_REF ADM1ALT1ES ADM1ALT2ES ADM0_ES
## 1 Region de Tarapaca <NA> <NA> Chile
## 2 Region de Antofagasta <NA> <NA> Chile
## 3 Region de Atacama <NA> <NA> Chile
## 4 Region de Coquimbo <NA> <NA> Chile
## 5 Region de Valparaiso <NA> <NA> Chile
## 6 Region del Libertador Bernardo O Higgins <NA> <NA> Chile
## ADM0_PCODE date validOn validTo
## 1 CL 2021/08/03 2021/10/08 <NA>
## 2 CL 2021/08/03 2021/10/08 <NA>
## 3 CL 2021/08/03 2021/10/08 <NA>
## 4 CL 2021/08/03 2021/10/08 <NA>
## 5 CL 2021/08/03 2021/10/08 <NA>
## 6 CL 2021/08/03 2021/10/08 <NA>## Subset La RegiĂłn del BiobĂo
biobio <- chile_shp[chile_shp$ADM1_ES == "RegiĂłn del BĂo-BĂo", ]
## Crop, mask y reproyectar
lc_biobio_crop<- crop(lc, biobio)
lc_biobio_mask<- mask(lc_biobio_crop,biobio)
lc_biobio_mask<-project(lc_biobio_mask, "EPSG:32719")## |---------|---------|---------|---------|========================================= ## Visualizar
plot(lc_biobio_mask)
##### Pausa: Discutir diferencia entre Crop y Mask#####
##### Pausa: Este es un paisaje? #####
## Cuenca
cuencas<-vect("SubsubcuencasBNA/Subsubcuencas_BNA.shp") #leer shape
head(cuencas) #visualizar primeras lineas## COD_CUEN COD_SUBC COD_SSUBC
## 1 010 0100 01000
## 2 012 0120 01200
## 3 012 0120 01201
## 4 012 0120 01202
## 5 010 0100 01001
## 6 011 0111 01110
## NOM_SSUBC Shape_Leng Shape_Area
## 1 Rios Uchusuma, Colpas, Putani y Cosapilla 147639.53 817235662
## 2 Rio Lluta bajo Rio Azufre 91675.53 420876494
## 3 Rio Lluta entre Rio Azufre y bajo Quebrada Huaylas 151160.51 991683107
## 4 Rio Lluta entre Quebrada Huaylas y Quebrada de Socoroma 146439.84 892613571
## 5 Rio Caquena hasta frontera (Rio Cosapilla) 112079.22 546015713
## 6 Quebradas Escritos y de La Concordia 182110.12 784799173## Subset Cuenca
# cuencas$NOM_SSUBC #como se llaman los subsubcuencas?
cuenca_biobio <- cuencas[cuencas$NOM_SSUBC == "Rio Bio-Bio Entre Estero Hualqui y Desembocadura", ] #subset
plot(cuenca_biobio) 
crs(cuenca_biobio) # ¿qué es el sistema de referencia de coordenadas?## [1] "PROJCRS[\"WGS 84 / UTM zone 19S\",\n BASEGEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 19S\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"Degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-69,\n ANGLEUNIT[\"Degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",10000000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",32719]]"## Crop + Mask otra vez
r2<- crop(lc_biobio_mask, cuenca_biobio) #crop
r2<-mask(r2, cuenca_biobio) #mask
## Plots
plot(r2)
plot(lc_biobio_mask, add=T)
plot(cuenca_biobio, add=T) #delimitaciĂłn administrativo
## Averigüe cuántas clases hay
unique_values <- unique(values(r2))
num_classes <- length(unique_values)
num_classes #mucho mucho!! (por la proyección también)## [1] 242850## Para que sea más pequeño, vamos a modificar la extensión:
## ExtensiĂłn (xmin, xmax, ymin, ymax)
extent_to_crop <- ext(123000, 125000, 5920000, 5924000) #Caleta Chome / Lenga
## Crop
r2_small <- crop(r2, extent_to_crop)
plot(r2_small)
check_landscape(r2_small)## 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 2300 âś–## Reclasificando
## Defina la matriz de reclasificaciĂłn
## Valores originales -> Valores nuevos
## Ejemplo: Grupo 1-2 -> 1, Grupo 3-4 -> 2, Grupo 5 -> 3
## Referencia original: Zhao et al 2016
reclass_matrix <- matrix(c(0, 99, 0, #los valores entre 0-90 --> 0
100, 199, 100,
200, 245, 200,
245, 299, 250,
300, 399, 300,
400, 499, 400,
500, 599, 500,
600, 699, 600,
700, 899, 800,
900, 999, 900,
1000, 1500, 1000),
ncol=3, byrow=TRUE)
## ReclasificaciĂłn
r2_small_r <- round(r2_small) #redondear valores
r2_small_r <- classify(r2_small_r, reclass_matrix, include.lowest=TRUE,brackets=TRUE)
unique_values <- unique(values(r2_small_r))
num_classes <- length(unique_values)
num_classes #ok, mejor!## [1] 10## Check landscape
check_landscape(r2_small_r) #ok! yahoo.## layer crs units class n_classes OK
## 1 1 projected m integer 9 âś”## Visualizar los parches
show_patches(r2_small_r)## $layer_1
## Ver la lista completa de los metricos que podemos calcular
## print(list_lsm(),n=133)
## P. ej. Calcular el área de cada uno de los parches
lsm_p_area(r2_small_r)## # A tibble: 266 Ă— 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 patch 0 1 area 0.414
## 2 1 patch 0 2 area 0.138
## 3 1 patch 0 3 area 0.621
## 4 1 patch 0 4 area 0.276
## 5 1 patch 0 5 area 0.0690
## 6 1 patch 0 6 area 0.207
## 7 1 patch 0 7 area 0.276
## 8 1 patch 0 8 area 0.345
## 9 1 patch 0 9 area 0.138
## 10 1 patch 0 10 area 0.621
## # â„ą 256 more rows## P. ej. Su area (cada clase)
lsm_c_ca(r2_small_r)## # A tibble: 9 Ă— 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 0 NA ca 5.86
## 2 1 class 100 NA ca 1.72
## 3 1 class 200 NA ca 63.1
## 4 1 class 250 NA ca 49.9
## 5 1 class 300 NA ca 97.5
## 6 1 class 400 NA ca 11.4
## 7 1 class 500 NA ca 0.759
## 8 1 class 600 NA ca 0.276
## 9 1 class 800 NA ca 0.414## Plot
## Land cover classes (Zhao et al, 2016), simplificado
cover2 <- c("sin data", "agricultura", "bosque", "plantaciĂłn", "pradera", "arbustos",
"humedales", "agua", "impermeable")
plot(r2_small_r, levels=cover2, type="class")