03-donnees_deforestation.qmd

---
output: html_document
editor_options: 
  chunk_output_type: console
---

# Couvert forestier {#sec-deforestation}

## Carvalho et al. (source MODIS)

Pour commencer, on récupère le travail réalisé par Carvalho et al. [-@carvalho_methods_2020] qui complète les informations physiques de Goodman et al. [-@goodman_les_2018] avec des données relatives au couvert forestier en 1996, 2006 et 2016 et la diversité d'espèces présentes.

```{r}
library(tidyverse)
library(readxl)
library(sf)
library(wdpar)
library(mapme.biodiversity)

# Voir le chapitre "Fondamentaux R" pour une aide à l'import.
sup2 <- read_xlsx("data/Carvalho2018sup2.xlsx", # Enplacement du fichier
                  skip = 8, # Premières lignes du tableau excel à ne pas lire
                  n_max = 101,  # on ne lit pas les dernières lignes (notes)
                  col_types = c("text", "text", "text", "text", 
        "numeric", "numeric", "numeric", "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", "numeric"))
sup4 <- read_xlsx("data/Carvalho2018sup4.xlsx", skip = 6,
                  col_types = c("text", "numeric", "text", "numeric", "numeric", 
        "numeric", "numeric", "numeric", "numeric", "numeric"))
        

# Carvalho et al. 2008 document in their supp. material 2: "The three parcels that made up
# Andohahela (Parcels I, II and III) comprised different types of dominant vegetation and
# associated animal species, and were exposed to distinct pressures. Andohahela was analysed
# in its entirety (site number 57), as well as separated"

sup2 <- sup2 %>% 
  mutate(PA = recode(PA, `Andohahela complete` = "Andohahela"),
         num_atlas_ = as.integer(`Site number`))

sup4 <- sup4 %>%
  filter(`Habitat type` == "TOTAL") %>%
  mutate(num_atlas_ = as.numeric(Parcel))


load("data/ch2_AP_Vahatra.rds")

AP_Vahatra <- AP_Vahatra %>%
  left_join(sup2, by = "num_atlas_") %>%
  relocate(PA, .after = nom) %>%
  left_join(sup4, by = "num_atlas_")

```

On complète cette information avec des données de couvert forestier.

## Mapme (exemple GFC)

La procédure de traitement de ces fichiers sur Mapme est analogue à celle employée dans la section @sec-caracteristiques.

```{r}
#| eval: false

# On charge les polygones travaillés au chapitre 1
WDPA_Mada <- wdpa_read("data/WDPA/WDPA_Oct2022_MDG-shapefile.zip") 

# On charge aussi le contour de Madagascar
load("data/contour_mada.rds")

# On charge les polygones travaillés au chapitre 1
WDPA_poly  <- WDPA_Mada %>% # On enlève les AP pour lesquelles on n'a que des points
  filter(st_geometry_type(.) == "MULTIPOLYGON") %>% 
  st_cast("POLYGON")

WDPA_poly <- init_portfolio(WDPA_poly,
                           years = 2000:2020,
                           outdir  = "out",
                           cores = 18,
                           add_resources = TRUE)


# Get GFW data
WDPA_poly  <- get_resources(x = WDPA_poly , 
                         resources = c("gfw_treecover", "gfw_lossyear"))

# Indicateurs de couvert forestier
WDPA_poly  <- calc_indicators(x = WDPA_poly,
                            indicators = "treecover_area", 
                            min_cover = 10, min_size = 1)

deforest_par_an <- WDPA_poly %>%
  unnest(treecover_area) %>%
  # filter(!is.na(years)) %>%
  pivot_wider(names_from = "years", values_from = "treecover", 
              names_prefix = "treecover_") %>%
  st_drop_geometry() %>%
  select(-assetid) %>%
  group_by(across(WDPAID:CONS_OBJ)) %>%
  summarise(across(starts_with("treecover"), sum, na.rm = TRUE))

# Stats pour AP_Vahatra ------------------------------------

# On charge les polygones travaillés au chapitre 2
load("data/ch2_AP_Vahatra.rds")
APV_poly  <-AP_Vahatra %>% # On enlève les AP pour lesquelles on n'a que des points
  filter(st_geometry_type(.) == "MULTIPOLYGON") %>% 
  st_make_valid() %>% 
  st_cast("POLYGON") %>% 
  st_make_valid()

APV_poly <- init_portfolio(APV_poly,
                           years = 2000:2020,
                           outdir  = "out_Vahatra",
                           cores = 18,
                           add_resources = TRUE)


# Get GFW data
APV_poly  <- get_resources(x = APV_poly , 
                         resources = c("gfw_treecover", "gfw_lossyear"))

# Indicateurs de couvert forestier
APV_poly  <- calc_indicators(x = APV_poly,
                            indicators = "treecover_area", 
                            min_cover = 10, min_size = 1)

APV_par_aire <- APV_poly %>%
  unnest(treecover_area) %>%
  # filter(!is.na(years)) %>%
  pivot_wider(names_from = "years", values_from = "treecover", 
              names_prefix = "treecover_") %>%
  st_drop_geometry() %>%
  select(-assetid) %>%
  group_by(across(!starts_with("treecover"))) %>%
  summarise(across(starts_with("treecover"), sum, na.rm = TRUE))


# Stats pour Mada --------------------------------------------------------

# On charge les polygones travaillés au chapitre 1
mada_poly  <- contour_mada %>% # On enlève les AP pour lesquelles on n'a que des points
  filter(st_geometry_type(.) == "MULTIPOLYGON") %>% 
  st_cast("POLYGON")

mada_poly <- init_portfolio(mada_poly ,
                           years = 2000:2020,
                           outdir  = "out_Mada",
                           cores = 18,
                           add_resources = TRUE)


# Get GFW data
mada_poly  <- get_resources(x = mada_poly , 
                         resources = c("gfw_treecover", "gfw_lossyear"))

# Indicateurs de couvert forestier
mada_poly  <- calc_indicators(x = mada_poly,
                            indicators = "treecover_area", 
                            min_cover = 10, min_size = 1)

mada_global <- mada_poly %>%
  unnest(treecover_area) %>%
  # filter(!is.na(years)) %>%
  pivot_wider(names_from = "years", values_from = "treecover", 
              names_prefix = "treecover_") %>%
  st_drop_geometry() %>%
  select(-assetid) %>%
  group_by(across(!starts_with("treecover"))) %>%
  summarise(across(starts_with("treecover"), sum, na.rm = TRUE))


library(writexl)
write_xlsx(list(WDPA = deforest_par_an,
             Vahatra = APV_par_aire, 
             Madagascar = mada_global),
        path = "couvert_forestier_10_1.xlsx")



```

Toutefois, en raison d'un problème liés à la gestion des calculs volumineux, les calculs pour certaines aires protégées renvoient des données aberrantes. Ce point sera mis à jour dans ce guide dès la résolution des erreurs rencontrées. 

A ce stade on se concentrera donc sur les données de TFM présentées plus bas.


## Google Earth Engine (exemple GFC)

La plateforme Google Earth Engine est un outil particulièrement pratique et performant pour mobiliser et traiter des données satellitaires. Google Earth Engine peut être utilisé :

-   en interrogeant son API, et notamment :
    -   en python, avec la librairie `gee` permet d'interroger l'API de Google Earth Engine.
    -   en R, au travers de la librairie `rgee`. Cette dernière est relativement facile d'usage, mais elle est difficile à configurer. Pour aller plus loin : https://r-earthengine.com/rgeebook/
-   directement sur la plateforme https://code.earthengine.google.com/

La consolde de codage de Google Earth Engine prend la forme suivante :

[![Diagramme des composants de la console Google Earth Engine](figs/Code_editor_diagram.png){fig-alt="Une capture annotée de l'interface de google earth engine"}](https://developers.google.com/earth-engine/guides/playground)

Le langage utilisé sur cet interface est du Javascript. Ci-dessous un exemple de code qui génère les surface (en hectares) de perte de couvert forestier. Pour fonctionner, ce code doit être collé dans un script sur la plateforme Google Earth Engine lancé en cliquant sur "Run", puis en cliquant sur "Tasks" pour exécuter le code.

```{js}
#| eval: false

 scale = 30
	
	// PREPARE DATA
	//look at tree cover, find the area
	var treeCover = gfc2021.select(['treecover2000']);
	var areaCover = treeCover.multiply(ee.Image.pixelArea())
					.divide(10000).select([0],["areacover"])
	// total loss area
	var loss = gfc2021.select(['loss']);
	var areaLoss = loss.gt(0).multiply(ee.Image.pixelArea())
				   .divide(10000).select([0],["arealoss"]);
	// total gain area
	var gain = gfc2021.select(['gain'])
	var areaGain = gain.gt(0).multiply(ee.Image.pixelArea())
				   .divide(10000).select([0],["areagain"]);
	// final image
	var total = gfc2021.addBands(areaCover)
				.addBands(areaLoss)
				.addBands(areaGain)

	// TOTAL COVER
	// Map cover area per feature
	var districtSums = areaCover.reduceRegions({
	  collection: testgu,
	  reducer: ee.Reducer.sum(),
	  scale: scale,
	});			
	
	var addVar = function(feature) {

	  // function to iterate over the sequence of years
	  var addVarYear = function(year, feat) {
		// cast var
		year = ee.Number(year).toInt()
		feat = ee.Feature(feat)

		// actual year to write as property
		var actual_year = ee.Number(2000).add(year)

		// filter year:
		// 1st: get mask
		var filtered = total.select("lossyear").eq(year)
		// 2nd: apply mask
		filtered = total.updateMask(filtered)

		// reduce variables over the feature
		var reduc = filtered.reduceRegion({
		  geometry: feature.geometry(),
		  reducer: ee.Reducer.sum(),
		  scale: scale,
		  maxPixels: 1e9
		})

		// get results
		var loss = ee.Number(reduc.get("arealoss"))
		var gain = ee.Number(reduc.get("areagain"))

		// set names
		var nameloss = ee.String("loss_").cat(actual_year)
		var namegain = ee.String("gain_").cat(actual_year)

		// set properties to the feature
		return feat.set(nameloss, loss, namegain, gain)
	  }

	  // iterate over the sequence
	  var newfeat = ee.Feature(years.iterate(addVarYear, feature));

	  // return feature with new properties
	  return newfeat
	}

	// Map over the FeatureCollection
	var areas = districtSums.map(addVar);
	
	// Export PA deforestation to a CSV file.
	Export.table.toDrive({
	  collection: areas,
	  description: 'forest_loss_WDPA_Madagascar',
	  fileFormat: 'CSV'
	});

```

## Python (exemple TMF)

Fichiers préparés en python (code à venir), directement sur les rasters.

```{r}
#| eval: false
# On charge les fichiers préparés par Marc en python (contient 2 feuilles)
tableur_tmf <- "data/TMFchangeYear_AP_Vahatra.xlsx"
# On commence par charger la feuille déforestation
tmf_vahatra_defor <- read_excel(tableur_tmf,
                              sheet = "TMFdeforestationYear") %>%
  select(nom, starts_with("TMF")) # on ne garde que les variables pertinentes
# On fait ensuite de même avec la feuille dégradation
tmf_vahatra_degrad <- read_excel(tableur_tmf,
                              sheet = "TMFdegradationYear") %>%
  select(nom, starts_with("TMF")) # Onn ne garde que les feuilles pertinentes

AP_Vahatra <- AP_Vahatra %>%
  left_join(tmf_vahatra_degrad, by = "nom") %>%
  left_join(tmf_vahatra_defor, by = "nom") 


TMF_ratios <- AP_Vahatra %>%
  st_drop_geometry() %>%
  select(nom, starts_with("Forest cover"), starts_with("TMF")) %>%
  pivot_longer(cols = starts_with("TMF"), 
               names_to = "variable", 
               values_to = "surface_ha") %>%
  mutate(an_valeur = str_extract(variable, "[:digit:]{4}"),
         an_valeur = as.numeric(an_valeur),
         surface_ratio = case_when(
           an_valeur < 2000 ~ surface_ha / `Forest cover (ha) in 2006`,
           an_valeur > 2009 ~ surface_ha / `Forest cover (ha) in 2016`,
           TRUE ~ surface_ha / `Forest cover (ha) in 2006`),
         variable = str_replace(variable, "HA", "ratio")) %>%
  select(nom, variable, surface_ratio) %>%
  pivot_wider(names_from = variable, values_from = surface_ratio) %>%
  select(nom, starts_with("TMF"))

AP_Vahatra <- AP_Vahatra %>%
  left_join(TMF_ratios, by = "nom")

save(AP_Vahatra, file = "data/ch3_AP_Vahatra.rds")
```

## Alternatives

Si on n'est pas à l'aise avec les outils mentionnés plus haut, [l'outil Geoquery d'AidData](http://geo.aiddata.org/query/ "Geoquery : un outil simple pour télécharger des statistiques élaborées à partir d'image satellitaires pour des aires administratives.") permet d'obtenir des statistiques par aire administrative. Il est également possible de formuler des demandes spécifiques pour d'autres polygones que des aires administratives [au travers d'un formulaire dédié](https://www.aiddata.org/geoquery/custom-requests "Un formulaire pour demander à AidData de réaliser des calculs sur des polygones autres que des aires administratives").