Forest Loss with rgee

Forest Loss with rgee

# packages
pacman::p_load('ps', "rgee", "tidyverse", "sf", "ggplot2", "patchwork",
               'reticulate')

ee_Initialize(email = '[email protected]')

# Hansen Global Forest Change
forestloss <- ee$Image("UMD/hansen/global_forest_change_2019_v1_7")
flossyear <- forestloss$select('lossyear')

# P Areas
wdpa_ee <- ee$FeatureCollection('WCMC/WDPA/current/polygons')
listpa = ee$List(c('33050', '351088', '342673'))
dulombi <- ee$Feature(wdpa_ee$filter(ee$Filter$eq("WDPA_PID", "33050")))
pas_ee <- ee$FeatureCollection(wdpa_ee$filter(ee$Filter$inList("WDPA_PID", listpa))) 
ee_national_parks <- pas_ee %>% ee_as_sf()

# Plot sf and get interactive map
plot(ee_national_parks['NAME'])

vis_lossyear = list(
  min = 0,
  max = 19,
  palette = c('yellow', 'orange', 'red', 'lightblue', 'darkblue')
)
Map$setCenter(-14.45, 11.65, 9)
Map$addLayer(flossyear,
             vis_lossyear,
             "vis_lossyear"
) +
  Map$addLayer(pas_ee, {}, "PA")

# Extract values from Hansen global_forest_change_2019_v1_7
list_losses <- list()
for(i in 1:nrow(ee_national_parks)){ #nrow(gnb_pa)
  igeom <- ee_national_parks[i, ] %>% select(WDPA_PID)
  igeom_fc <-  sf::st_make_valid(igeom) %>% sf_as_ee()
  vList <- flossyear$reduceRegion(
    reducer= ee$Reducer$toList(),
    maxPixels = 1e9,
    geometry =  igeom_fc
  )$values()$get(0)$getInfo()
  list_losses[[i]] <- tibble(yearloss = vList,
                             WDPA_PID = igeom$WDPA_PID)
}
losses_gnb_pa <- bind_rows(list_losses)

yseries <- tibble(yearloss = c(1:19),
                  data = seq.POSIXt(as.POSIXct('2001-01-01'),
                                    as.POSIXct('2019-01-01'), by = '1 year'))

parknames <- select(ee_national_parks, WDPA_PID, ORIG_NAME) %>%
  sf::st_drop_geometry()

pa_floss <- losses_gnb_pa %>% arrange(WDPA_PID, yearloss) %>%
  group_by(WDPA_PID, yearloss) %>% 
  summarise(
    n = n(),
    loss_ha = n*753.9/10000,
    .groups = 'drop') %>%
  group_by(WDPA_PID) %>%
  mutate(cumloss = cumsum(loss_ha)) %>%
  right_join(yseries, by = c('yearloss' = 'yearloss')) %>%
  left_join(parknames, by = c('WDPA_PID' = 'WDPA_PID'))

pa_floss %>% filter( WDPA_PID == 33050 )

ggplot(pa_floss) +
  geom_line(aes(x=as.Date(data), y=cumloss/1000, colour = ORIG_NAME), size = 2) +
  geom_vline(xintercept = as.Date('2013-01-01'), colour = 'white') +
  scale_x_date(breaks = '1 year', date_labels = "%Y") +
  scale_color_brewer('PA') +
  theme(
    legend.position="bottom",
    legend.text = element_text(size = 8, color = "white"),
    legend.background =  element_rect(colour = "grey20", fill = "#1e1e1e"),
    legend.title = element_text(size = 8, color = "white"),
    legend.key = element_blank(),
    strip.background =element_rect(fill="grey20"),
    strip.text = element_text(colour = 'white'),
    axis.text.y = element_text(size = 8, color = "white"),
    axis.text.x = element_text(angle = 90, vjust=.5, size = 8, color = "white"),
    axis.line = element_line(color = "grey20", size = .01),
    axis.title = element_text(size = 9, color = "white"),
    plot.title = element_text(size = 9, color = "white"),
    panel.grid.major.x = element_line(colour = "grey20",size=0.01),
    panel.grid.minor.x = element_line(colour = "grey20",size=0.01),
    panel.grid.major.y = element_line(colour = "grey20",size=0.01),
    panel.grid.minor.y = element_line(colour = "grey20",size=0.01),
    plot.subtitle = element_text(size=9, face="italic", color="white"),
    plot.caption = element_text(size=7, color="white"),
    panel.background = element_rect(colour = "grey20", fill = "#1e1e1e"),
    plot.background = element_rect(colour = "grey20", fill = "#1e1e1e")
  ) +
  labs(y="Cummulative Loss ( ha ) / 1000", x = NULL, 
       title = 'Forest loss 2001-2019', subtitle = 'No threshold applyed', colour = 'white',
       caption = '2013 onwards uses Landsat 8 OLI')