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#

In this example, we use additional dependencies included below below for convenience. However, it is strongly recommended to install requirements on a virtual environment.

Hide code cell source 
import os
import datetime

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.extract import bm_extract
from blackmarble.raster import bm_raster

%load_ext autoreload
%autoreload 2

plt.rcParams["figure.figsize"] = (18, 10)

Generate NASA Earthdata Token#

BlackMarblePy requires a NASA Earthdata bearer token. To generate a token, you may use Earthdata Login

  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, which may result in an error (i.e., OSError: Unable to synchronously open file (file signature not found))

  3. Copy and use the token with BlackMarblePy securely. For example, by setting secret or an environment variable.

     # An environment variable can obfuscate to secure a secret
 import os
 
 bearer = os.getenv("BLACKMARBLE_TOKEN")

    # If using Google Colab, use Secrets
from google.colab import userdata

bearer = userdata.get('BLACKMARBLE_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#

Define region of interest for which we will download nighttime lights data. We obtain the polygon from GADM.

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.

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", bearer=bearer
)
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/41)
    AlgorithmType:                     b'SCI'
    DataResolution:                    b'Moderate'
    DayNightFlag:                      b'Day'
    EastBoundingCoord:                 0.0
    EndTime:                           b'2021-02-05 23:59:59.000'
    GranuleDayNightFlag:               b'Day'
    ...                                ...
    VersionID:                         b'001'
    VerticalTileNumber:                b'08'
    WestBoundingCoord:                 -10.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0
Hide code cell source 
fig, ax = plt.subplots(figsize=(16, 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/b15372f7ee7f70b0441512d606f179d240aaf0b16a9ca5bfccd5b96941db6e29.png

Monthly#

# Monthly data: raster for October 2021
VNP46A3_202110 = bm_raster(
    gdf, product_id="VNP46A3", date_range="2021-10-01", bearer=bearer
)
VNP46A3_202110
Hide code cell output 
<xarray.Dataset> Size: 13MB
Dimensions:                        (x: 1071, y: 1545, time: 1)
Coordinates:
  * x                              (x) float64 9kB -3.259 -3.255 ... 1.193 1.197
  * y                              (y) float64 12kB 11.17 11.17 ... 4.746 4.742
  * 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/51)
    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'
    ...                                ...
    VersionID:                         b'001'
    VerticalTileNumber:                b'07'
    WestBoundingCoord:                 -10.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0
Hide code cell source 
fig, ax = plt.subplots(figsize=(16, 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/4bd7194e044dc486cdafd3a5ad4c390285c8e00293dd78d2b9c6c5dc01aa1630.png

Annual#

# Annual data: raster for 2021
VNP46A4_2021 = bm_raster(
    gdf, product_id="VNP46A4", date_range="2021-01-01", bearer=bearer
)
VNP46A4_2021
Hide code cell output 
<xarray.Dataset> Size: 13MB
Dimensions:                        (x: 1071, y: 1545, time: 1)
Coordinates:
  * x                              (x) float64 9kB -3.259 -3.255 ... 1.193 1.197
  * y                              (y) float64 12kB 11.17 11.17 ... 4.746 4.742
  * 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/51)
    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'
    ...                                ...
    VersionID:                         b'001'
    VerticalTileNumber:                b'07'
    WestBoundingCoord:                 0.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0
Hide code cell source 
fig, ax = plt.subplots(figsize=(16, 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/b38ae4f17112cc9198ce2cf362b644491127dc0fffa474058b7ec75d503e4992.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"),
    bearer=bearer,
)

<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/31)
    AlgorithmType:                     b'SCI'
    DataResolution:                    b'Moderate'
    DayNightFlag:                      b'Day'
    EastBoundingCoord:                 10.0
    GranuleDayNightFlag:               b'Day'
    GRingPointLatitude:                [ 0. 10. 10.  0.]
    ...                                ...
    VersionID:                         b'001'
    VerticalTileNumber:                b'08'
    WestBoundingCoord:                 0.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0
Hide code cell source 
# Create the figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# 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/8516df7842b1b8064d31c7fac4adfd6ed99c166574332ee3eeb8fd742702a62b.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"),
    bearer=bearer,
)

<xarray.Dataset> Size: 159MB
Dimensions:                        (x: 1071, y: 1545, time: 12)
Coordinates:
  * x                              (x) float64 9kB -3.259 -3.255 ... 1.193 1.197
  * y                              (y) float64 12kB 11.17 11.17 ... 4.746 4.742
  * 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/39)
    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'
    ...                                ...
    VersionID:                         b'001'
    VerticalTileNumber:                b'08'
    WestBoundingCoord:                 -10.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0
Hide code cell source 
# Create the figure and axis
fig, ax = plt.subplots(figsize=(10, 6))

# 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/a17068b998eaa64aa372b6d4d49462287b64a81174681cf10f658c709f277973.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"),
    bearer=bearer,
)

<xarray.Dataset> Size: 53MB
Dimensions:                        (x: 1071, y: 1545, time: 4)
Coordinates:
  * x                              (x) float64 9kB -3.259 -3.255 ... 1.193 1.197
  * y                              (y) float64 12kB 11.17 11.17 ... 4.746 4.742
  * 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/39)
    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'18'
    ...                                ...
    NorthBoundingCoord:                10.0
    publisher_email:                   b'modis-ops@lists.nasa.gov'
    publisher_name:                    b'LAADS'
    publisher_url:                     b'https://ladsweb.modaps.eosdis.nasa.gov'
    SouthBoundingCoord:                0.0
    WestBoundingCoord:                 0.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/ba867cd10bd90ac737fd775b5d87bd6fecb47316f27daf3d77b8f1654f42d36d.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.#

Downloading and Saving Black Marble Data Locally#

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.

bm_raster(
    gdf,
    product_id="VNP46A4",
    date_range=pd.date_range("2022-01-01", "2022-01-01", freq="YS"),
    output_directory="data/",
    output_skip_if_exists=True,  # default
    bearer=bearer,
)
Hide code cell output 
<xarray.Dataset> Size: 13MB
Dimensions:                        (x: 1071, y: 1545, time: 1)
Coordinates:
  * x                              (x) float64 9kB -3.259 -3.255 ... 1.193 1.197
  * y                              (y) float64 12kB 11.17 11.17 ... 4.746 4.742
  * time                           (time) datetime64[ns] 8B 2022-01-01
Data variables:
    NearNadir_Composite_Snow_Free  (time, y, x) float64 13MB nan nan ... nan nan
Attributes: (12/51)
    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'
    ...                                ...
    VersionID:                         b'001'
    VerticalTileNumber:                b'08'
    WestBoundingCoord:                 0.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0

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

bm_raster(
    gdf,
    product_id="VNP46A4",
    date_range=pd.date_range("2022-01-01", "2022-01-01", freq="YS"),
    output_directory="data/",
    output_skip_if_exists=False,
    bearer=bearer,
)
Hide code cell output 
<xarray.Dataset> Size: 13MB
Dimensions:                        (x: 1071, y: 1545, time: 1)
Coordinates:
  * x                              (x) float64 9kB -3.259 -3.255 ... 1.193 1.197
  * y                              (y) float64 12kB 11.17 11.17 ... 4.746 4.742
  * time                           (time) datetime64[ns] 8B 2022-01-01
Data variables:
    NearNadir_Composite_Snow_Free  (time, y, x) float64 13MB nan nan ... nan nan
Attributes: (12/51)
    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'
    ...                                ...
    VersionID:                         b'001'
    VerticalTileNumber:                b'07'
    WestBoundingCoord:                 -10.0
    AREA_OR_POINT:                     Area
    scale_factor:                      1.0
    add_offset:                        0.0

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, 8))

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/bd554e80da364e9be348d4a0c59958c082c84fa5906d9025ae3a6861394fdb54.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.#

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