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) # To hide progress bar

Introducció

El present treball s’emmarca en la labor continuada d’identificació i avaluació de fonts de dades que poden ser apropiades per a construir indicadors temàtics que mesurin l’acompliment dels Objectius de Desenvolupament Sostenible a l’àmbit local (municipal i comarcal), que duem a terme des de Gaia Consultoria i Dialòguia.

El fruit d’aquest treball es va plasmant en el visor multipanell Els mapes dels ODS als municipis de Catalunya i en altres instruments que hem creat per a les institucions per a les quals treballem.

En aquest marc doncs, hem estat utilitzant dades de tipus administratius (INE, IDESCAT, altres instituts d’estadística), dades web (com ara el visors de ECOEMBES), dades de d’organismes publics (Agència Catalana de l’Aigua, Institut Català d’Energia, el Ministerio de Interior etc.), i institucions privades (com ara l’Observatori Català de la Mobilitat, o el Fons Català per al Desenvolupament).

La necessitat d’obtenir dades amb granularitat local que puguin donar resposta a les institucions del territori d’àmbit municipal i comarcal ens exigeixen sondejar fonts que les poden aportar. En aquest sentit les dades generades pels satèl·lits poden ser útils per donar cos a més d’un indicador.

Pel que fa a aspectes de programació, a dalt facilitem l’enllaç per poder visualitzar el codi base. Respecte a això, devem molt a Milos Popovic, per la quantitat de saviesa que proporciona amb els seus articles i tutorials dedicats a la programació en R i, específicament, a la creació de cartografies a partir de dades obertes.

Mapejar dades de cobertura del territori amb Esri Sentinel-2

Sobre les dades

Els mapes de cobertura del sòl (LULC) són una eina cada cop més important per als qui prenen decisions en molts sectors industrials i nacions en desenvolupament d’arreu del món. La informació proporcionada per aquests mapes ajuda a informar les decisions polítiques i de gestió de la terra mitjançant una millor comprensió i quantificació dels impactes dels processos terrestres i de l’activitat humana.

En el nostre cas, hem accedit a les dades de Sentinel-2 Land Cover Explorer de ESRI.

Mapa de la cobertura territorial de 2022

Les dades de Sentinel-2 L2A tenen una enorme granularitat, mapejant cel·les de 10 m. Encara se’n està avaluant la precisió, tanmateix el potencial és molt gran pel detall que pot arribar a donar respecte a les dades administratives. Les dades del catastre poden contenir imprecisions per la poca freqüència amb què s’actualitzen; terres classificades com “bosc” sovint ja no ho són (incendis, urbanitzacions, canvis d’us, etc.). Les dades d’aquest origen ofereixen doncs, a més del nivell de detall, la màxima actualització.

D’altra banda s’identifiquen 9 diferents categories d’observació del satèl·lit: Water, Trees, Flooded vegetation, Crops, Built Area, Bare ground,Snow/Ice, Clouds, Rangeland. A nivel de Catalunya se’n identifiquen 7: “Aigua”, “Arbres”, “Aiguamolls”, “Cultius”, “Cobertura artificial”, “Terra nua”, “Pastures”.

Mapes directes de les dades

 # 2. DOWNLOAD RASTER FILES
#-------------------------
#==================================
# ADREÇA WEB DE SENTINEL 2
### https://livingatlas.arcgis.com/landcoverexplorer/#mapCenter=9.40600%2C45.75100%2C6&mode=swipe&timeExtent=2018%2C2022&renderingRule=0
#==================================

#url <- 
#  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2023/31T# _20230101-20240101.tif"

#  download.file(
#    url = url,
#    destfile = "31T_20230101-20240101.tif",
#    mode = "wb"
#  )

forest22_data <- terra::rast("31T_20230101-20240101.tif") # It loads raster files into R

# Reference of de mun_cat
crs(forest22_data) <- 'EPSG:25831' 

# CROPPIN RASTER INSIDE mun_cat

for22_region <- terra::crop(forest22_data,
            terra::vect(mun_cat))

# FIRST VISUALIZATION

terra::plot(
  for22_region,
)
plot(comarques_df,
  lwd=1.5, border= "white",
  add=TRUE
  )
legend("bottomright",
       inset=.02, title="Tipus de cobertura",
       legend=c("Arbres", "Artificial", "Cultius","Pastura"),
       fill= c("#287233", "#e51d2e","#ff7907","#f8e596"),
       cex=0.8, bg='lightblue', xpd = TRUE
       )
title(main = "Cobertura del terreny a Catalunya (Esri - Sentinel-2)",
      outer = TRUE, line = -1
      )
mtext("Esri, Impact Observatory", side=1, line=1, , adj=0.7, cex=1)

Aquest tipus de dades permet visualtzacions molt detallades, a nivell de municipi o de barri, com en aquest cas, on representem les diferent tipologies de cobertura al territori de l’Alt Empordà.

# MAPA ALT EMPORDÀ
#===========================

aemporda_sf <- mun_cat %>% 
  filter(NOMCOMAR== "Alt Empordà")

aemporda_region <- terra::crop(forest22_data,
            terra::vect(aemporda_sf))

bbox <- sf::st_bbox(
  aemporda_sf
)

bg_map <- maptiles::get_tiles(
  x = bbox,
  provider = "Esri.NatGeoWorldMap",
  zoom = 10,
  crop = TRUE,
  project = FALSE
)


p1 <- ggplot() +
  tidyterra::geom_spatraster_rgb( # tidyterra passes raster files to images
    data = bg_map,
    maxcell = terra::ncell( # set the cell number.
      bg_map
    ),
    alpha = 1
  )+
  tidyterra::geom_spatraster(
    data = as.factor(
      aemporda_region
    ),
    maxcell = terra::ncell(
      bg_map
    ),
    alpha = .6
  ) +
  scale_fill_manual(name= "",
               labels= c("Aigua",
                         "Arbres",
                         "Aiguamolls", 
                         "Cultius",
                         "Artificial",
                         "Terra nua",
                         "Pastura"),
               values= c(
                '1'="lightblue",
                '2'="#287233",
                '4'="lightgreen",
                '5'="#ff7907",
                '7'="#e51d2e",
                '8'="grey",
                '11'="#773813")) +
    geom_sf(
    data = aemporda_sf,
    color = "white",
    fill = "transparent",
    lwd= 0.75
    ) +
  ggtitle("Alt Empordà") +
  labs(caption = "Esri, Impact Observatory")+
  coord_sf(expand = FALSE)+
  theme(
    plot.title = element_text(size = 16),
    plot.caption = element_text(size = 12, vjust = 45, colour= "black"),
    legend.position = "bottom"
  )+
  guides(fill=guide_legend(title = NULL,nrow=1))

p1

Òbviament, si volem visualitzar les dades sense agrupar, potser l’opció millor és fer-ho mtjançant el visor Sentinel-2 Land Cover Explorer, excepte que necessitem visualitzar zones allunyades entre elles i les volem veure una al costat de l’altre com en el cas següent:

# MAPA BAIX LLOBREGAT
#=================================

bllobregat_sf <- mun_cat %>% 
  filter(NOMCOMAR== "Baix Llobregat")

bllobregat_region <- terra::crop(forest22_data,
            terra::vect(bllobregat_sf))

bbox <- sf::st_bbox(
  bllobregat_sf
)

bg_map <- maptiles::get_tiles(
  x = bbox,
  provider = "Esri.NatGeoWorldMap",
  zoom = 10,
  crop = TRUE,
  project = FALSE
)


p2 <- ggplot() +
  tidyterra::geom_spatraster_rgb( 
    data = bg_map,
    maxcell = terra::ncell(
      bg_map
    ),
    alpha = 1
  )+
  tidyterra::geom_spatraster(
    data = as.factor(
      bllobregat_region
    ),
    maxcell = terra::ncell(
      bg_map
    ),
    alpha = .6
  ) +
    scale_fill_manual(name= "",
               labels= c("Aigua",
                         "Arbres",
                         "Aiguamolls", 
                         "Cultius",
                         "Artificial",
                         "Terra nua",
                         "Pastura"),
               values= c(
                '1'="lightblue",
                '2'="#287233",
                '4'="lightgreen",
                '5'="#ff7907",
                '7'="#e51d2e",
                '8'="grey",
                '11'="#773813")) +

    geom_sf(
    data = bllobregat_sf,
    color = "white",
    fill = "transparent",
    lwd= 0.75
    ) +
  ggtitle("Baix Llobregat") +
  labs(caption = "Esri, Impact Observatory")+
  coord_sf(expand = FALSE)+
  theme(
    plot.title = element_text(size = 16),
    plot.caption = element_text(size = 12, vjust = 45, colour= "black"),
    legend.position = "bottom"
  )

library(patchwork)
p1 + p2 +

     plot_annotation(title = "Diferents cobertures del sòl a dues comarques de Catalunya",theme = theme(plot.title = element_text(size = 18, 
                                                                                                                     face = 'bold',
                                                                                                                     colour = 'dodgerblue3',
                                             hjust = 0.5)))

Mapes agregats per municipi de cobertura del sòl

En canvi, podem necessitar aquestes dades agrupades per poder-les visualitar i analitzar per divisions administratives, com ara els municipis, per tal de poder-ne obtenir indicadors útils a ser incorporats al nostre sistema d’indicadors sobre ODS.

Això es pot aconseguir amb unes poques manipulacions, generant mapes i altres tipus de gràfiques, com es pot veure en el mapa següent on visualitzem la cobertura arbrada/forestal a Catalunya l’any 2022 per municipi:

Mapa interactiu del % de cobertura forestal sobre la sup. del municipi 

# 1. OPCIÓ 1: GENERANT NOVES DADES (ÉS MOLT LENT)
#================================================

#=====================================
# EXTRACCIÓ DE LES DADES DEL RASTER
#=====================================

#for22 <- terra::crop(forest22_data,
#            terra::vect(mun_cat))

#lc <- exactextractr::exact_extract(
#      for22,
#      mun_cat,
#      function(df) {
#        df |>
#          dplyr::group_by(
#            CODIMUNI      
#          ) |>
#          dplyr::summarize(
#            area_km2 = sum(
#              coverage_area[value== '2'] / 1e6 # https://www.arcgis.com/home/item.html?id=cfcb7609de5f478eb7666240902d4d3d
#            )
#          )
#      },
#      summarize_df = T,      # It creates a dataframe with the extracted values
#      coverage_area = T,     # The altenative is to calculate coverage_fraction
#      include_cols = "CODIMUNI"
#    )

#===================================================
# 1. OPCIÓ 2: RECUPERANT DADES JA GENERADES (MÉS RÀPID)
#===================================================

#writexl::write_xlsx(lc, "lc.xlsx")

lc <- readxl::read_xlsx("lc.xlsx")

#===================================
# 2. DATAFRAME DE COBERTURA FORESTAL
#===================================

cat_forest_df <- mun_cat %>% 
  inner_join(lc, by= "CODIMUNI") %>% 
  mutate(frac_cover= round(area_km2 / AREAM5000 * 100,1))

#====================================
# 3. MAPA DE COBERTURA PER MUNICIPIS
#====================================
tmap_mode("view")

tm_shape(cat_forest_df)+
  tm_fill(col="frac_cover",
          id="NOMMUNI",
          #style = "fixed",
          #n=5,
          #breaks = c(0,50,100,1000,2000,5000),
          palette = c("papayawhip","darkgreen"),
          title = "Cobertura forestal <br>any 2022 (%)",
          legend.format = list(text.separator= '-'),
          textNA = 'Sense dades',
          popup.vars=c("Cobertura forestal (%): "= "frac_cover"),
          popup.format = list(frac_cover= list(digits = 1))
  )+
  tm_borders(col="grey40",lwd=0.5) +
  tmap_options(basemaps = "Esri.WorldImagery",
               basemaps.alpha = 1)

O en aquest altre mapa on representem la cobertura % de cobertura artificial del sòl als municipis de Catalunya:

Mapa interactiu del % de superfície amb cobertura artificial sobre la sup. total del municipi 

#=====================================================================
# 1. OPCIÓ 1: EXTRACCIÓ INICIAL (ÉS LLARGA: ES POT OBVIAR UN COP FETA)
#=====================================================================


#for22 <- terra::crop(forest22_data,
#            terra::vect(mun_cat))

#lc2 <- exactextractr::exact_extract(
#      for22,
#      mun_cat,
#      function(df) {
#        df |>
#          dplyr::group_by(
#            CODIMUNI
#          ) |>
#          dplyr::summarize(
#            area_km2 = sum(
#              coverage_area[value== '7'] / 1e6 # https://www.arcgis.com/home/item.html?id=cfcb7609de5f478eb7666240902d4d3d
#            )
#          )
#      },
#      summarize_df = T,      
#      coverage_area = T,
#      include_cols = "CODIMUNI"
#    )

#======================================================
# 1. OPCIÓ 2: RECUPERANT DADES JA GENERADES (MÉS RÀPID)
#======================================================

# writexl::write_xlsx(lc2, "lc2.xlsx")

lc2 <- readxl::read_xlsx("lc2.xlsx")

#===================================
# 2. DATAFRAME DE COBERTURA FORESTAL
#===================================

cat_constr_df <- mun_cat %>% 
  inner_join(lc2, by= "CODIMUNI") %>% 
  mutate(frac_cover= ifelse(round(area_km2 / AREAM5000 * 100,1)>100,100,
                            round(area_km2 / AREAM5000 * 100,1))
          )

#====================================
# 3. MAPA DE COBERTURA PER MUNICIPIS
#====================================
tmap_mode("view")

tm_shape(cat_constr_df)+
  tm_fill(col="frac_cover",
          id="NOMMUNI",
          #style = "fixed",
          #n=5,
          #breaks = c(0,50,100,1000,2000,5000),
          palette = c("papayawhip","darkred"),
          title = "% Superfície artificial <br>any 2022",
          alpha= 0.5,
          legend.format = list(text.separator= '-'),
          textNA = 'Sense dades',
          popup.vars=c("Superfície edificada (%): "= "frac_cover"),
          popup.format = list(frac_cover= list(digits = 1))
  )+
  tm_borders(col="grey40",lwd=0.5) +
  tmap_options(basemaps = "Esri.WorldImagery",
               basemaps.alpha = 1)

De l’observació de les visualitzacions, es desprèn que poden existir diferències respecte als registres administratius, sobre tot pel que fa a cobertura artificial. Primer, perquè es tracta d’aproximacions que realitza un algoritme sobre la base de les imatges rebudes pel satèl·lit. A més a més, cal tenir present que l’algoritme que classifica pot tenir un criteri diferent del que regeix un catastre: s’identifica com a cobertura artificial tot tipus de cobertura, fins i tot uns hibernacles de plàstic (això queda evident observant Vilassar de Mar, per exemple), uns terrenys d’urbanitzacions amb jardins també (fixeu-vos en el 40% de cobrtura artifical de Begur…), o terrenys dins un centre urbà que no tinguin prou dimensions (un petit parc…). I així altres aspectes. Per aquesta raó hem deixat un cert nivell de transparència en aquest darrer mapa per poder comprovar de molt a prop aquest indicador.

En qualsevol cas, es tracta de petites diferències que no creiem que incideixin de manera significativa en els valors generats. Al contrari, n’enriqueixen el significat i hi introdueixen una perspectiva diferent.

Conclusió

Aquesta font de dades ens permet tenir una alternativa als arxius administratius, quan aquests no ens acaben de satisfer o no ens aporten la informació que necessitem. Molt interessant per determinar la cobertura forestal real i actualitzada, així com la cobertura artificial dels municipis i comarques catalans.