Portfolio 1: Do Parks That Are Are Strictly Protected Have Less Deforestation?

setup:

library(sf)
library(rgee)
library(s2)
library(mapview)

ee_Initialize()

## -- rgee 1.0.5 ------ earthengine-api 0.1.233 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/kudoh 
## -----------------------------------------------

options(scipen = 999)

mapviewOptions(basemaps = c("CartoDB.DarkMatter", "Esri.WorldImagery"))

gathering data:

gfc2019 = ee$Image('UMD/hansen/global_forest_change_2019_v1_7')

tc = gfc2019$select("treecover2000")$updateMask(gfc2019$select("treecover2000")$gt(10))

loss = gfc2019$select("lossyear")$eq(19)$add(tc$gt(0))$eq(2)$selfMask()

# codes: https://en.wikipedia.org/wiki/List_of_FIPS_country_codes
uv = ee$Feature(ee$FeatureCollection("USDOS/LSIB_SIMPLE/2017")$filter(ee$Filter$eq('country_co', 'UV'))$first())

filtering protected areas:

pa = ee$FeatureCollection("WCMC/WDPA/current/polygons")$filter(ee$Filter$geometry(uv$geometry()))$map(function(feature) {
  return(feature$intersection(uv$geometry())) # necessary to use return() for functions, can't just leave last variable
})

strict = pa$filter(
  ee$Filter$And(
    ee$Filter$neq("STATUS", "Proposed"),
    ee$Filter$neq("IUCN_CAT", "Not Reported"),
    ee$Filter$neq("IUCN_CAT", "Not Applicable"),
    ee$Filter$neq("IUCN_CAT", "V"),
    ee$Filter$neq("IUCN_CAT", "VI")
  )
)

lenient = pa$filter(
  ee$Filter$Or(
    ee$Filter$eq("STATUS", "Proposed"),
    ee$Filter$eq("IUCN_CAT", "Not Reported"),
    ee$Filter$eq("IUCN_CAT", "Not Applicable"),
    ee$Filter$eq("IUCN_CAT", "V"),
    ee$Filter$eq("IUCN_CAT", "VI")
  )
)

.reduceRegion():

loss_area = loss$multiply(ee$Image$pixelArea())

tc_area = tc$gt(0)$multiply(ee$Image$pixelArea())

loss_pct = function(geom) {
  reducer = function(feature, geom) {
    return(
      feature$reduceRegion(
        reducer = ee$Reducer$sum(),
        geometry = geom$geometry(),
        scale = 30,
        maxPixels = 1e9
      )$getInfo()[[1]]
    )
  }
  return(reducer(loss_area, geom) / reducer(tc_area, geom) * 100)
}

print(paste('Strictly protected parks lost',
      paste0(loss_pct(strict), "%"),
      "of the tree cover."))

## [1] "Strictly protected parks lost 0.0000507434653888927% of the tree cover."

print(paste('Less protected parks lost',
            paste0(loss_pct(lenient), "%"),
            "of the tree cover."))

## [1] "Less protected parks lost 0.132170121931095% of the tree cover."

mapmaking:

cent = uv %>% 
  ee_as_sf %>% 
  st_geometry %>% 
  st_as_s2 %>% 
  s2_centroid %>% 
  st_as_sfc() %>% 
  st_coordinates()

Map$setCenter(cent[1],cent[2],7)

blank = ee$Image()$byte()
strict_bounds = blank$paint(strict, width = 1)
lenient_bounds = blank$paint(lenient, width = 1)

hansen_burkina = Map$addLayer(
  eeObject = tc$clip(uv),
  visParams = list(palette = c("lightgrey", "#006633")),
  name = "Tree Cover in 2000 (>10% canopy closure)",
  shown = F
) +
  Map$addLayer(loss$clip(uv),
               visParams = list(palette = "red"),
               name = 'Tree Cover Loss in 2019') +
  Map$addLayer(strict_bounds,
               visParams = list(palette = "#0131C3"),
               name = "Strictly Protected Parks") + 
  Map$addLayer(lenient_bounds,
               visParams = list(palette = "#09CDF7"),
               name = "Less Protected Parks")

hansen_burkina = ee_as_mapview(hansen_burkina)