---
title: "Quantificació de l'área residencial inundable amb PR=500 (ODS 11), als municipis de Catalunya amb dades del <i>Servicio de Zonas Inundables Asociadas a Periodos de Retorno</i> del MITECO i dades del SIOSE"
subtitle: "Notes tècniques 9"
author: '[Stefano Sanfilippo](https://linkedin.com/in/stefano-sanfilippo-2311881a2){target="_blank"}'
date: "`r Sys.Date()`"
output: 
  prettydoc::html_pretty:
    theme: cayman
    self_contained: no
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(warning = FALSE)
knitr::opts_chunk$set(cache = FALSE)


```

```{r}

pacman::p_load(
  raster,        # raster data
  terra,         # raster data
  exactextractr, # Extracció de raster data agregant-les per pol. sf
  sf,
  tmap,
  mapview,
  leafsync,
  tidyverse,
  paletteer,
  ggthemes,
  classInt)

terraOptions(progress=0) # Per a no veure la barra de progrés
```



```{r include=FALSE}
#=======================================
# EXTRACCIÓN DE DATOS
#=======================================

#===========================
# 1. MUNICIPIOS SF
#===========================

mun_cat <- st_read(
"../repo/divisions-administratives/divisions-administratives-v2r1-municipis-1000000-20240118.shp")
mun_cat <- sf::st_as_sf(mun_cat)
mun_cat <- st_transform(mun_cat, crs = 4326)



```


```{r include=FALSE}

#=================================================================
# EXTRACCIÓ DE LA DADES D'INUNDABILITAT PER PROVÍNCIA DE CATALUNYA
#=================================================================

#https://www.miteco.gob.es/es/cartografia-y-sig/ide/descargas/agua/zi-lamina.html

inunEsp <- st_read("../repo/laminaspb-q500/Q500_2Ciclo_PB_20240912_ETRS89.shp")
inunEsp <- sf::st_as_sf(inunEsp)
inunEsp <- st_transform(inunEsp, crs = 4326)

# Filtrado de los polígonos que coinciden con los de la prov. de Barcelona
cod.prov <- c('08','17','25','43')

sf_use_s2(FALSE)
list_inun <-lapply(cod.prov, function(x) {
  
  a <- st_intersects(inunEsp, mun_cat[mun_cat$CODIPROV== x, ])
  b <- inunEsp[lengths(a) > 0,]
})
sf_use_s2(TRUE)

list_inun_simpl <- lapply(list_inun, function(x) {
  rmapshaper::ms_simplify(x,keep = 0.1,keep_shapes = TRUE)
                          })



# lapply(list_inun_simpl, function(x) plot(st_geometry(x)))

```



```{r include=FALSE}

# LECTURA DADES DESCARREGADES DEL SIOSE 
#----------------------------------------

cob_cat <- st_read(
"../repo/SIOSE_Catalunya_2014_H31_GPKG/SIOSE_Catalunya_2014_H31.gpkg"
  )

# FILTRADO DE LAS ÁREAS RESIDENCIALES (VER MANUAL)
zona_residencial <- cob_cat %>% filter(CODIIGE %in% c(111:113))

# AJUSTE PROYECCION
zona_residencial <- st_transform(zona_residencial, crs = 4326)

#-----------------------------------------------
# Subset poligons segons províncies de Catalunya
#-----------------------------------------------

sf_use_s2(FALSE)
list_resid <-lapply(cod.prov, function(x) {
  
  a <- st_intersects(zona_residencial, mun_cat[mun_cat$CODIPROV== x, ])
  b <- zona_residencial[lengths(a) > 0,]
})
sf_use_s2(TRUE)



residencial_list <- lapply(list_resid, function(x) {
  x %>% 
  mutate(Cobertura= case_match(CODIIGE,
                               111~ "Casc urbà",
                               112~ "Eixample",
                               113~ "Discontinu"),
         Cobertura= factor(Cobertura,
                           levels= c(
                             "Casc urbà",
                             "Eixample",
                             "Discontinu")
                           ))
})

# lapply(residencial_list, function(x) plot(st_geometry(x)))

```


## Visualització de les zones inundables a les quatre províncies de Catalunya

```{r}

mapview(residencial_list[[1]], 
        map.types= c("OpenStreetMap","Esri.WorldImagery"),
        zcol= "Cobertura",
        col.regions= c("#7c0000", "#bc634f", "#e3b0a3"),
        popup = FALSE,
        label= FALSE,
        legend= TRUE,
        layer.name= "Cobertura residencial <br>Prov. de Barcelona") +
mapView(list_inun_simpl[[1]], color = "#9accff",
        popup= FALSE,
        label= FALSE,
        legend= TRUE,
        layer.name= "Superficie inundable <br>Prov. de Barcelona")
```

<p style= "margin-bottom:0; font-weight:bold; font-size: 16px; color:dodgerblue;">Visor del risc d'inundació amb màxim PR (500). Província de Girona</p>
```{r}

mapview(residencial_list[[2]], 
        map.types= c("OpenStreetMap","Esri.WorldImagery"),
        zcol= "Cobertura",
        col.regions= c("#7c0000", "#bc634f", "#e3b0a3"),
        popup = FALSE,
        label= FALSE,
        legend= TRUE,
        layer.name= "Cobertura residencial <br>Prov. de Girona") +
mapView(list_inun_simpl[[2]], color = "#9accff",
        popup= FALSE,
        label= FALSE,
        legend= TRUE,
        layer.name= "Superficie inundable <br>Prov. de Girona")

```

\

<p style= "margin-bottom:0; font-weight:bold; font-size: 16px; color:dodgerblue;">Visor del risc d'inundació amb màxim PR (500). Província de Lleida</p>
```{r}
mapview(residencial_list[[3]], 
        map.types= c("OpenStreetMap","Esri.WorldImagery"),
        zcol= "Cobertura",
        col.regions= c("#7c0000", "#bc634f", "#e3b0a3"),
        popup = FALSE,
        label= FALSE,
        legend= TRUE,
        layer.name= "Cobertura residencial <br>Prov. de Lleida") +
mapView(list_inun_simpl[[3]], color = "#9accff",
        popup= FALSE,
        label= FALSE,
        legend= TRUE,
        layer.name= "Superficie inundable <br>Prov. de Lleida")

```

<p style= "margin-bottom:0; font-weight:bold; font-size: 16px; color:dodgerblue;">Visor del risc d'inundació amb màxim PR (500). Província de Tarragona</p>
```{r}

mapview(residencial_list[[4]], 
        map.types= c("OpenStreetMap","Esri.WorldImagery"),
        zcol= "Cobertura",
        col.regions= c("#7c0000", "#bc634f", "#e3b0a3"),
        popup = FALSE,
        label= FALSE,
        legend= TRUE,
        layer.name= "Cobertura residencial <br>Prov. de Tarragona") +
mapView(list_inun_simpl[[4]], color = "#9accff",
        popup= FALSE,
        label= FALSE,
        legend= TRUE,
        layer.name= "Superficie inundable <br>Prov. de Tarragona")

```

## Càlcul i visualització de la Taxa de d'àrea residencial en risc d'inundació amb PR = 500

```{r}
#===================================================================
# CÀLCUL EXTENSIÓ/MUNICIPI DE LES ZONES RESIDENCIALS
#===================================================================
residencial_cat <- bind_rows(residencial_list) # Llista => dataframe

#--------------------------------------------------------------------
# Interseccio del df de les zones residencials amb la el df dels plígons dels municipis

sf_use_s2(FALSE)
resid_mun <- st_intersection(mun_cat,
                             st_geometry(residencial_cat)) %>% 
  st_as_sf() %>% 
  mutate(
    area= units::set_units(st_area(.), "ha")) %>% 
  st_drop_geometry() %>% 
  group_by(CODIMUNI, NOMMUNI) %>% 
  summarise(area_resid = sum(area))
sf_use_s2(TRUE)

#===================================================================
# CÀLCUL EXTENSIÓ/MUNICIPI DE LES ZONES RESIDENCIALS INUNDABLES
#===================================================================
#--------------------------------------------------------------------
# Interseccio del df de les zones residencials amb la llista de zones inundables => zones residencials inundables
sf_use_s2(FALSE)
list_zones_res_inun <- lapply(list_inun_simpl, function(x) {
  st_intersection(st_geometry(x), st_geometry(residencial_cat)) %>% 
  st_as_sf()
})
sf_use_s2(TRUE)

zones_res_inun <- bind_rows(list_zones_res_inun) # Llista => dataframe

#-------------------------------------------------------------
# Interseccio zones residencials inundables amb municipis =>
# => zones residencials inundables per municipi
sf_use_s2(FALSE)
mun_inund <- st_intersection((mun_cat), 
                             st_geometry(zones_res_inun)) %>% 
  st_as_sf() %>% 
  mutate(
    area= units::set_units(st_area(.), "ha")) %>% 
  st_drop_geometry()
sf_use_s2(TRUE)
#-------------------------------------------------------------
# S'agreguen els valors de zones inundables per municipi
area_inund <- mun_inund %>% 
  #mutate(CODIMUNI= str_sub(CODIMUNI,1,5)) %>% 
  group_by(CODIMUNI, NOMMUNI) %>%
  summarise(area_inund= sum(area))
  
mun_cat_inund <- mun_cat %>% 
  left_join(resid_mun[,c(1,3)], by= "CODIMUNI") %>%
  left_join(area_inund[,c(1,3)], by= "CODIMUNI") %>% 
  mutate(area_inund= ifelse(area_inund > area_resid,area_resid,
                            area_inund),
         pct_inun= round(area_inund / area_resid*100,2),
         pct_inun= ifelse(is.na(pct_inun),0,pct_inun))
```


<p style= "margin-bottom:0; font-weight:bold; font-size: 16px; color:dodgerblue;">Mapa de % d'àrea residencial en risc d'inundació amb màxim PR (500), per municipi</p>
```{r ig.height=8,fig.width=8}


tm_shape(mun_cat_inund)+
  tm_fill(col="pct_inun",
          id="NOMMUNI",
          style = "fixed",
          n=5,
          breaks = c(0,5,10,20,50,100),
          palette = c("white","darkblue"),
          title = "% d'Àrea residencial en risc <br>d'inundació (PR=500)",
          legend.format = list(text.separator= '-'),
          textNA = 'Sin datos',
          popup.vars=c("Superfície residencial del municipi: "=
                         "area_resid",
                       "% Àrea en risc d'inundació: "= "pct_inun"),
          popup.format = list(pct_inun= list(digits = 2),
                              area_resid= list(digits= 2))
  )+
  tm_borders(col="grey40",lwd=0.5) +
  tmap_options(basemaps = "Esri.WorldImagery",
               basemaps.alpha = 1) +
  tmap_mode("view") +
  tm_view(set.view = 7.5)



```


```{r}

mun_cat_inund %>% 
  select(2,14,15,17) %>%
  mutate(area_resid= round(area_resid,2),
         area_inund= round(area_inund,2),
         pct_inun= round(pct_inun,2)
         ) %>% 
  st_drop_geometry() %>% 
  head(n=15) %>% 
  DT::datatable(colnames = c("Municipi","Àrea residencial (ha)",
                           "Àrea resid. inundable (ha)", 
                           "Taxa d'inundabilitat (%)"))


```