Nighttime Lights Analysis

Nighttime lights have become a commonly used resource to estimate changes in local economic activity. This document analyzes nighttime lights within the country.

Data

We use nighttime lights data from VIIRS Black Marble. Raw nighttime lights data requires correction due to cloud cover and stray light, such as lunar light. The Black Marble dataset applies advanced algorithms to correct raw nighttime light values and calibrate data so that trends in lights over time can be meaningfully analyzed. From VIIRS Black Marble, we use data from January 2012 through present—where data is available at a 500-meter resolution.

Methodology

Within different units of analysis (e.g, administrative units) we use the sum of nighttime lights. Where relevant, we distinguish between lights generated from gas flaring and lights removing gas flaring. We use the World Bank’s Global Gas Flaring Tracker which indicates the location of gas flaring locations. When removing gas flaring lights, we remove lights within 10km of a gas flaring location; when looking at lights in gas flaring locations, we take the sum of lights within 10km of gas flaring locations.

Code Setup

#### MAIN PARAMS
proj_dir     <- "/Users/rmarty/Library/CloudStorage/OneDrive-SharedLibraries-WBG/Development Data Partnership - Algeria Economic Monitor"
iso3_code    <- "DZA"
iso2_code    <- "DZ"
pc_base_year <- 2019

#### Packages
library(tidyverse)
library(sf)
library(leaflet)
library(leaflet.providers)
library(ggpubr)
library(terra)
library(tidyterra)
library(gtools)
library(readxl)
library(janitor)
library(geodata)
library(exactextractr)
library(blackmarbler)
library(dplyr)
library(readxl)
library(janitor)
library(ggpubr)
library(WDI)
library(kableExtra)
library(elevatr)
library(ggnewscale)
library(pals)
library(haven)

#### Paths

## Define
data_dir       <- file.path(proj_dir, "Data")
ntl_bm_dir     <- file.path(data_dir, "Nighttime Lights BlackMarble")
gadm_dir       <- file.path(data_dir, "GADM")
gas_flare_dir  <- file.path(data_dir, "gas-flaring")

## Create
dir.create(proj_dir)
dir.create(data_dir)
dir.create(ntl_bm_dir)
dir.create(gadm_dir)
dir.create(gas_flare_dir)

#### Bearer
nasa_bearer <- read_csv("~/Dropbox/bearer_bm.csv") %>%
  pull(token)

#### Theme
theme_manual <- theme_classic2() +
  theme(strip.background = element_blank(),
        strip.text = element_text(face = "bold"),
        plot.title = element_text(face = "bold"))

#### GADM Boundaries
for(i in 0:3){
  for(iso3_i in iso3_code){
    gadm(iso3_code, level=i, path = file.path(gadm_dir))
  }
} 

#### Load GADM Functions
load_gadm <- function(i){
  roi_sf <- file.path(gadm_dir, 'gadm') %>%
    list.files(full.names = T) %>%
    str_subset(paste0("_",i,"_pk.rds")) %>%
    map_df(function(file_path){
      readRDS(file_path) %>%
        st_as_sf()
    })
}

#### Elevation
dir.create(file.path(data_dir, "elevation"))

OUT_PATH <- file.path(data_dir, "elevation", "elev.tiff")

if(!file.exists(OUT_PATH)){
  
  
  if(length(iso3_code) > 1){
    
    # elev_list <- lapply(iso3_code, function(iso3_code_i){
    #   elevation_3s(country=iso3_code_i, path=tempdir())
    # }) 
    # elev_r <- do.call(mosaic, c(elev_list, fun="mean"))
    
    elev_list <- lapply(iso3_code, function(iso3_code_i){
      
      adm0_sf <- gadm(iso3_code, level=0, path = tempdir())
      get_elev_raster(locations = adm0_sf, z = 8)
      
    })
    elev_r <- do.call(mosaic, c(elev_list, fun="mean"))
    
  } else{
    adm0_sf <- load_gadm(0)
    
    #elev_r <- elevation_3s(country=iso3_code, path=tempdir())
    elev_r <- get_elev_raster(locations = adm0_sf, z = 8)
    
  }
  
  writeRaster(elev_r, OUT_PATH, overwrite = T)
  
}

#### NTL Months to Download
ntl_months <- seq.Date(ymd("2012-01-01"), to = ymd("2025-09-01"), by = "month") %>%
  as.character()

#### NTL Days to Download
ntl_days <- seq.Date(ymd("2021-01-01"), to = Sys.Date(), by = "day") %>%
  as.character()

#### Colors
wbg_color_light <- "#1AA1DB"
wbg_color_dark <- "#02224B"

Prep Data

Download Cities

We use a dataset from Geonames which indicates the location of cities with over 1,000 people.

dir.create(file.path(data_dir, "cities"))

## Cities
OUT_FILE_ALL <- file.path(data_dir, "cities", "cities_all.Rds")

if(!file.exists(OUT_FILE_ALL)){
  # https://public.opendatasoft.com/explore/dataset/geonames-all-cities-with-a-population-1000/export/?flg=en-us&disjunctive.cou_name_en&sort=name
  city_df <- read_delim("https://public.opendatasoft.com/api/explore/v2.1/catalog/datasets/geonames-all-cities-with-a-population-1000/exports/csv?lang=en&timezone=America%2FNew_York&use_labels=true&delimiter=%3B", delim = ";")
  
  saveRDS(city_df, OUT_FILE_ALL)
} 

## Country Cities
OUT_FILE_COUNTRY_DF      <- file.path(data_dir, "cities", "cities_country_df.Rds")
OUT_FILE_COUNTRY_SF      <- file.path(data_dir, "cities", "cities_country_sf.Rds")
OUT_FILE_COUNTRY_5K_SF   <- file.path(data_dir, "cities", "cities_country_5km_sf.Rds")
OUT_FILE_COUNTRY_10K_SF  <- file.path(data_dir, "cities", "cities_country_10km_sf.Rds")
OUT_FILE_COUNTRY_BUFF_SF <- file.path(data_dir, "cities", "cities_country_buff_sf.Rds")

if(!file.exists(OUT_FILE_COUNTRY_SF)){
  
  ## Cleanup
  city_country_df <- city_df[city_df$`Country Code` %in% iso2_code,] %>%
    clean_names()
  
  city_country_df <- city_country_df %>%
    separate(coordinates, into = c("lat", "lon"), sep = ",")
  
  ## Spatially Define
  city_country_sf <- city_country_df %>%
    st_as_sf(coords = c("lon", "lat"),
             crs = 4326)
  
  ## Buffer
  city_country_5km_sf <- city_country_sf %>%
    st_buffer(dist = 5000)
  
  city_country_10km_sf <- city_country_sf %>%
    st_buffer(dist = 10000)
  
  ## Dynamic Buffer
  city_country_sf <- city_country_sf %>%
    dplyr::mutate(buffer = case_when(
      population >= 10000000 ~ 30,
      population >= 5000000 ~ 25,
      population >= 1000000 ~ 20,
      population >= 500000 ~ 15,
      population >= 100000 ~ 10,
      TRUE ~ 5
    ))
  
  city_sub_30km_buff_sf <- city_country_sf[city_country_sf$buffer %in% 30,] %>% st_buffer(dist = 30*1000)
  city_sub_25km_buff_sf <- city_country_sf[city_country_sf$buffer %in% 25,] %>% st_buffer(dist = 25*1000)
  city_sub_20km_buff_sf <- city_country_sf[city_country_sf$buffer %in% 20,] %>% st_buffer(dist = 20*1000)
  city_sub_15km_buff_sf <- city_country_sf[city_country_sf$buffer %in% 15,] %>% st_buffer(dist = 15*1000)
  city_sub_10km_buff_sf <- city_country_sf[city_country_sf$buffer %in% 10,] %>% st_buffer(dist = 10*1000)
  city_sub_5km_buff_sf  <- city_country_sf[city_country_sf$buffer %in% 5,]  %>% st_buffer(dist = 5*1000)
  
  city_buff_sf <- bind_rows(city_sub_30km_buff_sf,
                            city_sub_25km_buff_sf,
                            city_sub_20km_buff_sf,
                            city_sub_15km_buff_sf,
                            city_sub_10km_buff_sf,
                            city_sub_5km_buff_sf)
  
  
  ## Export
  saveRDS(city_country_df, OUT_FILE_COUNTRY_DF)
  
  saveRDS(city_country_5km_sf,  OUT_FILE_COUNTRY_5K_SF)
  saveRDS(city_country_10km_sf, OUT_FILE_COUNTRY_10K_SF)
  saveRDS(city_buff_sf,         OUT_FILE_COUNTRY_BUFF_SF)
  saveRDS(city_country_sf,      OUT_FILE_COUNTRY_SF)
  
} 

Download and Append gas flaring

We download and append gas flaring data from the World Bank Global Gas Flaring Database.

dir.create(file.path(gas_flare_dir, "rawdata"))
dir.create(file.path(gas_flare_dir, "finaldata"))

#### Download data
# https://datacatalog.worldbank.org/search/dataset/0037743
if(!file.exists(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2017.xlsx"))){
  download.file(url = "https://datacatalogfiles.worldbank.org/ddh-published/0037743/DR0045623/viirs_global_flaring_d.7_slope_0.029353_2017_web_v1.xlsx?versionId=2023-01-18T20:03:32.2273754Z", 
                destfile = file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2017.xlsx"), 
                mode = "wb")
}

if(!file.exists(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2018.xlsx"))){
  download.file(url = "https://datacatalogfiles.worldbank.org/ddh-published/0037743/DR0045622/viirs_global_flaring_d.7_slope_0.029353_2018_web.xlsx?versionId=2023-01-18T20:02:43.3965005Z", 
                destfile = file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2018.xlsx"), 
                mode = "wb")
}

if(!file.exists(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2019.xlsx"))){
  download.file(url = "https://datacatalogfiles.worldbank.org/ddh-published/0037743/DR0045621/viirs_global_flaring_d.7_slope_0.029353_2019_web_v20201114-3.xlsx?versionId=2023-01-18T20:03:09.2456111Z", 
                destfile = file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2019.xlsx"), 
                mode = "wb")
}

if(!file.exists(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2020.xlsx"))){
  download.file(url = "https://datacatalogfiles.worldbank.org/ddh-published/0037743/DR0084248/2020%20Global%20Gas%20Flaring%20Volumes.xlsx?versionId=2023-01-18T20:03:53.8309309Z", 
                destfile = file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2020.xlsx"), 
                mode = "wb")
}

if(!file.exists(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2021.xlsx"))){
  download.file(url = "https://datacatalogfiles.worldbank.org/ddh-published/0037743/DR0087112/2021%20Global%20Gas%20Flaring%20Volumes.xlsx?versionId=2023-01-18T20:02:21.4951166Z", 
                destfile = file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2021.xlsx"), 
                mode = "wb")
}

#### Append data
clean_data <- function(x){
  x %>% 
    clean_names() %>% 
    dplyr::filter(iso_code %in% iso3_code)
} 

df_2021 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2021.xlsx"), 2) %>% clean_data()

df_2020_1 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2020.xlsx"), 1) %>% clean_data()
df_2020_2 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2020.xlsx"), 2) %>% clean_data()
df_2020_3 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2020.xlsx"), 3) %>% clean_data()

df_2019 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2019.xlsx"), 1) %>% clean_data()

df_2018_4 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2018.xlsx"), 4) %>% clean_data()
df_2018_5 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2018.xlsx"), 5) %>% clean_data()
df_2018_6 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2018.xlsx"), 6) %>% clean_data()

df_2017_1 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2017.xlsx"), 1) %>% clean_data()
df_2017_2 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2017.xlsx"), 2) %>% clean_data()
df_2017_3 <- read_xlsx(file.path(gas_flare_dir, "rawdata", "viirs_global_flaring_2017.xlsx"), 3) %>% clean_data()

gs_df <- bind_rows(
  df_2021,
  df_2020_1,
  df_2020_2,
  df_2020_3,
  df_2019,
  df_2018_4,
  df_2018_5,
  df_2018_6,
  df_2017_1,
  df_2017_2,
  df_2017_3
)

if(nrow(gs_df) == 0){
  
  gs_5km_sf <- gs_df
  gs_10km_sf <- gs_df
  
} else{
  
  gs_df <- gs_df %>%
    dplyr::select(latitude, longitude) %>%
    distinct() %>%
    dplyr::mutate(uid = 1:n())
  
  gs_sf <- gs_df %>% 
    st_as_sf(coords = c("longitude", "latitude"),
             crs = 4326)
  
  gs_5km_sf <- gs_sf %>%
    st_buffer(dist = 5000)
  
  gs_10km_sf <- gs_sf %>%
    st_buffer(dist = 10000)
  
}

saveRDS(gs_df, file.path(gas_flare_dir, "finaldata", "gas_flare_locations.Rds"))
saveRDS(gs_5km_sf, file.path(gas_flare_dir, "finaldata", "gas_flare_locations_5km.Rds"))
saveRDS(gs_10km_sf, file.path(gas_flare_dir, "finaldata", "gas_flare_locations_10km.Rds"))

Download Nighttime Lights

#### NTL
roi_sf <- load_gadm(0)

dir.create(file.path(ntl_bm_dir, "rasters"))
dir.create(file.path(ntl_bm_dir, "rasters", "annual"))
dir.create(file.path(ntl_bm_dir, "rasters", "monthly"))
dir.create(file.path(ntl_bm_dir, "rasters", "daily"))
dir.create(file.path(ntl_bm_dir, "rasters", "monthly_quality"))

r <- bm_raster(roi_sf = roi_sf,
               product_id = "VNP46A4",
               date = 2012:2024,
               bearer = nasa_bearer,
               output_location_type = "file",
               file_dir = file.path(ntl_bm_dir, "rasters", "annual"))

r <- bm_raster(roi_sf = roi_sf,
               product_id = "VNP46A3",
               date = ntl_months,
               bearer = nasa_bearer,
               output_location_type = "file",
               file_dir = file.path(ntl_bm_dir, "rasters", "monthly"))

r <- bm_raster(roi_sf = roi_sf,
               product_id = "VNP46A3",
               date = ntl_months,
               bearer = nasa_bearer,
               output_location_type = "file",
               variable = "NearNadir_Composite_Snow_Free_Quality",
               file_dir = file.path(ntl_bm_dir, "rasters", "monthly_quality"))

Aggregate Nighttime Lights

# Individual Files -------------------------------------------------------------
dir.create(file.path(ntl_bm_dir, "aggregated_individual_files"))

#               "adm3",
# "cities5km",
#  "cities10km",
#  "citiesbuff"
#### Setup loops
for(unit in c("adm0", 
              "adm1", 
              "adm2")){
  for(time_type in c("annual", 
                     "monthly")){
    
    if(time_type == "annual"){
      time_vec <- 2012:2024
    }
    
    if(time_type == "monthly"){
      time_vec <- file.path(ntl_bm_dir, "rasters", "monthly") %>% 
        list.files() %>%
        str_replace_all(".tif", "") %>%
        str_replace_all("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t", "") %>%
        str_replace_all("_", "-") %>%
        paste0("-01")
    }
    
    if(time_type == "daily"){
      time_vec <- file.path(ntl_bm_dir, "rasters", "daily") %>% 
        list.files() %>%
        str_replace_all(".tif", "") %>%
        str_replace_all("VNP46A2_Gap_Filled_DNB_BRDF-Corrected_NTL_qflag_t", "") %>%
        str_replace_all("_", "-")
    }
    
    for(time_i in time_vec){
      
      #### Only process if file doesn't exist
      dir.create(file.path(ntl_bm_dir, "aggregated_individual_files", unit))
      dir.create(file.path(ntl_bm_dir, "aggregated_individual_files", unit, time_type))
      
      OUT_FILE <- file.path(ntl_bm_dir, "aggregated_individual_files", unit, time_type, 
                            paste0("ntl_", time_i, ".Rds"))
      
      if(!file.exists(OUT_FILE)){
        
        #### Load Unit
        if(unit == "adm0"){
          roi_sf <- load_gadm(0)
        }
        
        if(unit == "adm1"){
          roi_sf <- load_gadm(1)
        }
        
        if(unit == "adm2"){
          roi_sf <- load_gadm(2)
        }
        
        if(unit == "cities5km"){
          roi_sf <- readRDS(file.path(data_dir, "cities", "cities_country_5km_sf.Rds"))
        }
        
        if(unit == "cities10km"){
          roi_sf <- readRDS(file.path(data_dir, "cities", "cities_country_10km_sf.Rds"))
        }
        
        if(unit == "citiesbuff"){
          roi_sf <- readRDS(file.path(data_dir, "cities", "cities_country_buff_sf.Rds"))
        }
        
        #### Load Raster
        if(time_type == "annual"){
          bm_product <- "VNP46A4"
          time_raster_i <- time_i 
        }
        
        if(time_type == "monthly"){
          bm_product <- "VNP46A3"
          
          time_raster_i <- time_i %>% 
            substring(1, 7) %>%
            str_replace_all("-", "_")
        }
        
        if(time_type == "daily"){
          bm_product <- "VNP46A2"
          
          time_raster_i <- time_i %>% 
            str_replace_all("-", "_")
        }
        
        if(time_type %in% c("annual", "monthly")){
          ntl_r <- rast(file.path(ntl_bm_dir, "rasters", time_type, 
                                  paste0(bm_product,"_NearNadir_Composite_Snow_Free_qflag_t",time_raster_i,".tif")))
        } else{
          ntl_r <- rast(file.path(ntl_bm_dir, "rasters", "daily", 
                                  paste0("VNP46A2_Gap_Filled_DNB_BRDF-Corrected_NTL_qflag_t",time_raster_i,".tif")))
        }
        
        # Replace low values with 0
        ntl_r_1na <- ntl_r
        ntl_r_1na[ntl_r_1na < 1] <- 0
        
        ntl_r_0_5na <- ntl_r
        ntl_r_0_5na[ntl_r_0_5na < 0.5] <- 0
        
        ntl_r_0_1na <- ntl_r
        ntl_r_0_1na[ntl_r_0_1na < 0.1] <- 0
        
        #### Gas Flaring Mask
        gs_5km_sf  <- readRDS(file.path(gas_flare_dir, "finaldata", "gas_flare_locations_5km.Rds"))
        gs_10km_sf <- readRDS(file.path(gas_flare_dir, "finaldata", "gas_flare_locations_10km.Rds"))
        
        if(nrow(gs_5km_sf) > 0){
          ntl_gf_5km_r   <- mask(ntl_r, gs_5km_sf)
          ntl_nogf_5km_r <- mask(ntl_r, gs_5km_sf, inverse = T)
          
          ntl_gf_10km_r   <- mask(ntl_r, gs_10km_sf)
          ntl_nogf_10km_r <- mask(ntl_r, gs_10km_sf, inverse = T)
        }
        
        #### Extract NTL
        roi_sf$ntl_sum          <- exact_extract(ntl_r, roi_sf, "sum")
        roi_sf$ntl_mean         <- exact_extract(ntl_r, roi_sf, "mean")
        
        roi_sf$ntl_sum_l1na     <- exact_extract(ntl_r_1na, roi_sf, "sum")
        roi_sf$ntl_sum_l0_5na   <- exact_extract(ntl_r_0_5na, roi_sf, "sum")
        roi_sf$ntl_sum_l0_1na   <- exact_extract(ntl_r_0_1na, roi_sf, "sum")
        
        if(time_type == "monthly"){
          qual_r <- rast(file.path(ntl_bm_dir, "rasters", "monthly_quality", 
                                   paste0(bm_product,"_NearNadir_Composite_Snow_Free_Quality_qflag_t",time_raster_i,".tif")))
          
          tmp_r   <- ntl_r
          tmp_r[] <- 1
          tmp_r   <- mask(tmp_r, roi_sf)
          
          roi_sf$n_pixel  <- exact_extract(tmp_r, roi_sf, "sum")
          roi_sf$n_qual_0 <- exact_extract(qual_r == 0, roi_sf, "sum")
          roi_sf$n_qual_1 <- exact_extract(qual_r == 1, roi_sf, "sum")
          roi_sf$n_qual_2 <- exact_extract(qual_r == 2, roi_sf, "sum")
          
          roi_sf <- roi_sf %>%
            dplyr::mutate(prop_qual_0 = n_qual_0 / n_pixel,
                          prop_qual_1 = n_qual_1 / n_pixel,
                          prop_qual_2 = n_qual_2 / n_pixel) %>%
            dplyr::select(-c(n_pixel, prop_qual_0, prop_qual_1, prop_qual_2))
          
        }
        
        if(nrow(gs_5km_sf) > 0){
          roi_sf$ntl_gf_5km_sum   <- exact_extract(ntl_gf_5km_r, roi_sf, "sum")
          roi_sf$ntl_nogf_5km_sum <- exact_extract(ntl_nogf_5km_r, roi_sf, "sum")
          
          roi_sf$ntl_gf_10km_sum   <- exact_extract(ntl_gf_10km_r, roi_sf, "sum")
          roi_sf$ntl_nogf_10km_sum <- exact_extract(ntl_nogf_10km_r, roi_sf, "sum")
          
          
          
          roi_sf$ntl_gf_5km_mean   <- exact_extract(ntl_gf_5km_r, roi_sf, "mean")
          roi_sf$ntl_nogf_5km_mean <- exact_extract(ntl_nogf_5km_r, roi_sf, "mean")
          
          roi_sf$ntl_gf_10km_mean   <- exact_extract(ntl_gf_10km_r, roi_sf, "mean")
          roi_sf$ntl_nogf_10km_mean <- exact_extract(ntl_nogf_10km_r, roi_sf, "mean")
        }
        
        #### NA Values
        r_n_obs <- exact_extract(ntl_r, roi_sf, function(values, coverage_fraction)
          sum(!is.na(values)),
          progress = F)
        
        r_n_obs_poss <- exact_extract(ntl_r, roi_sf, function(values, coverage_fraction)
          length(values),
          progress = F)
        
        roi_sf$n_pixels           <- r_n_obs_poss
        roi_sf$n_non_na_pixels    <- r_n_obs
        roi_sf$prop_non_na_pixels <- roi_sf$n_non_na_pixels / roi_sf$n_pixels 
        
        #### Cleanup
        roi_sf$date <- time_i
        
        roi_df <- roi_sf %>%
          st_drop_geometry()
        
        #### Export
        saveRDS(roi_df, OUT_FILE)
        
      }
    }
  }
}

#### Append Files
dir.create(file.path(ntl_bm_dir, "aggregated"))

for(unit in c("adm0", 
              "adm1", 
              "adm2")){
  for(time_type in c("annual", "monthly")){
    
    ntl_df <- file.path(ntl_bm_dir, "aggregated_individual_files", unit, time_type) %>%
      list.files(full.names = T,
                 pattern = ".Rds") %>%
      map_df(readRDS)
    
    saveRDS(ntl_df,   file.path(ntl_bm_dir, "aggregated", paste0("ntl_", unit, "_", time_type, ".Rds")))
    write_csv(ntl_df, file.path(ntl_bm_dir, "aggregated", paste0("ntl_", unit, "_", time_type, ".csv")))
    write_dta(ntl_df, file.path(ntl_bm_dir, "aggregated", paste0("ntl_", unit, "_", time_type, ".dta")))
    
  }
}

Maps of Nighttime Lights

Interactive Map

#| warning: false
#| message: false

## Load boundaries
adm0_sf <- load_gadm(0)

## Load/prep raster
prep_r <- function(year_i){
  r <- rast(file.path(ntl_bm_dir, "rasters", "annual",
                      paste0("VNP46A4_NearNadir_Composite_Snow_Free_qflag_t",year_i,".tif")))
  r <- r %>% mask(adm0_sf)
  r[][r[] <= 0.05] <- NA
  r[] <- log(r[] + 1)
  #r[] <- log(r[] + 1)
  return(r)
}

r_2012 <- prep_r(2012)
r_2013 <- prep_r(2013)
r_2014 <- prep_r(2014)
r_2015 <- prep_r(2015)
r_2016 <- prep_r(2016)
r_2017 <- prep_r(2017)
r_2018 <- prep_r(2018)
r_2019 <- prep_r(2019)
r_2020 <- prep_r(2020)
r_2021 <- prep_r(2021)
r_2022 <- prep_r(2022)
r_2023 <- prep_r(2023)

## Make map
pal <- colorNumeric(c("black", "yellow", "red"), unique(c(r_2012[],
                                                          r_2013[],
                                                          r_2014[],
                                                          r_2015[],
                                                          r_2016[],
                                                          r_2017[],
                                                          r_2018[],
                                                          r_2019[],
                                                          r_2020[],
                                                          r_2021[],
                                                          r_2022[],
                                                          r_2023[])),
                    na.color = "transparent")

leaflet() %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addRasterImage(r_2012, colors = pal, opacity = 1, group = "2012") %>%
  addRasterImage(r_2013, colors = pal, opacity = 1, group = "2013") %>%
  addRasterImage(r_2014, colors = pal, opacity = 1, group = "2014") %>%
  addRasterImage(r_2015, colors = pal, opacity = 1, group = "2015") %>%
  addRasterImage(r_2016, colors = pal, opacity = 1, group = "2016") %>%
  addRasterImage(r_2017, colors = pal, opacity = 1, group = "2017") %>%
  addRasterImage(r_2018, colors = pal, opacity = 1, group = "2018") %>%
  addRasterImage(r_2019, colors = pal, opacity = 1, group = "2019") %>%
  addRasterImage(r_2020, colors = pal, opacity = 1, group = "2020") %>%
  addRasterImage(r_2021, colors = pal, opacity = 1, group = "2021") %>%
  addRasterImage(r_2022, colors = pal, opacity = 1, group = "2022") %>%
  addRasterImage(r_2023, colors = pal, opacity = 1, group = "2023") %>%
  addLayersControl(
    baseGroups = paste0(2012:2023),
    options = layersControlOptions(collapsed=FALSE)
  )

Static Map of Nighttime Lights

## Load boundaries
adm0_sf <- load_gadm(0)

## Elevation
elev_r <- rast(file.path(data_dir, "elevation", "elev.tiff"))
slope_r <- terrain(elev_r)
slope_r <- slope_r %>% mask(adm0_sf)

## Latest Year
r <- rast(file.path(ntl_bm_dir, "rasters", "annual",
                    paste0("VNP46A4_NearNadir_Composite_Snow_Free_qflag_t",2023,".tif")))
r <- r %>% mask(adm0_sf)

r[][r[] <= 1] <- NA
r[] <- log(r[] + 1)
r[] <- log(r[] + 1)

ocean_colors <- ocean.deep(100)
ocean_colors[98:100] <- "#FFFFFF"

ggplot() +
  geom_spatraster(data = slope_r, 
                  show.legend = FALSE,
                  maxcell = 9e+07) +
  scale_fill_gradient(low = "#49232E",
                      high = "black",
                      na.value = "transparent") +
  ggnewscale::new_scale_fill() +
  geom_spatraster(data = r,
                  aes(fill = t2023),
                  maxcell = 9e+07) +
  scale_fill_gradientn(colors = ocean_colors,
                       #midpoint = 2,
                       na.value = "transparent") +
  labs(fill = "Light Intensity",
       title = "Nighttime Lights",
       subtitle = 2023) +
  coord_sf() +
  theme_void() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
        plot.subtitle = element_text(face = "italic", hjust = 0.5, size = 16),
        legend.position = c(0.2, 0.15)) +
  guides(fill = guide_colorbar(direction = "horizontal",
                               label = FALSE,         
                               title.position = "top"))

Diagnostics

Nighttime Lights as Proxy for GDP

# Download GDP
dir.create(file.path(data_dir, "wdi"))
OUT_FILE <- file.path(data_dir, "wdi", "gdp.Rds")

if(!file.exists(OUT_FILE)){
  gdp_df <- WDI(indicator = c("NY.GDP.PCAP.CD",
                              "NY.GDP.MKTP.CD"), 
                country = iso3_code, 
                start = 2012, 
                end = 2023)
  
  saveRDS(gdp_df, OUT_FILE)
} else{
  gdp_df <- readRDS(OUT_FILE)
}

# Correlation with GDP
ntl_df <- readRDS(file.path(ntl_bm_dir, "aggregated", 
                            paste0("ntl_", "adm0", "_", "annual", ".Rds")))

ntl_df <- ntl_df %>%
  dplyr::rename(year = date)

# Filter
if(paste(iso3_code, collapse = ";") %in% "SSD;SDN"){
  ntl_df <- ntl_df[ntl_df$COUNTRY %in% "Sudan",]
  gdp_df <- gdp_df[gdp_df$country %in% "Sudan",]
}

ntl_gdp_df <- ntl_df %>%
  left_join(gdp_df, by = "year") %>%
  clean_names()

## Line Plots
ntl_gdp_df %>%
  dplyr::select(year, 
                ntl_sum,
                #ntl_gf_10km_sum,
                #ntl_nogf_10km_sum,
                ny_gdp_pcap_cd,
                ny_gdp_mktp_cd) %>%
  pivot_longer(cols = -year) %>%
  dplyr::mutate(name = case_when(
    name == "ntl_sum" ~ "Nighttime Lights",
    name == "ntl_gf_10km_sum" ~ "Nighttime Lights:\nGas Flares",
    name == "ntl_nogf_10km_sum" ~ "Nighttime Lights:\nNo Gas Flares",
    name == "ny_gdp_pcap_cd" ~ "GDP, Per Capita",
    name == "ny_gdp_mktp_cd" ~ "GDP"
  )) %>%
  dplyr::mutate(type = name %>% str_detect("GDP")) %>%
  ggplot(aes(x = year,
             y = value,
             color = type)) +
  geom_line(linewidth = 1) +
  facet_wrap(~name,
             scales = "free_y") +
  labs(x = NULL,
       y = NULL) +
  scale_color_manual(values = c(wbg_color_light, wbg_color_dark)) +
  theme_manual +
  theme(legend.position = "none")

## Correlation: GDP
ntl_gdp_df %>%
  dplyr::select(ntl_sum, 
                #ntl_gf_10km_sum,
                #ntl_nogf_10km_sum,
                ny_gdp_mktp_cd) %>%
  pivot_longer(cols = -c(ny_gdp_mktp_cd)) %>%
  dplyr::mutate(name = case_when(
    name == "ntl_sum" ~ "NTL"
  )) %>%
  
  ggplot(aes(x = value,
             y = ny_gdp_mktp_cd)) +
  geom_point() +
  geom_smooth(method=lm, 
              se=FALSE,
              color = wbg_color_light) +
  stat_cor() +
  labs(x = "NTL",
       y = "GDP",
       title = "Correlation of Nighttime Lights with GDP") +
  facet_wrap(~name,
             scales = "free") +
  theme_manual 

## Correlation: GDP
ntl_gdp_df %>%
  dplyr::select(ntl_sum, 
                #ntl_gf_10km_sum,
                #ntl_nogf_10km_sum,
                ny_gdp_pcap_cd) %>%
  pivot_longer(cols = -c(ny_gdp_pcap_cd)) %>%
  dplyr::mutate(name = case_when(
    name == "ntl_sum" ~ "NTL"
  )) %>%
  
  ggplot(aes(x = value,
             y = ny_gdp_pcap_cd)) +
  geom_point() +
  geom_smooth(method=lm, 
              se=FALSE,
              color = wbg_color_light) +
  stat_cor() +
  labs(x = "NTL",
       y = "GDP, per Capita",
       title = "Correlation of Nighttime Lights with per Capita GDP") +
  facet_wrap(~name,
             scales = "free") +
  theme_manual

Trends in Nighttime Lights

ADM0

ntl_df <- readRDS(file.path(ntl_bm_dir, "aggregated", 
                             paste0("ntl_", "adm0", "_", "annual", ".Rds")))

ggplot() +
  geom_col(data = ntl_df,
           aes(x = date, y = ntl_sum),
           fill = wbg_color_dark) +
  labs(x = NULL,
       y = "Nighttime Lights",
       color = NULL) +
  theme_manual

ntl_df <- readRDS(file.path(ntl_bm_dir, "aggregated", 
                             paste0("ntl_", "adm0", "_", "monthly", ".Rds"))) %>%
  dplyr::mutate(date = date %>% ymd())

ggplot() +
  geom_col(data = ntl_df,
           aes(x = date, y = ntl_sum),
           fill = wbg_color_dark) +
  labs(x = NULL,
       y = "Nighttime Lights",
       color = NULL) +
  theme_manual

ADM1

ntl_df <- readRDS(file.path(ntl_bm_dir, "aggregated", 
                             paste0("ntl_", "adm1", "_", "annual", ".Rds")))

ggplot() +
  geom_col(data = ntl_df,
           aes(x = date, y = ntl_sum),
           fill = wbg_color_dark) +
  labs(x = NULL,
       y = "Nighttime Lights",
       color = NULL) +
  theme_manual +
  facet_wrap(~NAME_1,
             scales = "free_y")

ntl_df <- readRDS(file.path(ntl_bm_dir, "aggregated", 
                             paste0("ntl_", "adm1", "_", "monthly", ".Rds"))) %>%
  dplyr::mutate(date = date %>% ymd())

ggplot() +
  geom_col(data = ntl_df,
           aes(x = date, y = ntl_sum),
           fill = wbg_color_dark) +
  labs(x = NULL,
       y = "Nighttime Lights",
       color = NULL) +
  theme_manual +
  facet_wrap(~NAME_1,
             scales = "free_y")

% Change Maps

ADM 1

## Load
adm_sf <- load_gadm(1)
ntl_df <- readRDS(file.path(ntl_bm_dir, "aggregated", "ntl_adm1_annual.Rds"))

## Prep data
ntl_wide_df <- ntl_df %>%
  dplyr::filter(date %in% c(pc_base_year:2023)) %>%
  dplyr::select(GID_1, date, ntl_sum) %>%
  dplyr::mutate(date = paste0("yr", date)) %>%
  pivot_wider(names_from = date,
              values_from = ntl_sum,
              id_cols = GID_1)

ntl_wide_df$base_year <- ntl_wide_df[[paste0("yr", pc_base_year)]]

ntl_wide_df <- ntl_wide_df %>%
  dplyr::mutate(pc20 = (yr2020 - base_year)/base_year * 100,
                pc21 = (yr2021 - base_year)/base_year * 100,
                pc22 = (yr2022 - base_year)/base_year * 100,
                pc23 = (yr2023 - base_year)/base_year * 100)

## Map
adm_data_sf <- adm_sf %>%
  left_join(ntl_wide_df, by = "GID_1") %>%
  pivot_longer(cols = c(pc20, pc21, pc22, pc23)) %>%
  dplyr::mutate(name = case_when(
    name == "pc20" ~ paste0("% Change: ",pc_base_year," to 2020"),
    name == "pc21" ~ paste0("% Change: ",pc_base_year," to 2021"),
    name == "pc22" ~ paste0("% Change: ",pc_base_year," to 2022"),
    name == "pc23" ~ paste0("% Change: ",pc_base_year," to 2023")
  ))

adm_data_sf$value[adm_data_sf$value >= 100] <- 100
adm_data_sf$value[adm_data_sf$value <= -100] <- -100

ggplot() +
  geom_sf(data = adm_data_sf,
          aes(fill = value)) +
  facet_wrap(~name) +
  labs(fill = "% Change") +
  scale_fill_gradient2(low = "red",
                       mid = "white",
                       high = "forestgreen",
                       midpoint = 0,
                       limits = c(-100, 100),
                       breaks = c(-100, -50, 0, 50, 100),
                       labels = c("< -100", "-50", "0", "50", "> 100")) +
  theme_manual +
  theme_void() +
  theme(strip.text = element_text(face = "bold"))

ADM 2

## Load
adm_sf <- load_gadm(2)
ntl_df <- readRDS(file.path(ntl_bm_dir, "aggregated", "ntl_adm2_annual.Rds"))

## Prep data
ntl_wide_df <- ntl_df %>%
  dplyr::filter(date %in% c(pc_base_year:2023)) %>%
  dplyr::select(GID_2, date, ntl_sum) %>%
  dplyr::mutate(date = paste0("yr", date)) %>%
  pivot_wider(names_from = date,
              values_from = ntl_sum,
              id_cols = GID_2)

ntl_wide_df$base_year <- ntl_wide_df[[paste0("yr", pc_base_year)]]

ntl_wide_df <- ntl_wide_df %>%
  dplyr::mutate(pc20 = (yr2020 - base_year)/base_year * 100,
                pc21 = (yr2021 - base_year)/base_year * 100,
                pc22 = (yr2022 - base_year)/base_year * 100,
                pc23 = (yr2023 - base_year)/base_year * 100)

## Map
adm_data_sf <- adm_sf %>%
  left_join(ntl_wide_df, by = "GID_2") %>%
  pivot_longer(cols = c(pc20, pc21, pc22, pc23)) %>%
  dplyr::mutate(name = case_when(
    name == "pc20" ~ paste0("% Change: ",pc_base_year," to 2020"),
    name == "pc21" ~ paste0("% Change: ",pc_base_year," to 2021"),
    name == "pc22" ~ paste0("% Change: ",pc_base_year," to 2022"),
    name == "pc23" ~ paste0("% Change: ",pc_base_year," to 2023")
  ))

adm_data_sf$value[adm_data_sf$value >= 100] <- 100
adm_data_sf$value[adm_data_sf$value <= -100] <- -100

ggplot() +
  geom_sf(data = adm_data_sf,
          aes(fill = value)) +
  facet_wrap(~name) +
  labs(fill = "% Change") +
  scale_fill_gradient2(low = "red",
                       mid = "white",
                       high = "forestgreen",
                       midpoint = 0,
                       limits = c(-100, 100),
                       breaks = c(-100, -50, 0, 50, 100),
                       labels = c("< -100", "-50", "0", "50", "> 100")) +
  theme_manual +
  theme_void() +
  theme(strip.text = element_text(face = "bold"))