Skip to content

Commit

Permalink
figures by earth mi
Browse files Browse the repository at this point in the history
  • Loading branch information
ramarty committed Oct 12, 2023
1 parent 1d3fc6a commit 2438d90
Show file tree
Hide file tree
Showing 34 changed files with 777 additions and 476 deletions.
152 changes: 152 additions & 0 deletions notebooks/ntl-analysis/01_clean_data/00_prep_mi_files.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,152 @@
# Prep MI Files

# Load data --------------------------------------------------------------------
mi_sf <- read_sf(file.path(data_dir, "earthquake_intensity", "shape", "mi.shp"))

country0_sf <- read_sf(file.path(admin_bnd_dir, "gadm41_MAR_shp", "gadm41_MAR_0.shp"))
country2_sf <- read_sf(file.path(admin_bnd_dir, "gadm41_MAR_shp", "gadm41_MAR_2.shp"))
mrk_sf <- country2_sf[country2_sf$NAME_2 %in% "Marrakech",]

# Subset -----------------------------------------------------------------------
mi_strong_3_sf <- mi_sf %>%
dplyr::filter(PARAMVALUE >= 3,
PARAMVALUE < 4) %>%
mutate(id = 1) %>%
group_by(id) %>%
dplyr::summarise(geometry = st_union(geometry)) %>%
ungroup()

mi_strong_4_sf <- mi_sf %>%
dplyr::filter(PARAMVALUE >= 4,
PARAMVALUE < 5) %>%
mutate(id = 1) %>%
group_by(id) %>%
dplyr::summarise(geometry = st_union(geometry)) %>%
ungroup()

mi_strong_5_sf <- mi_sf %>%
dplyr::filter(PARAMVALUE >= 5,
PARAMVALUE < 6) %>%
mutate(id = 1) %>%
group_by(id) %>%
dplyr::summarise(geometry = st_union(geometry)) %>%
ungroup()

mi_strong_6_sf <- mi_sf %>%
dplyr::filter(PARAMVALUE >= 6,
PARAMVALUE < 7) %>%
mutate(id = 1) %>%
group_by(id) %>%
dplyr::summarise(geometry = st_union(geometry)) %>%
ungroup()

mi_strong_7_sf <- mi_sf %>%
dplyr::filter(PARAMVALUE >= 7,
PARAMVALUE < 8) %>%
mutate(id = 1) %>%
group_by(id) %>%
dplyr::summarise(geometry = st_union(geometry)) %>%
ungroup()

mi_strong_8_sf <- mi_sf %>%
dplyr::filter(PARAMVALUE >= 8,
PARAMVALUE < 9) %>%
mutate(id = 1) %>%
group_by(id) %>%
dplyr::summarise(geometry = st_union(geometry)) %>%
ungroup()

mi_strong_3to5_sf <- mi_sf %>%
dplyr::filter(PARAMVALUE >= 3,
PARAMVALUE < 6) %>%
mutate(id = 1) %>%
group_by(id) %>%
dplyr::summarise(geometry = st_union(geometry)) %>%
ungroup()

mi_strong_4to6_sf <- mi_sf %>%
dplyr::filter(PARAMVALUE >= 4,
PARAMVALUE < 7) %>%
mutate(id = 1) %>%
group_by(id) %>%
dplyr::summarise(geometry = st_union(geometry)) %>%
ungroup()

mi_strong_7to9_sf <- mi_sf %>%
dplyr::filter(PARAMVALUE >= 7,
PARAMVALUE < 9) %>%
mutate(id = 1) %>%
group_by(id) %>%
dplyr::summarise(geometry = st_union(geometry)) %>%
ungroup()

# Remove Marrekesh -------------------------------------------------------------
mi_strong_3_sf <- st_intersection(mi_strong_3_sf, country0_sf)
mi_strong_4_sf <- st_intersection(mi_strong_4_sf, country0_sf)
mi_strong_5_sf <- st_intersection(mi_strong_5_sf, country0_sf)
mi_strong_6_sf <- st_intersection(mi_strong_6_sf, country0_sf)
mi_strong_7_sf <- st_intersection(mi_strong_7_sf, country0_sf)
mi_strong_8_sf <- st_intersection(mi_strong_8_sf, country0_sf)
mi_strong_3to5_sf <- st_intersection(mi_strong_3to5_sf, country0_sf)
mi_strong_4to6_sf <- st_intersection(mi_strong_4to6_sf, country0_sf)
mi_strong_7to9_sf <- st_intersection(mi_strong_7to9_sf, country0_sf)

if(st_intersects(mi_strong_3_sf, mrk_sf, sparse = F)[1]){
mi_strong_3_sf <- st_difference(mi_strong_3_sf, mrk_sf)
}

if(st_intersects(mi_strong_4_sf, mrk_sf, sparse = F)[1]){
mi_strong_4_sf <- st_difference(mi_strong_4_sf, mrk_sf)
}

if(st_intersects(mi_strong_5_sf, mrk_sf, sparse = F)[1]){
mi_strong_5_sf <- st_difference(mi_strong_5_sf, mrk_sf)
}

if(st_intersects(mi_strong_6_sf, mrk_sf, sparse = F)[1]){
mi_strong_6_sf <- st_difference(mi_strong_6_sf, mrk_sf)
}

if(st_intersects(mi_strong_7_sf, mrk_sf, sparse = F)[1]){
mi_strong_7_sf <- st_difference(mi_strong_7_sf, mrk_sf)
}

if(st_intersects(mi_strong_8_sf, mrk_sf, sparse = F)[1]){
mi_strong_8_sf <- st_difference(mi_strong_8_sf, mrk_sf)
}

if(st_intersects(mi_strong_3to5_sf, mrk_sf, sparse = F)[1]){
mi_strong_3to5_sf <- st_difference(mi_strong_3to5_sf, mrk_sf)
}

if(st_intersects(mi_strong_4to6_sf, mrk_sf, sparse = F)[1]){
mi_strong_4to6_sf <- st_difference(mi_strong_4to6_sf, mrk_sf)
}

if(st_intersects(mi_strong_7to9_sf, mrk_sf, sparse = F)[1]){
mi_strong_7to9_sf <- st_difference(mi_strong_7to9_sf, mrk_sf)
}

# Export -----------------------------------------------------------------------
OUT_PATH <- file.path(data_dir, "earthquake_intensity", "separate_files_by_intensity")

mi_indv_sf <- bind_rows(
mi_strong_3_sf %>% mutate(mi = 3),
mi_strong_4_sf %>% mutate(mi = 4),
mi_strong_5_sf %>% mutate(mi = 5),
mi_strong_6_sf %>% mutate(mi = 6),
mi_strong_7_sf %>% mutate(mi = 7),
mi_strong_8_sf %>% mutate(mi = 8)
)

write_sf(mi_indv_sf, file.path(OUT_PATH, "mi_indiv.geojson"), delete_dsn = T)

write_sf(mi_strong_3to5_sf, file.path(OUT_PATH, "mi_3to5.geojson"), delete_dsn = T)
write_sf(mi_strong_4to6_sf, file.path(OUT_PATH, "mi_4to6.geojson"), delete_dsn = T)
write_sf(mi_strong_7to9_sf, file.path(OUT_PATH, "mi_7to9.geojson"), delete_dsn = T)






Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ bearer <- file.path(ntl_dir, "blackmarble-bearer-token", "bearer_bm.csv") %>%
read_csv() %>%
pull("token")

# Load Lebanon boundaries
# Load boundaries
roi_sf <- read_sf(file.path(admin_bnd_dir, "gadm41_MAR_shp", "gadm41_MAR_0.shp"))

# Download data ----------------------------------------------------------------
Expand All @@ -28,7 +28,7 @@ bm_raster(roi_sf = roi_sf,

bm_raster(roi_sf = roi_sf,
product_id = "VNP46A2",
date = seq.Date(from = ymd("2023-08-01"), to = Sys.Date(), by = 1),
date = seq.Date(from = ymd("2023-08-01"), to = Sys.Date(), by = 1) %>% rev(),
bearer = bearer,
quality_flag_rm = c(255,2),
output_location_type = "file",
Expand Down
33 changes: 32 additions & 1 deletion notebooks/ntl-analysis/01_clean_data/02_extract_to_polygons.R
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,11 @@
# -- (c) Simulated DMSP (from VIIRS) [2014 - 2021]
# 3. Aggregate monthly data, VIIRS Black Marble only [2012 - present]

for(roi_name in c("adm0", "adm1", "adm2", "adm3", "adm4")){
OUT_PATH <- file.path(data_dir, "earthquake_intensity", "sepate_files_by_intensity")

for(roi_name in c("adm0", "adm1", "adm2", "adm3", "adm4",
"mi_indiv",
"mi_3to5", "mi_4to6", "mi_7to9")){

# Make Directories -------------------------------------------------------------
dir.create(file.path(ntl_dir, "aggregated-to-polygons", roi_name))
Expand Down Expand Up @@ -37,6 +41,11 @@ for(roi_name in c("adm0", "adm1", "adm2", "adm3", "adm4")){
roi_sf <- read_sf(file.path(admin_bnd_dir, "gadm41_MAR_shp", "gadm41_MAR_4.shp"))
}

if(str_detect(roi_name, "mi")){
roi_sf <- read_sf(file.path(data_dir, "earthquake_intensity", "separate_files_by_intensity",
paste0(roi_name, ".geojson")))
}

roi_sf$adm_id <- 1:nrow(roi_sf)

# Aggregate annual -------------------------------------------------------------
Expand All @@ -50,6 +59,14 @@ for(roi_name in c("adm0", "adm1", "adm2", "adm3", "adm4")){
bm_r <- raster(file.path(ntl_dir, "ntl-rasters", "blackmarble", "annual", paste0("VNP46A4_t",year,".tif")))

roi_sf$ntl_bm_mean <- exact_extract(bm_r, roi_sf, 'mean')
roi_sf$ntl_bm_q50 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.5)
roi_sf$ntl_bm_q60 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.6)
roi_sf$ntl_bm_q70 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.7)
roi_sf$ntl_bm_q80 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.8)
roi_sf$ntl_bm_q90 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.9)
roi_sf$ntl_bm_q95 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.95)
roi_sf$ntl_bm_q99 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.99)

roi_sf$year <- year

roi_df <- roi_sf
Expand Down Expand Up @@ -89,6 +106,13 @@ for(roi_name in c("adm0", "adm1", "adm2", "adm3", "adm4")){
bm_r <- raster(file.path(ntl_dir, "ntl-rasters", "blackmarble", "monthly", r_month_i))

roi_sf$ntl_bm_mean <- exact_extract(bm_r, roi_sf, 'mean')
roi_sf$ntl_bm_q50 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.5)
roi_sf$ntl_bm_q60 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.6)
roi_sf$ntl_bm_q70 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.7)
roi_sf$ntl_bm_q80 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.8)
roi_sf$ntl_bm_q90 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.9)
roi_sf$ntl_bm_q95 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.95)
roi_sf$ntl_bm_q99 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.99)

roi_sf$date <- month_i %>%
str_replace_all("_", "-") %>%
Expand Down Expand Up @@ -128,6 +152,13 @@ for(roi_name in c("adm0", "adm1", "adm2", "adm3", "adm4")){
bm_r <- raster(file.path(ntl_dir, "ntl-rasters", "blackmarble", "daily", r_day_i))

roi_sf$ntl_bm_mean <- exact_extract(bm_r, roi_sf, 'mean')
roi_sf$ntl_bm_q50 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.5)
roi_sf$ntl_bm_q60 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.6)
roi_sf$ntl_bm_q70 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.7)
roi_sf$ntl_bm_q80 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.8)
roi_sf$ntl_bm_q90 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.9)
roi_sf$ntl_bm_q95 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.95)
roi_sf$ntl_bm_q99 <- exact_extract(bm_r, roi_sf, 'quantile', quantiles = 0.99)

roi_sf$date <- day_i %>%
str_replace_all("_", "-") %>%
Expand Down
29 changes: 22 additions & 7 deletions notebooks/ntl-analysis/02_analysis/adm_tables.R
Original file line number Diff line number Diff line change
Expand Up @@ -83,9 +83,13 @@ for(roi_name in c("adm2", "adm3")){
# Export ---------------------------------------------------------------------
ntl_df <- ntl_df %>%
dplyr::select(-"2022") %>%
dplyr::mutate(pc_14_23pre_eq = (pre_2023 - `2014`) / `2014`,
pc_14_23post_eq = (post_2023 - `2014`) / `2014`) %>%
dplyr::rename(yr2014 = "2014",
yr2023_pre_eq = "pre_2023",
yr2023_post_eq = "post_2023")
ntl_df$pc_14_23pre_eq[ntl_df$pc_14_23pre_eq == Inf] <- NA
ntl_df$pc_14_23post_eq[ntl_df$pc_14_23post_eq == Inf] <- NA

write_csv(ntl_df, file.path(tables_dir,
"ntl_adm2_adm3",
Expand All @@ -94,19 +98,24 @@ for(roi_name in c("adm2", "adm3")){
# Make png table -------------------------------------------------------------
if(roi_name == "adm2"){
ntl_df %>%
dplyr::select(NAME_1, NAME_2, eq_mi, yr2014, yr2023_pre_eq, yr2023_post_eq) %>%
dplyr::select(NAME_1, NAME_2, eq_mi, yr2014, yr2023_pre_eq, yr2023_post_eq, pc_14_23pre_eq, pc_14_23post_eq) %>%
arrange(-eq_mi) %>%
head(30) %>%
dplyr::mutate(yr2014 = round(yr2014, 2),
#yr2022 = round(yr2022, 2),
yr2023_pre_eq = round(yr2023_pre_eq, 2),
yr2023_post_eq = round(yr2023_post_eq, 2)) %>%
yr2023_post_eq = round(yr2023_post_eq, 2),

pc_14_23pre_eq = paste0(round(pc_14_23pre_eq, 2),"% ."),
pc_14_23post_eq = paste0(round(pc_14_23post_eq, 2), "% .")) %>%
dplyr::rename("ADM 1" = NAME_1,
"ADM 2" = NAME_2,
"EQ: MI" = eq_mi,
"2014" = yr2014,
"2023, Pre Eq." = yr2023_pre_eq,
"2023, Post Eq." = yr2023_post_eq) %>%
"2023, Post Eq." = yr2023_post_eq,
"PC 2014 to 2023, Pre Eq." = pc_14_23pre_eq,
"PC 2014 to 2023, Post Eq." = pc_14_23post_eq) %>%
gt() %>%
# gt_color_rows(c(`PC: 19 to 22`,
# `PC: 20 to 22`,
Expand All @@ -118,27 +127,33 @@ for(roi_name in c("adm2", "adm3")){

if(roi_name == "adm3"){
ntl_df %>%
dplyr::select(NAME_1, NAME_2, NAME_3, eq_mi, yr2014, yr2023_pre_eq, yr2023_post_eq) %>%
dplyr::select(NAME_1, NAME_2, NAME_3, eq_mi, yr2014, yr2023_pre_eq, yr2023_post_eq, pc_14_23pre_eq, pc_14_23post_eq) %>%
arrange(-eq_mi) %>%
head(30) %>%
dplyr::mutate(yr2014 = round(yr2014, 2),
#yr2022 = round(yr2022, 2),
yr2023_pre_eq = round(yr2023_pre_eq, 2),
yr2023_post_eq = round(yr2023_post_eq, 2)) %>%
yr2023_post_eq = round(yr2023_post_eq, 2),

pc_14_23pre_eq = paste0(round(pc_14_23pre_eq, 2),"% ."),
pc_14_23post_eq = paste0(round(pc_14_23post_eq, 2), "% .")) %>%
dplyr::rename("ADM 1" = NAME_1,
"ADM 2" = NAME_2,
"ADM 3" = NAME_3,
"EQ: MI" = eq_mi,
"2014" = yr2014,
"2023, Pre Eq." = yr2023_pre_eq,
"2023, Post Eq." = yr2023_post_eq) %>%
"2023, Post Eq." = yr2023_post_eq,
"PC 2014 to 2023, Pre Eq." = pc_14_23pre_eq,
"PC 2014 to 2023, Post Eq." = pc_14_23post_eq) %>%
gt() %>%
# gt_color_rows(c(`PC: 19 to 22`,
# `PC: 20 to 22`,
# `PC: 21 to 22`),
# palette = "Reds",
# reverse = T) %>%
gtsave(filename = file.path(figures_dir, paste0("eq_table_",roi_name,".png")))
gtsave(filename = file.path(figures_dir,
paste0("eq_table_",roi_name,".png")))
}
}

Expand Down
63 changes: 63 additions & 0 deletions notebooks/ntl-analysis/02_analysis/daily_trends.R
Original file line number Diff line number Diff line change
Expand Up @@ -41,3 +41,66 @@ ntl_df %>%

ggsave(filename = file.path(figures_dir, "daily_trends_adm1.png"),
height = 5, width = 10)

# MI Zones ---------------------------------------------------------------------
ntl_indiv_df <- readRDS(file.path(ntl_dir, "aggregated-to-polygons", "mi_indiv",
"mi_indiv_daily_ntl.Rds"))

ntl_3to5_df <- readRDS(file.path(ntl_dir, "aggregated-to-polygons", "mi_3to5",
"mi_3to5_daily_ntl.Rds")) %>%
mutate(mi = "MI = 3 to 5")

ntl_4to6_df <- readRDS(file.path(ntl_dir, "aggregated-to-polygons", "mi_4to6",
"mi_4to6_daily_ntl.Rds")) %>%
mutate(mi = "MI = 4 to 6")

ntl_7to9_df <- readRDS(file.path(ntl_dir, "aggregated-to-polygons", "mi_7to9",
"mi_7to9_daily_ntl.Rds")) %>%
mutate(mi = "MI = 7 to 9")

ntl_range_df <- bind_rows(ntl_3to5_df,
ntl_4to6_df,
ntl_7to9_df)

ntl_indiv_df %>%
mutate(mi = paste0("MI = ", mi)) %>%
ggplot() +
geom_col(aes(x = date,
y = ntl_bm_mean),
fill = "gray30") +
geom_vline(xintercept = ymd("2023-09-08"),
color = "red") +
labs(x = NULL,
y = "Nighttime Lights",
title = "Daily Nighttime Lights: by Earthquake Intensity",
subtitle = "Marrakech excluded from analysis") +
theme_classic2() +
facet_wrap(~mi) +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(size = 6))

ggsave(filename = file.path(figures_dir, "daily_trends_eq_indiv.png"),
height = 3, width = 5.5)



ntl_range_df %>%
ggplot() +
geom_col(aes(x = date,
y = ntl_bm_mean),
fill = "gray30") +
geom_vline(xintercept = ymd("2023-09-08"),
color = "red") +
labs(x = NULL,
y = "Nighttime Lights",
title = "Daily Nighttime Lights: by Earthquake Intensity",
subtitle = "Marrakech excluded from analysis") +
theme_classic2() +
facet_wrap(~mi) +
theme(strip.background = element_blank(),
strip.text = element_text(face = "bold"),
axis.text.x = element_text(size = 6))

ggsave(filename = file.path(figures_dir, "daily_trends_eq_range.png"),
height = 2.5, width = 5)
Loading

0 comments on commit 2438d90

Please sign in to comment.