### Installation

Source: https://worldbank.github.io/blackmarbler/index.html

Instructions for installing the blackmarbler package from CRAN and the development version from GitHub.

```R
install.packages("blackmarbler")
```

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

--------------------------------

### Make raster

Source: https://worldbank.github.io/blackmarbler/index.html

Example of creating a raster from nighttime lights data for a specific date.

```R
r <- bm_raster(roi_sf = roi_sf,
               product_id = "VNP46A3",
               date = "2021-10-01",
               bearer = bearer)
```

--------------------------------

### Examples

Source: https://worldbank.github.io/blackmarbler/reference/download_h5_files.html

Example of how to download h5 files for Ghana for October 3, 2021, and create a raster using the downloaded files.

```R
if (FALSE) { # \dontrun{ 
# Define bearer token
bearer <- "BEARER-TOKEN-HERE"

# sf polygon of Ghana
library(geodata)
roi_sf <- gadm(country = "GHA", level=0, path = tempdir()) %>% st_as_sf()

# h5 files for Ghana for October 3, 2021
wget_h5_files(roi_sf = roi_sf,
              product_id = "VNP46A2",
              date = "2021-10-03",
              h5_dir = getwd(),        
              bearer = bearer)

# Make raster using h5_files
ken_202103_r <- bm_raster(roi_sf = roi_sf,
                          product_id = "VNP46A3",
                          date = "2021-03-01",
                          bearer = bearer,
                          h5_dir = getwd())

} # } 

```

--------------------------------

### Setup

Source: https://worldbank.github.io/blackmarbler/index.html

Loading necessary R packages, defining the NASA bearer token, and setting up a region of interest (ROI) using sf and geodata packages.

```R
#### 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()) 

```

--------------------------------

### Workflow to update data

Source: https://worldbank.github.io/blackmarbler/index.html

Example code for a workflow to regularly update daily nighttime lights data.

```R
# 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"))
```

--------------------------------

### Prep data

Source: https://worldbank.github.io/blackmarbler/index.html

Example of masking a raster to a region of interest and applying a log transformation.

```R
r <- r |> terra::mask(roi_sf)

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

--------------------------------

### Daily data in March 2021

Source: https://worldbank.github.io/blackmarbler/index.html

Example of how to retrieve daily nighttime lights data for March 2021.

```R
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)
```

--------------------------------

### Programmatically retrieve token

Source: https://worldbank.github.io/blackmarbler/index.html

Example of programmatically retrieving the NASA bearer token using username and password.

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

--------------------------------

### Extract Black Marble Data

Source: https://worldbank.github.io/blackmarbler/reference/bm_extract.html

Examples of extracting daily, monthly, and annual Black Marble raster data using the `bm_extract` function.

```R
if (FALSE) { # \dontrun{
# Define bearer token
bearer <- "BEARER-TOKEN-HERE"

# sf polygon of Ghana
library(geodata)
roi_sf <- gadm(country = "GHA", level=1, path = tempdir()) %>% st_as_sf()

# Daily data: raster for October 3, 2021
ken_20210205_r <- bm_extract(roi_sf = roi_sf,
                            product_id = "VNP46A2",
                            date = "2021-10-03",
                            bearer = bearer)

# Monthly data: raster for March 2021
ken_202103_r <- bm_extract(roi_sf = roi_sf,
                          product_id = "VNP46A3",
                          date = "2021-03-01",
                          bearer = bearer)

# Annual data: raster for 2021
ken_2021_r <- bm_extract(roi_sf = roi_sf,
                        product_id = "VNP46A4",
                        date = 2021,
                        bearer = bearer)
} # } 

```

--------------------------------

### Make raster of nighttime lights

Source: https://worldbank.github.io/blackmarbler/index.html

Examples of creating daily, monthly, and annual nighttime lights rasters for a specified region of interest.

```R
### 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)

```

--------------------------------

### Extract annual data

Source: https://worldbank.github.io/blackmarbler/index.html

Example of extracting annual nighttime lights data for Ghana's first administrative divisions.

```R
ntl_df <- bm_extract(roi_sf = roi_sf,
                     product_id = "VNP46A4",
                     date = 2012:2022,
                     bearer = bearer)
```

--------------------------------

### Map

Source: https://worldbank.github.io/blackmarbler/index.html

Example of creating a map of nighttime lights using ggplot2.

```R
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")
```

--------------------------------

### Monthly aggregated data in 2021 and 2022

Source: https://worldbank.github.io/blackmarbler/index.html

Example of how to retrieve monthly aggregated nighttime lights data for 2021 and 2022.

```R
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

Source: https://worldbank.github.io/blackmarbler/index.html

Example of how to retrieve yearly aggregated nighttime lights data for 2012 and 2021.

```R
r_annual <- bm_raster(roi_sf = roi_sf,
                      product_id = "VNP46A4",
                      date = 2012:2021,
                      bearer = bearer)
```

--------------------------------

### Trends over time

Source: https://worldbank.github.io/blackmarbler/index.html

Example of plotting trends in annual average nighttime lights for Ghana's administrative divisions.

```R
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"))
```

--------------------------------

### Retrieving Black Marble Raster Data

Source: https://worldbank.github.io/blackmarbler/reference/bm_raster.html

Examples demonstrating how to use the bm_raster function to fetch daily, monthly, and annual raster data for a specified region and date range. It includes setting up a bearer token and defining a spatial region of interest using sf polygons.

```r
if (FALSE) { # \dontrun{
# Define bearer token
bearer <- "BEARER-TOKEN-HERE"

# sf polygon of Ghana
library(geodata)
roi_sf <- gadm(country = "GHA", level=0, path = tempdir()) %>% st_as_sf()

# Daily data: raster for October 3, 2021
ken_20210205_r <- bm_raster(roi_sf = roi_sf,
                            product_id = "VNP46A2",
                            date = "2021-10-03",
                            bearer = bearer)

# Monthly data: raster for March 2021
ken_202103_r <- bm_raster(roi_sf = roi_sf,
                          product_id = "VNP46A3",
                          date = "2021-03-01",
                          bearer = bearer)

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

```

--------------------------------

### Usage

Source: https://worldbank.github.io/blackmarbler/reference/get_nasa_token.html

Fetch a NASA Earthdata bearer token from using the Earthdata API. If none exist, this will create one, or if one already exists it will fetch that one instead.

```R
get_nasa_token(username, password)
```

--------------------------------

### Download quality data for nighttime lights

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Downloads quality data for nighttime lights using the NearNadir_Composite_Snow_Free_Quality variable.

```r
quality_r <- bm_raster(roi_sf = roi_sf,
                       product_id = "VNP46A3",
                       date = "2023-01-01",
                       bearer = bearer,
                       variable = "NearNadir_Composite_Snow_Free_Quality")
```

--------------------------------

### Download monthly nighttime lights data

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Downloads nighttime lights data for January 2023 using the NearNadir_Composite_Snow_Free variable.

```r
ntl_r <- bm_raster(roi_sf = roi_sf,
                   product_id = "VNP46A3", 
                   date = "2023-01-01",
                   bearer = bearer,
                   variable = "NearNadir_Composite_Snow_Free")
```

--------------------------------

### Map of temporal gap for high quality pixels

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Prepares and visualizes the 'Latest_High_Quality_Retrieval' data as a map, showing the temporal gap in days.

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

##### Map 
ggplot() +
  geom_spatraster(data = ntl_tmp_gap_r) +
  scale_fill_distiller(palette = "Spectral",
                       na.value = "transparent") +
  coord_sf() + 
  theme_void() +
  labs(fill = "Temporal\nGap\n(Days)",
       title = "Temporal gap between date (Jan 1, 2023)\nand date of high quality pixel used") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))
```

--------------------------------

### Download 'Latest_High_Quality_Retrieval' data

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Downloads the 'Latest_High_Quality_Retrieval' variable, which indicates the temporal gap for gap filling.

```r
ntl_tmp_gap_r <- bm_raster(roi_sf = roi_sf,
                           product_id = "VNP46A2",
                           date = "2023-01-01",
                           bearer = bearer,
                           variable = "Latest_High_Quality_Retrieval")
```

--------------------------------

### Download number of observations for nighttime lights

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Downloads data to determine the number of observations used for nighttime light values, using the NearNadir_Composite_Snow_Free_Num variable.

```r
cf_r <- bm_raster(roi_sf = roi_sf,
                  product_id = "VNP46A3",
                  date = "2023-01-01",
                  bearer = bearer,
                  variable = "NearNadir_Composite_Snow_Free_Num")
```

--------------------------------

### Map of non-gap-filled nighttime lights

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Prepares and visualizes the non-gap-filled nighttime lights data as a map, applying a logarithmic scale.

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

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


```

--------------------------------

### Map of gap-filled nighttime lights

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Prepares and visualizes the gap-filled nighttime lights data as a map, applying a logarithmic scale and a specific color gradient.

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

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

##### Map 
ggplot() +
  geom_spatraster(data = ntl_m_r) +
  scale_fill_gradient2(low = "black",
                       mid = "yellow",
                       high = "red",
                       midpoint = 4,
                       na.value = "transparent") +
  coord_sf() + 
  theme_void() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        legend.position = "none")
```

--------------------------------

### Download nighttime lights, removing poor quality and gap-filled observations

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Downloads nighttime lights data, filtering out poor quality (1) and gap-filled (2) observations using the quality_flag_rm parameter.

```r
ntl_good_qual_r <- bm_raster(roi_sf = roi_sf,
                             product_id = "VNP46A3", 
                             date = "2023-01-01",
                             bearer = bearer,
                             variable = "NearNadir_Composite_Snow_Free",
                             quality_flag_rm = c(1,2)) # 1 = poor quality; 2 = gap filled based on historical data
```

--------------------------------

### Map monthly nighttime lights

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Prepares and maps the downloaded nighttime lights data, applying a logarithmic transformation and a diverging color scale.

```r
ntl_r <- ntl_r |> terra::mask(roi_sf) 

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

##### Map 
ggplot() +
  geom_spatraster(data = ntl_r) +
  scale_fill_gradient2(low = "black",
                       mid = "yellow",
                       high = "red",
                       midpoint = 4,
                       na.value = "transparent") +
  coord_sf() +
  theme_void() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        legend.position = "none")
```

--------------------------------

### Load packages and obtain a polygon for a region of interest (Switzerland)

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

This code block sets up the necessary R packages and defines a spatial object for Switzerland as the region of interest.

```r
library(blackmarbler)
library(geodata)
library(sf)
library(raster)
library(ggplot2)
library(tidyterra)
library(dplyr)
library(exactextractr)
library(lubridate)
library(tidyr)
library(geodata)
library(knitr)

bearer <- "BEARER-TOKEN-HERE"
```

--------------------------------

### Make raster of nighttime lights across multiple time periods

Source: https://worldbank.github.io/blackmarbler/index.html

Demonstrates how to extract nighttime lights data for multiple dates, months, or years, resulting in a SpatRaster object with multiple bands.

```R
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.

```

--------------------------------

### Download daily non-gap-filled nighttime lights data

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Downloads the 'DNB_BRDF-Corrected_NTL' variable, which represents nighttime light values directly from the selected date without gap filling.

```r
ntl_r <- bm_raster(roi_sf = roi_sf,
                   product_id = "VNP46A2",
                   date = "2023-01-01",
                   bearer = bearer,
                   variable = "DNB_BRDF-Corrected_NTL")
```

--------------------------------

### Map quality data

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Prepares and maps the quality data, assigning meaningful labels to quality levels and using a Spectral color palette.

```r
quality_r <- quality_r |> terra::mask(roi_sf) 

qual_levels <- data.frame(id=0:2, cover=c("0: Good quality", 
                                  "1: Poor quality", 
                                  "2: Gap filled"))

levels(quality_r) <- qual_levels

##### Map 
ggplot() +
  geom_spatraster(data = quality_r) +
  scale_fill_brewer(palette = "Spectral", 
                    direction = -1,
                    na.value = "transparent") + 
  labs(fill = "Quality") +
  coord_sf() +
  theme_void() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))
```

--------------------------------

### Map number of observations

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Prepares and maps the number of observations data using a viridis color scale.

```r
cf_r <- cf_r |> terra::mask(roi_sf) 

##### Map 
ggplot() +
  geom_spatraster(data = cf_r) +
  scale_fill_viridis_c(na.value = "transparent") + 
  labs(fill = "Number of\nObservations") +
  coord_sf() +
  theme_void() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))
```

--------------------------------

### Download daily gap-filled nighttime lights data

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Downloads the 'Gap_Filled_DNB_BRDF-Corrected_NTL' variable for January 1st, 2023, for the specified region of interest.

```r
ntl_r <- bm_raster(roi_sf = roi_sf,
                   product_id = "VNP46A2",
                   date = "2023-01-01",
                   bearer = bearer,
                   variable = "Gap_Filled_DNB_BRDF-Corrected_NTL")
```

--------------------------------

### Map nighttime lights for good quality observations

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

Prepares and maps the nighttime lights data filtered for good quality observations, applying a logarithmic transformation and a diverging color scale.

```r
ntl_good_qual_r <- ntl_good_qual_r |> terra::mask(roi_sf) 

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

##### Map 
ggplot() +
  geom_spatraster(data = ntl_good_qual_r) +
  scale_fill_gradient2(low = "black",
                       mid = "yellow",
                       high = "red",
                       midpoint = 4,
                       na.value = "transparent") +
  coord_sf() +
  theme_void() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        legend.position = "none")
```

--------------------------------

### Usage

Source: https://worldbank.github.io/blackmarbler/reference/download_h5_files.html

The usage of the download_h5_files function.

```R
download_h5_files(
  roi_sf = NULL,
  product_id,
  date,
  h5_dir,
  bearer,
  download_method = "httr"
)
```

--------------------------------

### Map Quality Flags

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

This code generates a map visualizing the quality flags for nighttime light data.

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

qual_levels <- data.frame(id=0:2, cover=c("0: High-quality, persistent", 
                                  "1: High-quality, ephemeral", 
                                  "2: Poor-quality"))

levels(quality_r) <- qual_levels

##### Map 
ggplot() +
  geom_spatraster(data = quality_r) +
  scale_fill_brewer(palette = "Spectral", 
                    direction = -1,
                    na.value = "transparent") + 
  labs(fill = "Quality") +
  coord_sf() +
  theme_void() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5))
```

--------------------------------

### Calculate Pixel Counts

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

This code defines functions to calculate the total number of pixels and non-NA pixels within a region, and then applies them to determine the proportion of non-NA nighttime light pixels.

```r
n_pixel <- function(values, coverage_fraction){
  length(values)
}

n_non_na_pixel <- function(values, coverage_fraction){
  sum(!is.na(values))
}

n_pixel_num        <- exact_extract(ntl_r, roi_sf, n_pixel)
n_non_na_pixel_num <- exact_extract(ntl_r, roi_sf, n_non_na_pixel)

print(n_pixel_num)
#> [1] 282934
print(n_non_na_pixel_num)
#> [1] 152285
print(n_non_na_pixel_num / n_pixel_num)
#> [1] 0.5382351
```

--------------------------------

### Map Nighttime Lights for Good Quality Observations

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

This code generates a map of nighttime lights, displaying only data from pixels with good quality flags.

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

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

##### Map 
ggplot() +
  geom_spatraster(data = ntl_good_qual_r) +
  scale_fill_gradient2(low = "black",
                       mid = "yellow",
                       high = "red",
                       midpoint = 4,
                       na.value = "transparent") +
  coord_sf() +
  theme_void() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        legend.position = "none")
```

--------------------------------

### License Text

Source: https://worldbank.github.io/blackmarbler/LICENSE-text.html

The license text for the blackmarbler project.

```text
YEAR: 2023
COPYRIGHT HOLDER: Robert Marty
```

--------------------------------

### Extract Nighttime Lights with Quality Filtering

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

This code extracts nighttime light data, filtering out pixels based on a specified quality flag removal value.

```r
ntl_good_qual_r <- bm_raster(roi_sf = roi_sf,
                             product_id = "VNP46A2", 
                             date = "2023-01-01",
                             bearer = bearer,
                             variable = "DNB_BRDF-Corrected_NTL",
                             quality_flag_rm = 2)
```

--------------------------------

### Extract Nighttime Light Data

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

This code uses the bm_extract function to retrieve nighttime light data for a specified date range and region.

```r
ntl_df <- bm_extract(roi_sf = roi_sf,
                     product_id = "VNP46A2",
                   date = seq.Date(from = ymd("2023-01-01"),
                                   to = ymd("2023-01-10"),
                                   by = 1),
                   bearer = bearer,
                   variable = "DNB_BRDF-Corrected_NTL")

knitr::kable(ntl_df)
```

--------------------------------

### Usage

Source: https://worldbank.github.io/blackmarbler/reference/bm_raster.html

Function signature for bm_raster.

```R
bm_raster(
  roi_sf,
  product_id,
  date,
  bearer,
  variable = NULL,
  quality_flag_rm = NULL,
  check_all_tiles_exist = TRUE,
  interpol_na = FALSE,
  output_location_type = "memory",
  file_dir = NULL,
  file_prefix = NULL,
  file_skip_if_exists = TRUE,
  file_return_null = FALSE,
  h5_dir = NULL,
  download_method = "httr",
  quiet = FALSE,
  ...
)
```

--------------------------------

### Plot Trends in Nighttime Lights and Data Coverage

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

This code generates a plot showing trends in average nighttime lights and the proportion of the country with nighttime light data over time.

```r
ntl_df %>%
  dplyr::select(date, ntl_sum, prop_non_na_pixels) %>%
  pivot_longer(cols = -date) %>%
  
  ggplot(aes(x = date,
             y = value)) + 
  geom_line() + 
  facet_wrap(~name,
             scales = "free")
```

--------------------------------

### bm_extract Function Usage

Source: https://worldbank.github.io/blackmarbler/reference/bm_extract.html

This snippet shows the usage of the bm_extract function, detailing its parameters for extracting and aggregating Black Marble data.

```R
bm_extract(
  roi_sf,
  product_id,
  date,
  bearer,
  aggregation_fun = c("sum"),
  add_n_pixels = TRUE,
  variable = NULL,
  quality_flag_rm = NULL,
  check_all_tiles_exist = TRUE,
  interpol_na = FALSE,
  output_location_type = "memory",
  file_dir = NULL,
  file_prefix = NULL,
  file_skip_if_exists = TRUE,
  file_return_null = FALSE,
  h5_dir = NULL,
  download_method = "httr",
  quiet = FALSE,
  ...
)
```

--------------------------------

### Extract Mandatory Quality Flag Raster

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

This code extracts the Mandatory Quality Flag raster data for a specific date and region.

```r
quality_r <- bm_raster(roi_sf = roi_sf,
                       product_id = "VNP46A2",
                       date = "2023-01-01",
                       bearer = bearer,
                       variable = "Mandatory_Quality_Flag")
```

--------------------------------

### Define region of interest (ROI)

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

This code snippet defines the region of interest using the gadm function for Switzerland.

```r
roi_sf <- gadm(country = "CHE", level=0, path = tempdir()) |> st_as_sf()
```

--------------------------------

### Map Nighttime Lights

Source: https://worldbank.github.io/blackmarbler/articles/assess-quality.html

This code snippet generates a map of nighttime lights using ggplot2 and geom_spatraster.

```r
ggplot() +
  geom_spatraster(data = ntl_m_r) +
  scale_fill_gradient2(low = "black",
                       mid = "yellow",
                       high = "red",
                       midpoint = 4,
                       na.value = "transparent") +
  coord_sf() +
  theme_void() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5),
        legend.position = "none")
```

--------------------------------

### Citation in DESCRIPTION format

Source: https://worldbank.github.io/blackmarbler/authors.html

The citation for the blackmarbler package as it might appear in a DESCRIPTION file.

```r
@Manual{
  title = {blackmarbler: Black Marble Data and Statistics},
  author = {Robert Marty and Gabriel {Stefanini Vicente}},
  year = {2025},
  note = {R package version 0.2.6},
  url = {https://worldbank.github.io/blackmarbler/},
}
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.