Introducción

¡Hola! Este sencillo script ayudará a echar un vistazo al análisis espacial de cuencas con R.

### Ejercicio 1: Explorando las Categorías de Uso de Suelo
# 1. **Objetivo:** Identificar y explorar las diferentes categorías de uso de suelo en tu ráster.
# 2. **Pasos:**
#   - Carga tu archivo ráster de uso de suelo.
# - Utiliza la función `freq()` para identificar las categorías de uso de suelo únicas y sus frecuencias.
# - Muestra los resultados.

## Carga la librería terra
# install.packages("terra") #si es necesario
library(terra)
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.7.78
# Carga tu archivo ráster de uso de suelo, para Biobio
lc <- rast("reproj_lc_biobio.tif")

# Plotear el raster
plot(lc)

# 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)

#Plot
plot(lc_biobio_mask)

# Obtén la frecuencia de cada categoría de uso de suelo
uso_suelo_freq <- freq(lc_biobio_mask)

# Muestra las categorías de uso de suelo y sus frecuencias
print(uso_suelo_freq)
##    layer value   count
## 1      1     0   77910
## 2      1   110     179
## 3      1   120     804
## 4      1   130  872590
## 5      1   140    7221
## 6      1   150  139520
## 7      1   211  928767
## 8      1   212 8220662
## 9      1   222      42
## 10     1   241 3955419
## 11     1   251 5518601
## 12     1   311 2387526
## 13     1   312       1
## 14     1   320 3452888
## 15     1   330   45503
## 16     1   410 5476576
## 17     1   420  170797
## 18     1   450     524
## 19     1   510    9730
## 20     1   610  217088
## 21     1   620     319
## 22     1   630  230753
## 23     1   640   28430
## 24     1   800  163327
## 25     1   910      23
## 26     1   920   15839
## 27     1   931  224394
## 28     1   932  758177
## 29     1  1010   23205
## 30     1  1020       1
## 31     1  1210     172
### Ejercicio 2: Reclasificando Categorías de Uso de Suelo
# 1. **Objetivo:** Reclasificar el ráster de uso de suelo en categorías más amplias.
# 2. **Pasos:**
# - Define una matriz de reclasificación (por ejemplo, combinando diferentes usos de suelo en categorías más amplias como "Urbano," "Agrícola," "Bosque").
# - Usa la función `classify()` para reclasificar el ráster.
# - Grafica el ráster original y el reclasificado.

# Define la matriz de reclasificación (ejemplo: códigos de uso de suelo a categorías más amplias)
# 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 
matriz_reclas <- matrix(c(0, 99, 0, #los valores entre 0-90 --> 0
                           100, 499, 100,
                           500, 999, 500,
                           1000, 1499, 1000),
                         ncol=3, byrow=TRUE)

# Reclasifica el ráster
lc_reclas <- classify(lc_biobio_mask, rcl=matriz_reclas)

# Grafica los rásters original y reclasificado
par(mfrow=c(1,2))
plot(lc_biobio_mask, main="Uso de Suelo Original")
plot(lc_reclas, main="Uso de Suelo Reclasificado")

# 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]]"
# reproyección
lc_biobio_mask<-project(lc_biobio_mask, "EPSG:32719") #UTM 19S
## |---------|---------|---------|---------|=========================================                                          
# Crop + Mask otra vez
r2<- crop(lc_biobio_mask, cuenca_biobio) #crop
r2<- mask(r2, cuenca_biobio) #mask

# Plots
par(mfrow=c(1,1))

plot(r2)
plot(lc_biobio_mask, add=T) #intentar cambiar titulo, etc.