BlackMarbleR

CRAN_Status_Badge downloads GitHub Repo stars activity License: MIT R-CMD-check

BlackMarbleR is a R package that provides a simple way to use nighttime lights data from NASA’s Black Marble. Black Marble is a NASA Earth Science Data Systems (ESDS) project that provides a product suite of daily, monthly and yearly global nighttime lights. This package automates the process of downloading all relevant tiles from the NASA LAADS DAAC to cover a region of interest, converting and mosaicing the raw files (in HDF5 format) to georeferenced rasters.

Installation

The package can be installed via CRAN.

install.packages("blackmarbler")

To install the development version from Github:

# install.packages("devtools")
devtools::install_github("worldbank/blackmarbler")

Bearer Token

To obtain a bearer token, you’ll need to have a registered NASA Earth Data account. On the webpage, click “login” and create a username/password if needed.

After an account is created, the NASA Bearer Token can be retrieved either using the get_nasa_token function or manually (see below).

Programmatically retrieve token

The NASA Bearer Token can also be programmatically retrieved using the get_nasa_token() function. After making an account, the get_nasa_token() function uses your username and password to retrieve the Bearer token.

bearer <- get_nasa_token(username = "USERNAME-HERE", 
                         password = "PASSWORD-HERE")

Manually retrieve token

The function requires using a Earthdata Download Bearer Token; to obtain a token, follow the below steps:

  1. Go to the NASA LAADS Archive
  2. Click “Login” (button on top right), and click “Earthdata Login”; create an account if needed.
  3. Enter your username and password.
  4. Click “Login”, then “Generate Token” on the dropdown.
  5. Click the “Generate Token” tab, then click the green button to generate a token
  6. Click the blue “Show token” button; this is your bearer token. It will be a long string of text (over 500 characters).

NASA LAADS Bearer Token

NASA LAADS Bearer Token

NASA LAADS Bearer Token

Usage

Setup

Before downloading and extracting Black Marble data, we first load packages, define the NASA bearer token, and define a region of interest.

#### Setup
# Load packages
library(blackmarbler)
library(geodata)
library(sf)
library(terra)
library(ggplot2)
library(tidyterra)
library(lubridate)

#### Define NASA bearer token
bearer <- "BEARER-TOKEN-HERE"

### ROI
# Define region of interest (roi). The roi must be (1) an sf polygon and (2)
# in the WGS84 (epsg:4326) coordinate reference system. Here, we use the
# getData function to load a polygon of Ghana
roi_sf <- gadm(country = "GHA", level=1, path = tempdir())

Make raster of nighttime lights

The below example shows making daily, monthly, and annual rasters of nighttime lights for Ghana.

### Daily data: raster for February 5, 2021
r_20210205 <- bm_raster(roi_sf = roi_sf,
                        product_id = "VNP46A2",
                        date = "2021-02-05",
                        bearer = bearer)

### Monthly data: raster for October 2021
r_202110 <- bm_raster(roi_sf = roi_sf,
                      product_id = "VNP46A3",
                      date = "2021-10-01", # The day is ignored
                      bearer = bearer)

### Annual data: raster for 2021
r_2021 <- bm_raster(roi_sf = roi_sf,
                    product_id = "VNP46A4",
                    date = 2021,
                    bearer = bearer)

Make raster of nighttime lights across multiple time periods

To extract data for multiple time periods, add multiple time periods to date. The function will return a SpatRaster object with multiple bands, where each band corresponds to a different date. The below code provides examples getting data across multiple days, months, and years.

#### Daily data in March 2021
r_daily <- bm_raster(roi_sf = roi_sf,
                     product_id = "VNP46A2",
                     date = seq.Date(from = ymd("2021-03-01"), to = ymd("2021-03-31"), by = "day"),
                     bearer = bearer)

#### Monthly aggregated data in 2021 and 2022
r_monthly <- bm_raster(roi_sf = roi_sf,
                       product_id = "VNP46A3",
                       date = seq.Date(from = ymd("2021-01-01"), to = ymd("2022-12-01"), by = "month"),
                       bearer = bearer)

#### Yearly aggregated data in 2012 and 2021
r_annual <- bm_raster(roi_sf = roi_sf,
                      product_id = "VNP46A4",
                      date = 2012:2021,
                      bearer = bearer)

Map of nighttime lights

Using one of the rasters, we can make a map of nighttime lights

#### Make raster
r <- bm_raster(roi_sf = roi_sf,
               product_id = "VNP46A3",
               date = "2021-10-01",
               bearer = bearer)

#### Prep data
r <- r |> terra::mask(roi_sf)

## Distribution is skewed, so log
r[] <- log(r[] + 1)

##### Map
ggplot() +
  geom_spatraster(data = r) +
  scale_fill_gradient2(low = "black",
                       mid = "yellow",
                       high = "red",
                       midpoint = 4.5,
                       na.value = "transparent") +
  labs(title = "Nighttime Lights: October 2021") +
  coord_sf() +
  theme_void() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
  legend.position = "none")

Nighttime Lights Map

We can use the bm_extract function to observe changes in nighttime lights over time. The bm_extract function leverages the exactextractr package to aggregate nighttime lights data to polygons. Below we show trends in annual nighttime lights data across Ghana’s first administrative divisions.

#### Extract annual data
ntl_df <- bm_extract(roi_sf = roi_sf,
                     product_id = "VNP46A4",
                     date = 2012:2022,
                     bearer = bearer)

#### Trends over time
ntl_df |>
  ggplot() +
  geom_col(aes(x = date,
  y = ntl_mean),
  fill = "darkorange") +
  facet_wrap(~NAME_1) +
  labs(x = NULL,
       y = "NTL Luminosity",
       title = "Ghana Admin Level 1: Annual Average Nighttime Lights") +
  scale_x_continuous(labels = seq(2012, 2022, 4),
                     breaks = seq(2012, 2022, 4)) +
  theme_minimal() +
  theme(strip.text = element_text(face = "bold"))

Nighttime Lights Trends

Workflow to update data

Some users may want to monitor near-real-time changes in nighttime lights. For example, daily Black Marble nighttime lights data is updated regularly, where data is available roughly on a week delay; same use cases may require examining trends in daily nighttime lights data as new data becomes available. Below shows example code that could be regularly run to produce an updated daily dataset of nighttime lights.

The below code produces a dataframe of nighttime lights for each date, where average nighttime lights for Ghana’s 1st administrative division is produced. The code will check whether data has already been downloaded/extracted for a specific date, and only download/extract new data.

# Create directories to store data
dir.create(file.path(getwd(), "bm_files"))
dir.create(file.path(getwd(), "bm_files", "daily"))

# Extract daily-level nighttime lights data for Ghana's first administrative divisions.
# Save a separate dataset for each date in the `"~/Desktop/bm_files/daily"` directory.
# The code extracts data from January 1, 2023 to today. Given that daily nighttime lights
# data is produced on roughly a week delay, the function will only extract data that exists;
# it will skip extracting data for dates where data has not yet been produced by NASA Black Marble.
bm_extract(roi_sf = roi_sf,
           product_id = "VNP46A2",
           date = seq.Date(from = ymd("2023-01-01"), to = Sys.Date(), by = 1),
           bearer = bearer,
           output_location_type = "file",
           file_dir = file.path(getwd(), "bm_files", "daily"))

# Append daily-level datasets into one file
file.path(getwd(), "bm_files", "daily") |>
  list.files(pattern = "*.Rds",
  full.names = T) |>
  map_df(readRDS) |>
  saveRDS(file.path(getwd(), "bm_files", "ntl_daily.Rds"))

Functions and arguments

Functions

The package provides two functions.

Both functions take the following arguments:

Required arguments

Optional arguments

If output_location_type = "file", the following arguments can be used:

Argument for bm_extract only

Black Marble Resources

For more information on NASA Black Marble, see: