¿Cómo analizar coberturas y realizar planes de conservación?

Se trata de hacer un análisis sobre las tendencias históricas de pérdida y ganancia de la cobertura forestal desde 1999 hasta 2019. El modelo que vamos a crear está hecho por criterios globales, lo que significa que algunos tipos de bosques han sido incluidos en él, pero no significa que haya bosques primarios en absoluto.Por tanto algunas cosas quedaran dentro del criterio de uds los expertos.

Atención: ¿Qué tal crear una Área de Conservación para tu Polígono? ¿Lo has hecho ya? Por favor, ten cuidado, no podrás avanzar sin tenerlo.

Antes de empezar: creación del espacio de trabajo

Antes de empezar, es necesario instalar los paquetes de R requeridos. Por favor, realiza una sola ejecución de instalación por computadora. Recuerda que la instalación se realiza con la función install.packages() y el nombre del paquete entre comillas. Por ejemplo, install.packages(“maptools”). Luego, puedes cargar los paquetes en tu sesión de R utilizando la función library(), como se muestra en el siguiente ejemplo: library(maptools). A continuación los paquetes.

library(rgdal)
library(sp)
library(raster)
library(gdalUtilities)
library(parallel)
library(sf)
library(mapview)
library(gfcanalysis)

Dado nuestros conocimientos previos, ahora podremos avanzar rápidamente en los siguientes pasos sin mayor explicación.

nota: Debes tener una carpeta llamada Reserva dentro de “clase_conservacion”.

Seleccionando nuestro polígono en el espacio de trabajo

Recuerda que para esta altura ya debes de haber creado tu polígono en QGIS. Entonces, vamos a importarlo. Debes ir punto a punto. Lastimosamente el siguiente codigo es muy delicado y como resultado no podemos dividirlo o explicarlo del todo.

setwd("D:/biato/wikiaves/clase_conservacion/Reserva")
getwd()
## [1] "D:/biato/wikiaves/clase_conservacion/Reserva"
park <- st_read("D:/biato/wikiaves/clase_conservacion/Reserva/Corredor_gigante.shp")
## Reading layer `Corredor_gigante' from data source 
##   `D:\biato\wikiaves\clase_conservacion\Reserva\Corredor_gigante.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.04547 ymin: 3.621559 xmax: -76.01371 ymax: 3.653911
## Geodetic CRS:  WGS 84
st <- "+proj=longlat +datum=WGS84 +no_defs"
park <- st_transform(park, crs = st)


mapview( park, layer = "Admin. level 0", color = "000000")
output_folder <- setwd("D:/biato/wikiaves/clase_conservacion/Reserva")

Establecimos algunos parámetros de grabado y ahora realizaremos la descarga de la información de cobertura correspondiente a nuestra área.

tiles <- calc_gfc_tiles(park) 

tiles
## class       : SpatialPolygons 
## features    : 1 
## extent      : -80, -70, 0, 10  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
plot(tiles, col = "#EF2929", lwd = 1, main = "Plot Area + Tiles", xlab = "Long", ylab = "Lat") ## please input 1 or 2 depends on the features result - it is the number of tiles needed
plot(park,  col = "#555753", add = TRUE)

box()

####### Descarga de rastres##### Solo se descarga una vez
download_tiles(
  tiles,
  output_folder,
  images = c("treecover2000", "lossyear", "gain", "datamask"),
  dataset = "GFC-2019-v1.7"
)

#######Corre haste arriba #####
gfc_extract <- extract_gfc(park, output_folder, dataset = "GFC-2019-v1.7", overwrite=TRUE, to_UTM= F, filename="Park_extract.tif")
#Going on

A continuación tendremos los valores del cambio de vegetación en el área dependiendo nuestro criterio de porcentaje de cobertura. Vamos a exportarlos de la siguiente manera. Ten en cuenta que se graban en .csv.

forest_threshold <- 75 

gfc75 <- threshold_gfc(gfc_extract, forest_threshold=forest_threshold, filename="Park_ex_extract_2.tif", overwrite=TRUE)

gfc75_stats <- gfc_stats(park, gfc75, dataset = "GFC-2019-v1.7")

write.csv(gfc75_stats$loss_table, file="park_GFC75_losstable.csv", row.names=FALSE)
write.csv(gfc75_stats$gain_table, file="park_GFC75_gaintable.csv", row.names=FALSE)

#Video showing the "devastation"

gfc75_annual_stack <- annual_stack(gfc75, dataset = "GFC-2019-v1.7")
writeRaster(gfc75_annual_stack, filename="park_GFC75_annual.tif", overwrite=TRUE, bylayer=TRUE)

park$label <- "ZOI"

animate_annual(park, gfc75_annual_stack, out_dir='output_folder', site_name='MNA_GFC75',type = "html")

nota: en tu carpeta encontraras las carperas: park_GFC75_losstable.csv = cantidad en hectáreas de bosque pedido y park_GFC75_gaintable.csv = cantidad de bosque ganado en el área de estudio.

nota 2: aparecerá en la carpeta un archivo llamado “gfc_animation.html” le das doble click y listo. También, aparecerá una carpeta llamada “gfc_animation_imgs” donde podrás encontrar todas las imágenes o mapas que crearon el video.

Exportar un output en forma raster para QGIS

Ahora deseamos exportar la informacion de tal forma que sea util para otras plataformas de SIG, en nuestro caso seran exportadas en formato raster.

Forest<- raster::raster("Park_ex_extract_2.tif", band = 1)
Lost<- raster::raster("Park_ex_extract_2.tif", band = 2)
gained<- raster::raster("Park_ex_extract_2.tif", band = 3)
fondo_1<- raster::raster("Park_ex_extract_2.tif", band = 4)
fondo_2<- raster::raster("Park_ex_extract_2.tif", band = 5)

writeRaster(Forest,
            "D:/biato/wikiaves/clase_conservacion/Reserva/bosque1.tiff",
            format="GTiff", overwrite=T)

writeRaster(Lost,
            "D:/biato/wikiaves/clase_conservacion/Reserva/perdida.tiff",
            format="GTiff", overwrite=T)

writeRaster(gained,
            "D:/biato/wikiaves/clase_conservacion/Reserva/ganado.tiff",
            format="GTiff", overwrite=T)


mapview(Forest, layer = "Admin. level 0")  ###please I encourage your to type every raster made before, making an idea of them visualization.

nota: donde, “bosque1.tiff” = Cobertura boscosa total, “perdida.tiff” = total de bosque perdido hasta el 2020 y “ganado.tiff” = total de bosque ganado en el area.

Análisis adicionales de vegetación (Toma de decisiones finales)

Ahora analizaremos algunos otros aspectos de la cobertura boscosa, cosa que nos brindara información para toma de decisiones más asertiva.

Primero instalaremos algunos nuevos paquetes. Para esta altura ud debe contar con el conocimiento para entender que el siguiente código es solo ejecución del paquete y debe realizar una única instalación previa.

library(landscapemetrics)
library(landscapetools)
library(ggplot2)

Ahora importaremos de nuevo el raster de bosque creado de nuestra area

primario <- raster("D:/biato/wikiaves/clase_conservacion/Reserva/bosque1.tiff")

### 
primario_Park <- crop(x = primario, y = park) ## it is to cut it off
plot(primario_Park)

A continuación realizaremos algunos análisis de la cobertura,

sr <- "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs" 
primario_Park2 <- projectRaster(primario_Park, crs = sr, method ="ngb")

#Pulling your weight forward

check_landscape(primario_Park)
##   layer        crs   units   class n_classes OK
## 1     1 geographic degrees integer         2  ✖
sr <- "+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs" 
primario_Park2 <- projectRaster(primario_Park, crs = sr, method ="ngb")

check_landscape(primario_Park2)
##   layer       crs units   class n_classes OK
## 1     1 projected     m integer         2  ✔
#extract some metrics

metrics_ <-calculate_lsm(primario_Park2, #a raster object
                         what=c("lsm_c_ai", # metric name
                                "lsm_c_frac_mn", #  metric name
                                "lsm_l_division",# metric name
                                "lsm_l_pd")) # metric name

head(metrics_)
## # A tibble: 6 × 6
##   layer level     class    id metric    value
##   <int> <chr>     <int> <int> <chr>     <dbl>
## 1     1 class         0    NA ai       85.9  
## 2     1 class         1    NA ai       88.1  
## 3     1 class         0    NA frac_mn   1.04 
## 4     1 class         1    NA frac_mn   1.03 
## 5     1 landscape    NA    NA division  0.639
## 6     1 landscape    NA    NA pd       19.3
metrics_names<-lsm_abbreviations_names #store the names
tail(metrics_names) #show the last 6 observations
## # A tibble: 6 × 5
##   metric name                      type                 level     function_name
##   <chr>  <chr>                     <chr>                <chr>     <chr>        
## 1 sidi   simpson's diversity index diversity metric     landscape lsm_l_sidi   
## 2 siei   simspon's evenness index  diversity metric     landscape lsm_l_siei   
## 3 split  splitting index           aggregation metric   landscape lsm_l_split  
## 4 ta     total area                area and edge metric landscape lsm_l_ta     
## 5 tca    total core area           core area metric     landscape lsm_l_tca    
## 6 te     total edge                area and edge metric landscape lsm_l_te
class_metrics<-calculate_lsm(list(primario_Park2), # add multiple maps 
                             level = c("class"), #patch- landscape other options
                             type="aggregation metric") #get all the metrics for the given level

Estos análisis nos otorgan una variedad de resultados, pero para nuestro ejercicio solo vamos a necesitar los siguientes:

Resultado 1 : número de parches de bosque que están en el área. Graficando el paisaje con los parches etiquetados con el id de parche correspondiente, esto quiere decir que cuando un área se encuentra del mismo color se consideran unidas y pertenecientes al mismo parche de bosque.

patches <- show_patches(primario_Park2,labels = FALSE,class = 1) #shows the patches of class 3
patches
## $layer_1

Resultado 2 : Tamaño de los parches de bosque, entre más amarillo, mayor si área.

show_lsm(primario_Park2,labels = F,class = 1,what = "lsm_p_area") #show the patches in colours based on the area, good to observe where that largest patch is located  
## $layer_1

Resultado 3 : Borde de bosque efectivo a lo largo de toda la cobertura.

show_cores(primario_Park2,labels = F,class = 1)
## $layer_1