pacman::p_load(
rgeoboundaries,
terra, # For working with raster data
exactextractr, # To extract raster data grouping it by sf polygons
sf,
maptiles, # To get tile maps
tmap,
tidyverse,
tidyterra) # To map raster files as they are without turning
# them into dataframes
terraOptions(progress=0) # Per a no veure la barra de progrés#=======================================
# EXTRACCIÓ DE LES DADES
#=======================================
#===========================
# 1. MUNICIPIS SF
#===========================
mun_cat <- st_read(
"MapaMun/divisions-administratives-v2r1-municipis-1000000-20240118.shp")## Reading layer `divisions-administratives-v2r1-municipis-1000000-20240118' from data source `C:\DataPortatil\Laboratori_R\Mapes\deforestation\MapaMun\divisions-administratives-v2r1-municipis-1000000-20240118.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 947 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 260160.2 ymin: 4488767 xmax: 527397.6 ymax: 4747976
## Projected CRS: ETRS89 / UTM zone 31Nmun_cat <- st_transform(mun_cat, crs = 4326) # Projecció del raster
mun_emp <- mun_cat %>%
filter(NOMCOMAR== "Alt Empordà")
#plot(sf::st_geometry(mun_emp))
comar_cat <- st_read( # Mapes comarques de Catalunya
"MapaComarc/shapefiles_catalunya_comarcas.shp")## Reading layer `shapefiles_catalunya_comarcas' from data source
## `C:\DataPortatil\Laboratori_R\Mapes\deforestation\MapaComarc\shapefiles_catalunya_comarcas.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 41 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 0.1594132 ymin: 40.52303 xmax: 3.332542 ymax: 42.8615
## Geodetic CRS: WGS 84comar_cat <- sf::st_transform(comar_cat, crs = 'EPSG:25831')
#===============================
# 2. FOREST DATA
#===============================
# ----------------
# From https://glad.umd.edu/
# https://glad.earthengine.app/view/global-forest-change#bl=off;old=off;dl=1;lon=-50.0283874235115;lat=27.905882522562944;zoom=4;
#----------------
# COBERTURA ARBRES L'ANY 2000
# https://storage.googleapis.com/earthenginepartners-hansen/GFC-2023-v1.11/Hansen_GFC-2023-v1.11_treecover2000_50N_000E.tif
# PERDUES DE COBERTURA DES DE 2000
# https://storage.googleapis.com/earthenginepartners-hansen/GFC-2023-v1.11/Hansen_GFC-2023-v1.11_lossyear_50N_000E.tif
var <- c(
"treecover2000", "lossyear"
)
urls <- paste0(
"https://storage.googleapis.com/earthenginepartners-hansen/GFC-2023-v1.11/Hansen_GFC-2023-v1.11_",
var,
"_50N_000E.tif"
)
forest_data <- lapply(
urls, terra::rast # It loads raster files into R
)
#==============================
# 3. CROP DATA: ENCAIXAR LES DADES DEL RASTER EN LES COORDENADES DE MAPA
# VECTORIAL DE CATALUNYA
#==============================
# Crop the forest data only for Catalonia
# terra::crop takes the raster files and crop them for the extent we choose
cat_data <- lapply(
forest_data,
function(x) {
terra::crop(
x,
terra::vect(mun_cat)
)
}
)El Global Land Analysis and Discovery (GLAD)
El laboratori Global Land Analysis and Discovery (GLAD) del Departament de Ciències Geogràfiques de la Universitat de Maryland investiga mètodes, causes i impactes del canvi global de la superfície terrestre. Les imatges d’observació de la Terra són la font principal de dades i la cobertura terrestre els seus canvies en són el tema d’interès principal.
Global Forest Change 2000-2023
El laboratori Global Land Analysis and Discovery (GLAD) de la Universitat de Maryland, en col·laboració amb Global Forest Watch (GFW), proporciona dades de pèrdua de boscos a escala global actualitzades anualment, derivades mitjançant imatges de sèries temporals de Landsat. Aquestes dades, són un indicador relatiu de les tendències espaciotemporals de la dinàmica de pèrdua de boscos a nivell mundial.
Els responsebles de les dades alerten d’algunes inconsistències que poden contenir, degudes a canvis tecnològics que van haver-hi durant el període que en tot cas mantenen la seva validesa en la identificació de tendències duran el període.
Els arbres es defineixen com a vegetació de més de 5 m d’alçada i s’expressen com a percentatge per cel·la de la quadrícula de sortida com a “percentuale de cobertura d’arbres de l’any 2000”.
La “pèrdua de coberta forestal” es defineix com una pertorbació de substitució, o un canvi d’un estat forestal a un estat no forestal, durant el període 2001-2023.
Les dades de Global Forest Change per a Catalunya
Com ja hem fet al capítol anterior, assaigem per tant aquestes dades per al territori català generant diferents mapes temàtics.
Superfície forestal a Catalunya l’any 2000
cat_data[[1]][cat_data[[1]] < 1 ] = NA
# aggregate() to decrease size of the raster in this case one half ('fact=2' by default)
bosc_2000 <- terra::aggregate(cat_data[[1]], fun="mean")
tmap_mode("view")
tm_shape(bosc_2000)+
tm_raster(alpha = .7,
n= 5,
breaks = c(1,20,40,60,80,100),
palette = c("grey90","darkgreen"),
legend.format = list(text.separator= '-'),
title = "Densitat forestal l'any 2000 (%)") +
tm_shape(comar_cat) +
tm_polygons(id= "nom_comar",
alpha = 0,
popup.vars=c("Superfície comarcal (km2): "= "sup_comar"),
popup.format = list(sup_comar= list(big.mark = ".",
decimal.mark= ",",
digits= 1)),
border.col = "grey35",
lwd= 2.2) +
tmap_options(basemaps = "CartoDB",
alpha= 1)+
tmap_options(max.raster = 1e8)La primera dada que ens facilita aquesta font és la cobertura forestal l’any 2000. El primer mapa ubica la cobertura forestal/arbrada cuantificant-ne la densitat (en %), l’any 2000 que és el punt 0 per a determinar els canvis posteriors.
La segona dada és la ubicació de la pèrdua i l’any que es produeix. Així en el segon mapa, ubiquem la cobertura forestal de 2000, sobreposant-hi les zones on la cobertura s’ha perdut durant el període 2000-2023.
Hi hauria una tercera dada que se’ns facilita que seria les recuperacions de cobertura que s’esdevenen, però només les d’entre 2001 i 2012. No l’utilitzem per l’evident desequilibti que hi hauria amb les dades de pèrdues que cobreixen tot el període.
Visualitarem bosc inicial i pèrdues utilitzant el mode interactiu que permet apropar la imatge per identificar millor les ubicacions i valorar millor l’abast de la perdua. Si juguem amb les capes (a dalt a l’esquerra) veurem com les cel·les de les pèrdues se sobreposen a ubicacions de bosc.
Superfície forestal 2000 i pèrdues en el període 2001-2023
# 4. DEFORESTATION YEARS IN ONE GROUP (1-23)
#---------------------------------------
cat_data[[2]][cat_data[[2]] < 1 ] = NA
# aggregate() to decrease size of the raster in this case one half (fact=2)
bosc_perdut <- terra::aggregate(cat_data[[2]], fun="mean")
tmap_mode("view")
tm_shape(bosc_2000) +
tm_raster(
alpha = 0.6,
n= 2,
breaks = c(1,100),
palette = c("#287233"),
labels = "Any 2000",
legend.format = list(text.separator= '-'),
title = "Presència de bosc"
)+
tm_shape(bosc_perdut)+
tm_raster(
alpha = 0.75,
n= 1,
breaks = c(1,23),
palette = c("#fe0000"),
labels = c("2001-2023"),
legend.format = list(text.separator= '-'),
title = "Pèrdues forestals"
) +
tm_shape(mun_cat)+
tm_polygons(id= "NOMMUNI",
alpha = 0,
popup.vars=c("Comarca: "= "NOMCOMAR",
"Superfície municipi (km2): "= "AREAM5000"),
popup.format = list(AREAM5000= list(big.mark = ".",
decimal.mark= ",",
digits= 1)),
border.col = "grey70",
lwd= .8)+
tm_shape(comar_cat)+
tm_borders(col = "steelblue",
lwd=1.6)