In impact-initiatives-geospatial/surveyGEER: Bringing relevant GEE Data Into Your Assessment

knitr::opts_chunk$set(echo = TRUE)

R Markdown

Code is currently a mix of tidyrgee and easyrgee workflows. tidyrgee is under active development and this script utilizes the virtual-dataframes branch.

library(rgee)
library(tidyverse)
library(sf)
library(tidyrgee)
ee_Initialize()
devtools::load_all()
country_code <- c("ssd","eth","som")[2]
baseline_years <- c(2000:2015)

load input data

con <- DBI::dbConnect(RPostgres::Postgres(),
                      dbname = "global_gdb",
                      user      = rstudioapi::askForPassword("Database user"),
                      password      = rstudioapi::askForPassword("Database password"),
                      port     = 5432)

adm1 <- st_read(con,paste0(country_code,"_adm1"))

# make grid
# - reproject and dissolve - basically return adm0 of country 
adm1_dissolved_utm <- adm1 |> 
  reach_reproject_utm(country_code) |> 
  summarise() 

# make a grid based on the extent of the abovve
hex_grid <- adm1_dissolved_utm |> 
  st_make_grid(cellsize = 12000,square = F) |> 
  st_sf() |> 
  mutate(grid_id=row_number()) |> 
  rename(geometry=1)


# the grid now covers the bbox of the adm... 
# too make extraction faster later lets clip the grid
# lets give some buffer on the admin first to clip

# buffer
adm1_dissolved_buffered<- adm1_dissolved_utm |> 
  st_buffer(dist = 50000)

hex_grid_clipped <- hex_grid[adm1_dissolved_buffered,] 

earth engine imagery

modis_link <- "MODIS/006/MOD13Q1"
modisIC <- ee$ImageCollection(modis_link)

Rescale NDVI

modis_ndvi <- modisIC$
  select("NDVI")$
  map(
    ee_utils_pyfunc(
      function(x){x$
          multiply(0.0001)$
          copyProperties(x,x$propertyNames())
          }
    )
  )

calculate monthly stats for baseline

modis_ndvi_tidy <- as_tidyee(modis_ndvi)

monthly_baseline <- modis_ndvi_tidy |> 
  filter(year %in% baseline_years) |>
  group_by(month) |> 
  summarise(stat=list("mean","sd","median","min","max"))

composite recent monthly NDVI and combine with historical

ndvi_recent_monthly <- modis_ndvi_tidy |> 
  filter(year %in% c(2021:2022)) |>
  group_by(year,month) |> 
  summarise(
    stat="mean"
  )

ndvi_recent_renamed <- ndvi_recent_monthly |> 
  select(NDVI="NDVI_mean")

ndvi_recent_and_baseline<- inner_join(x = ndvi_recent_renamed,
                                     y = monthly_baseline,
                                     by = "month")

Calculate Indicators

• NDVI Standard Scores
• NDVI % median historical
• VCI

ndvi_recent_baseline_imageCol <- ndvi_recent_and_baseline |> 
  as_ee()

ndvi_zscore<- ndvi_recent_baseline_imageCol$map(
  function(img){
    zscore<- img$expression(
      "float((NDVI-NDVI_mean)/(NDVI_stdDev))",
      opt_map= list(NDVI= img$select("NDVI"),
                    NDVI_mean= img$select("NDVI_mean"),
                    NDVI_stdDev= img$select("NDVI_stdDev")
      )
    )$rename("NDVI_z_score")
    img$select("NDVI","NDVI_mean","NDVI_median","NDVI_min","NDVI_max")$addBands(zscore)
  }

) 
ndvi_pct_median<- ndvi_zscore$map(
  function(img){
    NDVI_pct_median<- img$expression(
      "float((NDVI)/(NDVI_median))",
      opt_map= list(NDVI= img$select("NDVI"),
                    NDVI_median= img$select("NDVI_median")

      )
    )$rename("NDVI_pct_median")
    img$select("NDVI","NDVI_min","NDVI_max","NDVI_z_score")$addBands(NDVI_pct_median)
  }

) 

vci<- ndvi_pct_median$map(
  function(img){
    vci<- img$expression(
      "float((NDVI-NDVI_min)/(NDVI_max-NDVI_min))",
      opt_map= list(NDVI= img$select("NDVI"),
                    NDVI_min= img$select("NDVI_min"),
                    NDVI_max= img$select("NDVI_max")
      )
    )$rename("VCI")
    img$select("NDVI","NDVI_pct_median","NDVI_z_score")$addBands(vci)
  }

) 

select just indicators of interest and time period of interest

ndvi_indicators <- as_tidyee(vci)

ndvi_indicators_pre_processed <- ndvi_indicators |>
  filter(year>=2021) |>
  select(
    "NDVI_z_score",
    "NDVI_pct_median",
    "VCI"
    )

Zonal Statistics

Here we aggregated to the hex grid - this could take some time, be patient.

note: We have set he via argument to "drive". You can read more about the via methods available in the rgee documentation. However, I have found "drive" to be the simplest option for large zonal reductions. However, there is a limit where GEE will send a timeout error. This limit is dynamic so there is not a defined rule of thumb yet. This being said, I was able to run the reductions on 12 km hex grids for the entire countries of SSD and ETH with no problem. Nonetheless, if you hit that error you can circumvent it but splitting up the polygon layer and looping through it... This is what you see happening with the if statements and the for-loop below. I set and arbitrary max # of grid cells to 16000. If this is exceeded the grid file is split and half and looped through then binded.

# if we have a ton of grid cells we might need to loop through them....
if(nrow(hex_grid_clipped)>16000){
  hex_grid_split <- hex_grid_clipped |> 
    mutate(
      splitter= grid_id %in% (1:(hex_grid_clipped |> nrow())/2)
    ) |> 
    group_split(splitter)
  # remamber you need a "geometry" column

  extract_grid_list <- list()
  for(i in seq_along(hex_grid_split)){
    temp_grid <- hex_grid_split[[i]]
    extract_grid_list[[i]] <- ndvi_indicators_pre_processed |> 
      ee_extract_tidy(y = temp_grid,stat = "mean",scale = 250,via = "drive")


  } 
  cat("extraction complete, binding outputs")
  grid_with_ndvi_indicators <- bind_rows(extract_grid_list)  
}

if(nrow(hex_grid_clipped)<=16000){
 grid_with_ndvi_indicators <-  ndvi_indicators_pre_processed |> 
      ee_extract_tidy(y = hex_grid_clipped,stat = "mean",scale = 250,via = "drive")
}

Pivot each indicator into it's own column

grid_with_ndvi_pcts <- grid_with_ndvi_indicators |> 
  pivot_wider(names_from = "parameter",values_from = "value") |> 
  mutate(
    NDVI_pct_median=round(NDVI_pct_median*100,1),
    VCI = round(VCI*100,1)
    )

grid_with_ndvi_yrmo <- grid_with_ndvi_pcts |> 
  mutate(yr_mo = paste0(lubridate::month(date,label=T,abbr = F)," ",lubridate::year(date))) 

grid_with_ndvi_yrmo$date |> range()


hex_grid_clipped|> 
  st_transform(crs=4326) |> 
  left_join(grid_with_ndvi_yrmo, by="grid_id") |> 
  st_write(
    "C:\\Users\\zack.arno\\Documents\\GeoCrunch\\Impact\\projects\\impact-initiatives-geospatial\\REACH_SOM\\som_hex_climate_maps.gpkg",
    paste0(country_code,"_hex_ndvi_base2015"), delete_layer=T,overwrite=T)

Set QGIS automation

create victim file

som_adm0_dissolved <- st_read(con,"som_adm1") |> 
  summarise() |> 
  mutate(uid=1)

atlas_victim <- som_adm0_dissolved |> 
  left_join(
  grid_with_ndvi |> 
  mutate(yr_mo = paste0(lubridate::month(date,label=T,abbr = F)," ",lubridate::year(date))) |> 
  distinct(date, yr_mo) |> 
  mutate(uid=1)   

  ) 

atlas_victim |> print(n=17)

atlas_victim|> 
  st_write(
    "C:\\Users\\zack.arno\\Documents\\GeoCrunch\\Impact\\projects\\impact-initiatives-geospatial\\REACH_SOM\\som_hex_climate_maps.gpkg",
    "som_atlas_coverage",delete_layer=T, overwrite=T
  )

impact-initiatives-geospatial/surveyGEER documentation built on Feb. 4, 2023, 12:13 p.m.