Nighttime Lights Analysis

Nighttime lights have become a commonly used resource to estimate changes in local economic activity. This section shows where nighttime lights are concentrated across Syria and observed changes over time.

Data

We use nighttime lights data from the VIIRS Black Marble dataset. 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.

The following code downloads and processes data: _main.R Defines filepaths and loads packages 01_download_ntl.R Downloads NTL data 02_extract_data.R Extracts/aggregates nighttime light data to different polygons

Map of nighttime lights

We first show a map of nighttime lights. Most of the country is dark, with lights concentrated within cities.

Code

## Load boundaries
adm0_sf <- read_sf(file.path(boundaries_dir, 
                             "ner_adm_ignn_20230720_ab_shp",
                             "NER_admbnda_adm0_IGNN_20230720.shp"))

## Load/prep raster
prep_r <- function(year_i){
  r <- rast(file.path(ntl_dir, "individual_rasters", "annually",
                      paste0("VNP46A4_NearNadir_Composite_Snow_Free_qflag_t",year_i,".tif")))
  r <- r %>% mask(adm0_sf)
  r[][r[] == 0] <- 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("yellow", "orange", "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)
  )

Association of NTL with GDP

Existing research demonstrates that nighttime lights correlates with GDP, serving as a proxy for economic activity (for example, see here). Below we show the association between nighttime lights and GDP for Niger at the annual and country level.

Code

## NTL
adm0_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm0_annually.Rds"))

## GDP
gdp_df <- WDI(country = "NER",
              indicator = c("NY.GDP.MKTP.KD", "NY.GDP.PCAP.KD"),
              start = 1992,
              end = 2023,
              extra = TRUE)

gdp_long_df <- gdp_df %>%
  dplyr::select(year, NY.GDP.MKTP.KD, NY.GDP.PCAP.KD) %>%
  dplyr::rename("GDP" = "NY.GDP.MKTP.KD",
                "GDP, Per Capita" = "NY.GDP.PCAP.KD") %>%
  pivot_longer(cols = -year) %>%
  dplyr::rename(date = year) %>%
  clean_names()

## Annual trend figures
data_long_df <- bind_rows(
  gdp_long_df %>%
    dplyr::select(date, value, name),
  
  adm0_df %>%
    dplyr::mutate(value = ntl_sum,
                  name = "Nighttime Lights") %>%
    dplyr::select(date, value, name)
)

data_long_df %>%
  dplyr::filter(date >= 2012) %>%
  ggplot(aes(x = date,
             y = value)) +
  geom_line() +
  facet_wrap(~name,
             scales = "free_y",
             ncol = 1) +
  theme_classic2() +
  theme_custom +
  labs(x = "Year",
       y = NULL)

Code

## Correlation Figures
data_wide_df <- adm0_df %>%
  left_join(gdp_long_df, by = "date")

data_wide_df %>%
  dplyr::filter(name == "GDP") %>%
  ggplot(aes(x = value,
             y = ntl_sum)) +
  geom_point() +
  geom_smooth(se = F,
              method = "lm",
              color = "darkorange") +
  stat_cor() +
  theme_classic2() +
  theme_custom +
  labs(x = "GDP",
       y = "Nighttime lights",
       title = "GDP vs Nighttime Lights")

Code

data_wide_df %>%
  dplyr::filter(name == "GDP, Per Capita") %>%
  ggplot(aes(x = value,
             y = ntl_sum)) +
  geom_point() +
  geom_smooth(se = F,
              method = "lm",
              color = "darkorange") +
  stat_cor() +
  theme_classic2() +
  theme_custom +
  labs(x = "GDP, Per Capita",
       y = "Nighttime lights",
       title = "Per Capita GDP vs Nighttime Lights")

Trends in nighttime lights within past 5 years

Below we show trends in nighttime lights from January 2019 through present at the monthly level

Code

adm0_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm0_monthly.Rds"))

adm0_df %>%
  filter(date >= ymd("2019-01-01")) %>%
  ggplot() +
  geom_vline(xintercept = ymd("2023-07-26"),
             color = "red") +
  geom_col(aes(x = date,
               y = ntl_sum)) +
  geom_line(aes(x = date,
                y = ntl_sum_3m_ma,
                color = "3 Month\nMoving\nAverage"),
            linewidth = 1) +
  scale_color_manual(values = "darkorange") +
  theme_classic2() +
  labs(x = NULL,
       y = "Luminosity",
       title = "Trends in Nighttime Lights Across Niger",
       color = NULL)

Code

adm1_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm1_monthly.Rds"))

adm1_df %>%
  filter(date >= ymd("2019-01-01")) %>%
  ggplot() +
  geom_vline(xintercept = ymd("2023-07-26"),
             color = "red") +
  geom_col(aes(x = date,
               y = ntl_sum)) +
  geom_line(aes(x = date,
                y = ntl_sum_3m_ma,
                color = "3 Month\nMoving\nAverage"),
            linewidth = 1) +
  scale_color_manual(values = "darkorange") +
  theme_classic2() +
  theme_custom +
  labs(x = NULL,
       y = "Luminosity",
       color = NULL,
       title = "Trends in Nighttime Lights Across Regions") +
  facet_wrap(~ADM1_FR,
             scales = "free_y")

Code

city_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "city_monthly.Rds"))

city_df %>%
  filter(date >= ymd("2019-01-01")) %>%
  ggplot() +
  geom_vline(xintercept = ymd("2023-07-26"),
             color = "red") +
  geom_col(aes(x = date,
               y = ntl_sum)) +
  geom_line(aes(x = date,
                y = ntl_sum_3m_ma,
                color = "3 Month\nMoving\nAverage"),
            linewidth = 1) +
  scale_color_manual(values = "darkorange") +
  theme_classic2() +
  theme_custom +
  labs(x = NULL,
       color = NULL,
       y = "Luminosity",
       title = "Trends in Nighttime Lights Across Top Cities") +
  facet_wrap(~name,
             scales = "free_y")

Maps of change in nighttime lights

The below section shows the percent change in nighttime lights from August to December 2023 relative to August to December 2022. We show the percent change at the pixel level and at the ADM3 level.

Pixel level

Code

#### 2023
r_22_08 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                          paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2022_08",".tif")))
r_22_09 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                          paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2022_09",".tif")))
r_22_10 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                          paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2022_10",".tif")))
r_22_11 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                          paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2022_11",".tif")))
r_22_12 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                          paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2022_12",".tif")))

#### 2024
r_23_08 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                          paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_08",".tif")))
r_23_09 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                          paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_09",".tif")))
r_23_10 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                          paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_10",".tif")))
r_23_11 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                          paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_11",".tif")))
r_23_12 <- rast(file.path(ntl_dir, "individual_rasters", "monthly",
                          paste0("VNP46A3_NearNadir_Composite_Snow_Free_qflag_t","2023_12",".tif")))

#### Percent Change
r_22 <- r_22_08
r_23 <- r_23_08

r_22[] <- (r_22_08[] + r_22_09[] + r_22_10[] + r_22_11[] + r_22_12[]) / 5
r_23[] <- (r_23_08[] + r_23_09[] + r_23_10[] + r_23_11[] + r_23_12[]) / 5

r_23_pc <- r_23
r_23_pc[] <- (r_23[] - r_22[]) / r_22[] * 100

r_23_pc[][r_22[] <= 0.5] <- NA
r_23_pc[][r_23_pc[] >= 100] <- 100
r_23_pc[][r_23_pc[] <= -100] <- -100

r_23_pc <- r_23_pc %>% crop(adm0_sf) %>% mask(adm0_sf)

## Make map
r_values <- r_23_pc[] %>%
  unlist() %>%
  as.vector() %>%
  unique()

pal <- colorNumeric(c("red", "white", "green"), r_values,
                    na.color = "transparent")

leaflet() %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addRasterImage(r_23_pc, colors = pal, opacity = 1) %>%
  addLegend("bottomright", pal = pal, values = r_values,
            title = "% Change in NTL<br>Aug-Dec 2023<br>Relative to<br>Aug-Dec 2022",
            opacity = 1,
            labels = c("< -100", "-50", "0", "50", "> 100")
  )

ADM level

Static Map

Code

## Polygon
adm3_sf <- read_sf(file.path(boundaries_dir, 
                             "ner_adm_ignn_20230720_ab_shp",
                             "NER_admbnda_adm3_IGNN_20230720.shp"))

adm3_sf <- adm3_sf %>%
  select(-date)

## Data
adm3_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm3_monthly.Rds"))

adm3_sum_df <- adm3_df %>%
  dplyr::mutate(base_end = case_when(
    date %in% c(ymd("2022-08-01"),
                ymd("2022-09-01"),
                ymd("2022-10-01"),
                ymd("2022-11-01"),
                ymd("2022-12-01")) ~ "baseline",
    date %in% c(ymd("2023-08-01"),
                ymd("2023-09-01"),
                ymd("2023-10-01"),
                ymd("2023-11-01"),
                ymd("2023-12-01")) ~ "endline"
  )) %>%
  dplyr::filter(!is.na(base_end)) %>%
  group_by(ADM3_PCODE, base_end) %>%
  dplyr::summarise(ntl_sum = mean(ntl_sum)) %>%
  ungroup() %>%
  pivot_wider(id_cols = ADM3_PCODE,
              names_from = base_end,
              values_from = ntl_sum) %>%
  dplyr::mutate(pc = (endline - baseline) / baseline * 100)

adm3_sum_df$pc[adm3_sum_df$pc >= 100] <- 100

adm3_data_sf <- adm3_sf %>%
  left_join(adm3_sum_df, by = "ADM3_PCODE")

## Static Map
ggplot() +
  geom_sf(data = adm3_data_sf,
          aes(fill = pc)) +
  theme_void() +
  scale_fill_gradient2(low = "red",
                       high = "forestgreen") +
  labs(fill = "% Change",
       title = "Percent change in nighttime lights from Aug-Dec 2023 relative to Aug-Dec 2022")

Interactive Map

Click an administrative division to see the percent change

Code

## Interactive map
pal <- colorNumeric(c("red", "white", "green"), adm3_data_sf$pc,
                    na.color = "transparent")

adm3_data_sf <- adm3_data_sf %>%
  mutate(popup = paste0(ADM3_FR, "<br>% Change: ", round(pc, 2), "%"))

leaflet() %>%
  addProviderTiles(providers$CartoDB) %>%
  addPolygons(data = adm3_data_sf, 
              fillColor = ~pal(pc), 
              fillOpacity = 1,
              color = "black",
              weight = 1,
              popup = ~popup) %>%
  addLegend("bottomright", pal = pal, values = r_values,
            title = "% Change in NTL<br>Aug-Dec 2023<br>Relative to<br>Aug-Dec 2022",
            opacity = 1,
            labels = c("< -100", "-50", "0", "50", "> 100")
  )

Percent change in nighttime lights from 2022 to 2023

Below we show the percent change in nighttime lights for each month in 2023 relative to the same month in 2022

Code

## Load data
adm0_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm0_monthly.Rds"))
adm1_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm1_monthly.Rds"))
city_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "city_monthly.Rds"))

## Add dummpy ADM1 variable so function will more easily work on both datasets
adm0_df$ADM1_FR <- ""
city_df$ADM1_FR <- city_df$name

## Function to prep data as percent changes
prep_month_pc <- function(adm_df){
  
  adm_long_df <- adm_df %>%
    filter(date >= ymd("2022-01-01"),
           date <= ymd("2023-12-01")) %>%
    select(ADM1_FR, date, ntl_sum) %>%
    mutate(month = date %>% month,
           year = date %>% year) %>%
    pivot_wider(id_cols = c(ADM1_FR, month),
                values_from = ntl_sum,
                names_from = year) %>%
    #filter(!is.na(`2024`)) %>%
    mutate(pc = (`2023` - `2022`) / `2022` * 100)
  
  return(adm_long_df)
}

adm0_pc_df <- prep_month_pc(adm0_df)
adm1_pc_df <- prep_month_pc(adm1_df)
city_pc_df <- prep_month_pc(city_df)

## Maps
adm0_pc_df %>%
  ggplot() +
  geom_col(aes(x = month,
               y = pc)) +
  geom_vline(xintercept = 7, color = "red") +
  theme_classic2() +
  labs(x = "Month",
       y = "Percent change",
       title = "Percent change in nighttime lights across Niger from 2022 to 2023")

Code

adm1_pc_df %>%
  ggplot() +
  geom_col(aes(x = month,
               y = pc)) +
  geom_vline(xintercept = 7, color = "red") +
  theme_classic2() +
  theme_custom +
  labs(x = "Month",
       y = "Percent change",
       title = "Percent change in nighttime lights across regions from 2022 to 2023") +
  facet_wrap(~ADM1_FR,
             scales = "free_y")

Code

city_pc_df %>%
  ggplot() +
  geom_col(aes(x = month,
               y = pc)) +
  geom_vline(xintercept = 7, color = "red") +
  theme_classic2() +
  theme_custom +
  labs(x = "Month",
       y = "Percent change",
       title = "Percent change in nighttime lights across regions from 2023 to 2024") +
  facet_wrap(~ADM1_FR,
             scales = "free_y")

Percent change in nighttime lights from 2023 to 2024

Below we show the percent change in nighttime lights for each month in 2024 relative to the same month in 2023.

Code

## Load data
adm0_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm0_monthly.Rds"))
adm1_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "adm1_monthly.Rds"))
city_df <- readRDS(file.path(ntl_dir, "aggregated_appended", "city_monthly.Rds"))

## Add dummpy ADM1 variable so function will more easily work on both datasets
adm0_df$ADM1_FR <- ""
city_df$ADM1_FR <- city_df$name

## Function to prep data as percent changes
prep_month_pc <- function(adm_df){
  
  adm_long_df <- adm_df %>%
    filter(date >= ymd("2023-01-01"),
           date <= ymd("2024-06-01")) %>%
    select(ADM1_FR, date, ntl_sum) %>%
    mutate(month = date %>% month,
           year = date %>% year) %>%
    pivot_wider(id_cols = c(ADM1_FR, month),
                values_from = ntl_sum,
                names_from = year) %>%
    filter(!is.na(`2024`)) %>%
    mutate(pc = (`2024` - `2023`) / `2023` * 100)
  
  return(adm_long_df)
}

adm0_pc_df <- prep_month_pc(adm0_df)
adm1_pc_df <- prep_month_pc(adm1_df)
city_pc_df <- prep_month_pc(city_df)

## Maps
adm0_pc_df %>%
  ggplot() +
  geom_col(aes(x = month,
               y = pc)) +
  theme_classic2() +
  labs(x = "Month",
       y = "Percent change",
       title = "Percent change in nighttime lights across Niger from 2023 to 2024")

Code

adm1_pc_df %>%
  ggplot() +
  geom_col(aes(x = month,
               y = pc)) +
  theme_classic2() +
  theme_custom +
  labs(x = "Month",
       y = "Percent change",
       title = "Percent change in nighttime lights across regions from 2023 to 2024") +
  facet_wrap(~ADM1_FR,
             scales = "free_y")

Code

city_pc_df %>%
  ggplot() +
  geom_col(aes(x = month,
               y = pc)) +
  theme_classic2() +
  theme_custom +
  labs(x = "Month",
       y = "Percent change",
       title = "Percent change in nighttime lights across regions from 2023 to 2024") +
  facet_wrap(~ADM1_FR,
             scales = "free_y")