Using BlackMarblePy

Using BlackMarblePy#

Open In Colab Binder

This Jupyter notebook provides a guided exploration of the BlackMarblePy package, showcasing its capabilities for downloading, visualizing, and analyzing NASA Black Marble nighttime lights data. Through interactive examples, you’ll learn how to:

  • Download daily, monthly, and yearly data for specific dates, regions, or bounding boxes.

  • Visualize downloaded data in various forms, including maps and bar charts.

  • Save visualizations and analysis results for further use.

Requirements#

Before downloading or extracting NASA Black Marble data, make sure the following requirements are met:

📦 Install blackmarblepy#

To use this library, install the blackmarblepy package. It is highly recommended to install it inside a Python virtual environment to avoid dependency conflicts.

BlackMarblePy is available on PyPI as blackmarblepy and can installed using pip:

pip install blackmarblepy
Hide code cell content 
import os
import datetime
import logging

import colorcet as cc
import contextily as cx
import geopandas
import matplotlib.pyplot as plt
import pandas as pd
from bokeh.plotting import figure, output_notebook, show
from bokeh.models import HoverTool, Title

from blackmarble import BlackMarble

from dotenv import load_dotenv

logging.getLogger("blackmarblepy").setLevel(logging.INFO)
load_dotenv()

%load_ext autoreload
%load_ext watermark
%autoreload 2
%watermark -v -u -n -p blackmarble

Last updated: Tue Jun 10 2025

Python implementation: CPython
Python version       : 3.12.9
IPython version      : 9.3.0

blackmarble: 2025.6.3.dev3+g830d70c.d20250611

âś… Set Up NASA Earthdata Token#

BlackMarblePy requires a valid, unexpired NASA Earthdata bearer token, which you can retrieve from your Earthdata profile. For ease and security, we recommend setting this as the BLACKMARBLE_TOKEN environment variable on your system.

  1. Access Earthdata Login. In case you haven’t already, you must register.

    ../_images/nasa_earthdata_profile.png

  2. Select Generate Token. If the token is expired or you are in need of one, click the Generate Token button.

    ../_images/nasa_earthdata_generate_token.png

    Caution

    Please be aware that the “Affiliation” information on your Earthdata profile is mandatory. Without this information, the NASA Earthdata token will be invalid. In case of a download error, try visiting the URL which is failing, as you may be prompted to grant permissions.

  3. Use your token securely. We recommend setting it as a secret or an environment variable rather than hardcoding it. For example, you can set the BLACKMARBLE_TOKEN environment variable to store your token safely in Unix-like systems:

    export BLACKMARBLE_TOKEN=<your_nasa_earthdata-token>

    Important

    Using a secret token securely in Python code involves several practices to ensure the token is not exposed unintentionally. For instance, storing the secret token in an environment variable, in configuration files or using secret management services. In this example, we set up an environment variable.

    See also

🌍 Define Region of Interest#

You must define a region of interest as a geopandas.GeoDataFrame. This should represent the geographic area you want to extract data for.

For example, we obtain the polygon below from GADM for Ghana.

Hide code cell source 
gdf = geopandas.read_file(
    "https://geodata.ucdavis.edu/gadm/gadm4.1/json/gadm41_GHA_1.json.zip"
)
gdf.explore(tiles="CartoDB dark_matter")
Make this Notebook Trusted to load map: File -> Trust Notebook
../images/favicon.ico

Fig. 1 This map of Ghana displays the administrative boundaries as obtained from the Global Administrative Areas (GADM) database. The map highlights the regional divisions and key cities, providing a detailed geographic overview of the country. Data source: GADM, version 4.0.#

Examples#

In this section, we will demonstrate how to use BlackMarblePy to download and manipulate NASA Black Marble data.

BlackMarblePy offers two main interfaces:

  • A class-based interface (BlackMarble) that retains configuration details like the bearer token and output directory across multiple calls.

  • A procedural API (bm_extract, bm_raster) that provides the same functionality without requiring an instantiated object, making it convenient for one-off operations.

Both approaches give you flexible options depending on your workflow. We instantiate BlackMarble to use on the examples below.

# Class-based interface example

# Initialize the BlackMarble interface.
bm = BlackMarble()

# Optional: explicitly configure all options
# bm = BlackMarble(
#     bearer="your_token_here",                     # NASA Earthdata bearer token (can also be set via env var)
#     check_all_tiles_exist=True,                   # Skip dates if any tile is missing
#     drop_values_by_quality_flag=[255],            # Mask out invalid data (e.g., fill value)
#     output_directory="data",                      # Directory to save outputs
#     output_skip_if_exists=True                    # Skip downloading if file already exists
# )

Create raster of nighttime lights#

In this section, we show examples of creating daily, monthly, and annual rasters of nighttime lights for the Region of Interest selected.

Daily#

# Daily data: raster for February 5, 2021
VNP46A2_20210205 = bm.raster(
    gdf,
    product_id="VNP46A2",
    date_range="2021-02-05",
)
VNP46A2_20210205
Hide code cell output 
<xarray.Dataset> Size: 13MB
Dimensions:                            (x: 1071, y: 1544, time: 1)
Coordinates:
  * x                                  (x) float64 9kB -3.26 -3.256 ... 1.198
  * y                                  (y) float64 12kB 11.17 11.17 ... 4.744
  * time                               (time) datetime64[ns] 8B 2021-02-05
Data variables:
    Gap_Filled_DNB_BRDF-Corrected_NTL  (time, y, x) float64 13MB nan nan ... nan
Attributes: (12/42)
    AlgorithmType:                     b'SCI'
    DataResolution:                    b'Moderate'
    DayNightFlag:                      b'Day'
    EastBoundingCoord:                 0.0
    EndTime:                           b'2021-02-05 23:59:59.000'
    GranuleDayNightFlag:               b'Day'
    ...                                ...
    VerticalTileNumber:                b'07'
    WestBoundingCoord:                 -10.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0
    _FillValue:                        nan
Hide code cell source 
fig, ax = plt.subplots(figsize=(8, 8))

VNP46A2_20210205["Gap_Filled_DNB_BRDF-Corrected_NTL"].sel(
    time="2021-02-05"
).plot.pcolormesh(
    ax=ax,
    cmap=cc.cm.bmy,
    robust=True,
)
cx.add_basemap(ax, crs=gdf.crs.to_string())

ax.text(
    0,
    -0.1,
    "Source: NASA Black Marble VNP46A2",
    ha="left",
    va="center",
    transform=ax.transAxes,
    fontsize=10,
    color="black",
    weight="normal",
)
ax.set_title("Ghana: NTL Radiance on Feb 5 2021", fontsize=20);
../_images/b456bf338b5564c625df6c7abb6377822fe66392c82954da4f243679283b62be.png

Monthly#

# Monthly data: raster for October 2021
VNP46A3_202110 = bm.raster(
    gdf,
    product_id="VNP46A3",
    date_range="2021-10-01",
)
VNP46A3_202110
Hide code cell output 
<xarray.Dataset> Size: 13MB
Dimensions:                        (x: 1071, y: 1544, time: 1)
Coordinates:
  * x                              (x) float64 9kB -3.26 -3.256 ... 1.194 1.198
  * y                              (y) float64 12kB 11.17 11.17 ... 4.748 4.744
  * time                           (time) datetime64[ns] 8B 2021-10-01
Data variables:
    NearNadir_Composite_Snow_Free  (time, y, x) float64 13MB nan nan ... nan nan
Attributes: (12/52)
    AlgorithmType:                     b'SCI'
    AlgorithmVersion:                  b'NPP_PR46A3 2.0.0'
    Conventions:                       b'CF-1.6'
    creator_email:                     b'modis-ops@lists.nasa.gov'
    creator_name:                      b'VIIRS Land SIPS Processing Group'
    creator_url:                       b'https://ladsweb.modaps.eosdis.nasa.gov'
    ...                                ...
    VerticalTileNumber:                b'08'
    WestBoundingCoord:                 0.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0
    _FillValue:                        nan
Hide code cell source 
fig, ax = plt.subplots(figsize=(8, 8))

VNP46A3_202110["NearNadir_Composite_Snow_Free"].sel(time="2021-10-01").plot.pcolormesh(
    ax=ax,
    cmap=cc.cm.bmy,
    robust=True,
)
cx.add_basemap(ax, crs=gdf.crs.to_string())

ax.text(
    0,
    -0.1,
    "Source: NASA Black Marble VNP46A3",
    ha="left",
    va="center",
    transform=ax.transAxes,
    fontsize=10,
    color="black",
    weight="normal",
)
ax.set_title("Ghana: NTL Radiance in Oct 2021", fontsize=20);
../_images/150446735a255f3a944e97e56504ee07bd22d0fd92060d4e9a415204397618de.png

Annual#

# Annual data: raster for 2021
VNP46A4_2021 = bm.raster(
    gdf,
    product_id="VNP46A4",
    date_range="2021-01-01",
)
VNP46A4_2021
Hide code cell output 
<xarray.Dataset> Size: 13MB
Dimensions:                        (x: 1071, y: 1544, time: 1)
Coordinates:
  * x                              (x) float64 9kB -3.26 -3.256 ... 1.194 1.198
  * y                              (y) float64 12kB 11.17 11.17 ... 4.748 4.744
  * time                           (time) datetime64[ns] 8B 2021-01-01
Data variables:
    NearNadir_Composite_Snow_Free  (time, y, x) float64 13MB nan nan ... nan nan
Attributes: (12/52)
    AlgorithmType:                     b'SCI'
    AlgorithmVersion:                  b'NPP_PR46A3 2.0.0'
    Conventions:                       b'CF-1.6'
    creator_email:                     b'modis-ops@lists.nasa.gov'
    creator_name:                      b'VIIRS Land SIPS Processing Group'
    creator_url:                       b'https://ladsweb.modaps.eosdis.nasa.gov'
    ...                                ...
    VerticalTileNumber:                b'07'
    WestBoundingCoord:                 0.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0
    _FillValue:                        nan
Hide code cell source 
fig, ax = plt.subplots(figsize=(8, 8))

VNP46A4_2021["NearNadir_Composite_Snow_Free"].sel(time="2021-01-01").plot.pcolormesh(
    ax=ax,
    cmap=cc.cm.bmy,
    robust=True,
)
cx.add_basemap(ax, crs=gdf.crs.to_string())

ax.text(
    0,
    -0.1,
    "Source: NASA Black Marble VNP46A4",
    ha="left",
    va="center",
    transform=ax.transAxes,
    fontsize=10,
    color="black",
    weight="normal",
)
ax.set_title("Ghana: NTL Radiance in 2021", fontsize=20);
../_images/7068d60d9be5a9bc078d47184eee228c126c93b5a90a04de41c5295591489baa.png

Create a raster stack of nighttime lights across multiple time periods#

In this section, we illustrate how to retrieve and extract NASA Black Marble data for multiple time periods. The function will return a raster stack, where each raster band corresponds to a different date. The following code snippet provides examples of getting data across multiple days, months, and years. For each example, we define a date range using pd.date_range.

# Raster stack of daily data
r_daily = bm.raster(
    gdf,
    product_id="VNP46A2",
    date_range=pd.date_range("2022-01-01", "2022-03-31", freq="D"),
)

<xarray.Dataset> Size: 1GB
Dimensions:                            (x: 1071, y: 1544, time: 90)
Coordinates:
  * x                                  (x) float64 9kB -3.26 -3.256 ... 1.198
  * y                                  (y) float64 12kB 11.17 11.17 ... 4.744
  * time                               (time) datetime64[ns] 720B 2022-01-01 ...
Data variables:
    Gap_Filled_DNB_BRDF-Corrected_NTL  (time, y, x) float64 1GB nan nan ... nan
Attributes: (12/32)
    AlgorithmType:                     b'SCI'
    DataResolution:                    b'Moderate'
    DayNightFlag:                      b'Day'
    EastBoundingCoord:                 0.0
    GranuleDayNightFlag:               b'Day'
    GRingPointLatitude:                [ 0. 10. 10.  0.]
    ...                                ...
    VerticalTileNumber:                b'08'
    WestBoundingCoord:                 -10.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0
    _FillValue:                        nan
Hide code cell source 
# Create the figure and axis
fig, ax = plt.subplots(figsize=(16, 8))

# Plot the mean NTL radiance over the dimensions x and y
r_daily["Gap_Filled_DNB_BRDF-Corrected_NTL"].mean(dim=["x", "y"]).plot(ax=ax)

# Add the data source text
ax.text(
    0,
    -0.2,
    "Source: NASA Black Marble VNP46A2",
    ha="left",
    va="center",
    transform=ax.transAxes,
    fontsize=10,
    color="black",
    weight="normal",
)

# Set the title with appropriate fontsize
ax.set_title("Ghana: Daily Average NTL Radiance", fontsize=20, weight="bold")

# Add labels to the axes
ax.set_xlabel("Date", fontsize=12)
ax.set_ylabel("Radiance (nW/cm²/sr)", fontsize=12)

# Adjust layout for better spacing
fig.tight_layout()

# Show the plot
plt.show()
../_images/808c750ca678629478fba7f33cbecc0c24e7fb9dac6fbb64d90e9c76d0eb40dc.png
../images/favicon.ico

Fig. 2 This figures describes the daily average nighttime lights radiance data plotted over time. The data reflects fluctuations in radiance levels due to varying cloud cover, affecting the accuracy of the measurements#

# Raster stack of monthly data
r_monthly = bm.raster(
    gdf,
    product_id="VNP46A3",
    date_range=pd.date_range("2022-01-01", "2022-12-31", freq="MS"),
)

<xarray.Dataset> Size: 159MB
Dimensions:                        (x: 1071, y: 1544, time: 12)
Coordinates:
  * x                              (x) float64 9kB -3.26 -3.256 ... 1.194 1.198
  * y                              (y) float64 12kB 11.17 11.17 ... 4.748 4.744
  * time                           (time) datetime64[ns] 96B 2022-01-01 ... 2...
Data variables:
    NearNadir_Composite_Snow_Free  (time, y, x) float64 159MB nan nan ... nan
Attributes: (12/40)
    AlgorithmType:                     b'SCI'
    AlgorithmVersion:                  b'NPP_PR46A3 2.0.0'
    Conventions:                       b'CF-1.6'
    creator_email:                     b'modis-ops@lists.nasa.gov'
    creator_name:                      b'VIIRS Land SIPS Processing Group'
    creator_url:                       b'https://ladsweb.modaps.eosdis.nasa.gov'
    ...                                ...
    VerticalTileNumber:                b'07'
    WestBoundingCoord:                 0.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0
    _FillValue:                        nan
Hide code cell source 
# Create the figure and axis
fig, ax = plt.subplots(figsize=(16, 8))

# Plot the mean NTL radiance over the dimensions x and y
r_monthly["NearNadir_Composite_Snow_Free"].mean(dim=["x", "y"]).plot(ax=ax)

# Add the data source text
ax.text(
    0,
    -0.2,
    "Source: NASA Black Marble VNP46A3",
    ha="left",
    va="center",
    transform=ax.transAxes,
    fontsize=10,
    color="black",
    weight="normal",
)

# Set the title with appropriate fontsize
ax.set_title("Ghana: Monthly Average NTL Radiance", fontsize=20, weight="bold")

# Add labels to the axes
ax.set_xlabel("Date", fontsize=12)
ax.set_ylabel("Radiance (nW/cm²/sr)", fontsize=12)

# Adjust layout for better spacing
fig.tight_layout()

# Show the plot
plt.show()
../_images/6dc67dab984e5a6c264a85404caf949aab89dfa6da0953468e78b494d48eb2c5.png
../images/favicon.ico

Fig. 3 This figures describes the monthly average nighttime lights radiance data plotted over time. The data reflects fluctuations in radiance levels due to varying cloud cover, affecting the accuracy of the measurements#

# Raster stack of yearly data
r_yearly = bm.raster(
    gdf,
    product_id="VNP46A4",
    date_range=pd.date_range("2019-01-01", "2022-01-01", freq="YS"),
)

<xarray.Dataset> Size: 53MB
Dimensions:                        (x: 1071, y: 1544, time: 4)
Coordinates:
  * x                              (x) float64 9kB -3.26 -3.256 ... 1.194 1.198
  * y                              (y) float64 12kB 11.17 11.17 ... 4.748 4.744
  * time                           (time) datetime64[ns] 32B 2019-01-01 ... 2...
Data variables:
    NearNadir_Composite_Snow_Free  (time, y, x) float64 53MB nan nan ... nan nan
Attributes: (12/40)
    AlgorithmType:                     b'SCI'
    AlgorithmVersion:                  b'NPP_PR46A3 2.0.0'
    Conventions:                       b'CF-1.6'
    DataResolution:                    b'15 arc second'
    DayNightFlag:                      b'Night'
    HorizontalTileNumber:              b'17'
    ...                                ...
    NorthBoundingCoord:                20.0
    publisher_email:                   b'modis-ops@lists.nasa.gov'
    publisher_name:                    b'LAADS'
    publisher_url:                     b'https://ladsweb.modaps.eosdis.nasa.gov'
    SouthBoundingCoord:                10.0
    WestBoundingCoord:                 -10.0
Hide code cell source 
# Set up the figure and subplots
fig, axs = plt.subplots(2, 2, figsize=(16, 16))

for i, t in enumerate(r_yearly["time"]):
    row = i // 2
    col = i % 2
    ax = axs[row, col]

    r_yearly["NearNadir_Composite_Snow_Free"].sel(time=t).plot.pcolormesh(
        ax=ax,
        cmap=cc.cm.bmw,
        robust=True,
        vmax=50,
    )
    cx.add_basemap(ax, crs=gdf.crs.to_string())

plt.show()
../_images/92de36d710411337a651e53e0e003c4aac1aa3779021b2b7e4ce168e648f1d7b.png
../images/favicon.ico

Fig. 4 Temporal variation of NearNadir_Composite_Snow_Free mapped over multiple years. Each subplot represents a different time period overlaid with a basemap.#

Visualizing difference in radiance year over year#

Lastly, we calculate the increase/decrease in nighttime lights radiance levels.

Hide code cell source 
fig, ax = plt.subplots(figsize=(16, 16))

delta = (
    (
        (
            r_yearly["NearNadir_Composite_Snow_Free"].sel(time="2022-01-01")
            - r_yearly["NearNadir_Composite_Snow_Free"].sel(time="2019-01-01")
        )
        / r_yearly["NearNadir_Composite_Snow_Free"].sel(time="2019-01-01")
    )
    # .drop("time")
    .plot.pcolormesh(ax=ax, cmap="Spectral", robust=True)
)
cx.add_basemap(ax, crs=gdf.crs.to_string(), source=cx.providers.CartoDB.DarkMatter)

ax.text(
    0,
    -0.1,
    "Source: NASA Black Marble VNP46A3",
    ha="left",
    va="center",
    transform=ax.transAxes,
    fontsize=10,
    color="black",
    weight="normal",
)
ax.set_title("Ghana: NTL Radiance Increase/Decrease (2019-2022)", fontsize=16);
../_images/d1c9e5248916290bd2618b236b1e0530b6fa48181969f7bc582094d827de514d.png
../images/favicon.ico

Fig. 5 This figure displays the percentage change in radiance for Ghana between 2019 and 2022. The data, sourced from NASA’s Black Marble VNP46A3 product is visualized with basemap.#

Downloading and Storing Black Marble data#

In this section, we provide a guide on using BlackMarblePy to download NASA Black Marble data and save it to a specified local directory. You can use the output_directory parameter to designate the directory for saving the files. By default, files that have already been downloaded will not be re-downloaded in subsequent executions.

BlackMarble(
    output_directory="data/",  # Save files to local directory
    output_skip_if_exists=True,  # Set to skip if exists
).raster(
    gdf,
    product_id="VNP46A4",
    date_range=pd.date_range("2022-01-01", "2022-01-01", freq="YS"),
)

Alternatively, set output_skip_if_exists=False to force the redownload of files from NASA.

References#

1

Miguel O. Román, Zhuosen Wang, Qingsong Sun, Virginia Kalb, Steven D. Miller, Andrew Molthan, Lori Schultz, Jordan Bell, Eleanor C. Stokes, Bhartendu Pandey, Karen C. Seto, Dorothy Hall, Tomohiro Oda, Robert E. Wolfe, Gary Lin, Navid Golpayegani, Sadashiva Devadiga, Carol Davidson, Sudipta Sarkar, Cid Praderas, Jeffrey Schmaltz, Ryan Boller, Joshua Stevens, Olga M. Ramos González, Elizabeth Padilla, José Alonso, Yasmín Detrés, Roy Armstrong, Ismael Miranda, Yasmín Conte, Nitza Marrero, Kytt MacManus, Thomas Esch, and Edward J. Masuoka. Nasa's black marble nighttime lights product suite. Remote Sensing of Environment, 210:113–143, 2018. URL: https://www.sciencedirect.com/science/article/pii/S003442571830110X, doi:https://doi.org/10.1016/j.rse.2018.03.017.

2

Gabriel Stefanini Vicente and Robert Marty. BlackMarblePy: Georeferenced Rasters and Statistics of Nighttime Lights from NASA Black Marble. https://worldbank.github.io/blackmarblepy, 2023. URL: https://worldbank.github.io/blackmarblepy, doi:10.5281/zenodo.10667907.

Contents