Ctrl+K

Examples using BlackMarblePy

Examples using BlackMarblePy#

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 colorcet as cc
import contextily as cx
import geopandas
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns

from blackmarble.extract import bm_extract
from blackmarble.raster import bm_raster

%load_ext autoreload
%autoreload 2

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

Create NASA Earthdata Bearer Token#

BlackMarblePy requires using NASA Earthdata bearer token. To obtain a token, please follow the steps below:

  1. Go to the NASA LAADS Archive

  2. Click “Login” (bottom on top right); create an account if needed.

  3. Click “See wget Download Command” (bottom near top, in the middle)

  4. After clicking, you will see text that can be used to download data. The “Bearer” token will be a long string in red.

Tip

After logging in, the below will show the bearer token in red instead of INSERT_DOWNLOAD_TOKEN_HERE. Sometimes, after logging in, the NASA website will redirect to another part of the website. To obtain the bearer token, just navigate to the NASA LAADS Archive after logging in.

../_images/nasa_laads_login.png

Define NASA bearer token#

For instructions on obtaining a NASA bearer token, please see above.

bearer = os.getenv("BLACKMARBLE_TOKEN")

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()
Make this Notebook Trusted to load map: File -> Trust Notebook

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
r_20210205 = bm_raster(
    gdf, product_id="VNP46A2", date_range="2021-02-05", bearer=bearer
)

<xarray.Dataset>
Dimensions:                            (x: 1071, y: 1544, time: 1)
Coordinates:
  * x                                  (x) float64 -3.26 -3.256 ... 1.194 1.198
  * y                                  (y) float64 11.17 11.17 ... 4.748 4.744
  * time                               (time) datetime64[ns] 2021-02-05
Data variables:
    Gap_Filled_DNB_BRDF-Corrected_NTL  (time, y, x) float64 nan nan ... nan nan
Attributes: (12/41)
    AlgorithmType:                     b'SCI'
    AREA_OR_POINT:                     Area
    DataResolution:                    b'Moderate'
    DayNightFlag:                      b'Day'
    EastBoundingCoord:                 0.0
    EndTime:                           b'2021-02-05 23:59:59.000'
    ...                                ...
    TileID:                            b'61017008'
    VersionID:                         b'001'
    VerticalTileNumber:                b'08'
    WestBoundingCoord:                 -10.0
    scale_factor:                      1.0
    add_offset:                        0.0
Hide code cell source 
fig, ax = plt.subplots()

r_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/1b41deae12e7ef541074da05caf34d57c7adca1052f3091ff8dab6b0511fe190.png

Monthly#

# Monthly data: raster for October 2021
r_202110 = bm_raster(gdf, product_id="VNP46A3", date_range="2021-10-01", bearer=bearer)
Hide code cell source 
fig, ax = plt.subplots()

r_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/6374abb13455caced55f4e2fd90be7d4dd609104210249d5f993e2c042b33449.png

Annual#

### Annual data: raster for 2021
r_2021 = bm_raster(gdf, product_id="VNP46A4", date_range="2021-01-01", bearer=bearer)
Hide code cell output
Hide code cell source 
fig, ax = plt.subplots()

r_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/de55de609202a1231d1e057437ecca4b3f0c5a3a7ccb31575accd83fc2a4069c.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
date_range = pd.date_range("2022-01-01", "2022-03-31", freq="D")

r_daily = bm_raster(
    gdf,
    product_id="VNP46A2",
    date_range=date_range,
    bearer=bearer,
)
Hide code cell source 
fig, ax = plt.subplots()

r_daily["Gap_Filled_DNB_BRDF-Corrected_NTL"].mean(dim=["x", "y"]).plot(ax=ax)

ax.text(
    0,
    -0.2,
    "Source: NASA Black Marble VNP46A2",
    ha="left",
    va="center",
    transform=ax.transAxes,
    fontsize=10,
    color="black",
    weight="normal",
)
ax.set_title("Ghana: NTL Radiance", fontsize=20);

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

#### Raster stack of annual data
r_annual = bm_raster(
    gdf,
    product_id="VNP46A4",
    date_range=pd.date_range("2019-01-01", "2022-01-01", freq="YS"),
    bearer=bearer,
)
Hide code cell source 
# Set up the figure and subplots
fig, axs = plt.subplots(1, 4, figsize=(20, 6))

for i, t in enumerate(r_annual["time"]):
    ax = axs[i]
    r_annual["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()

Visualizing difference in radiance year over year#

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

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

delta = (
    (
        (
            r_annual["NearNadir_Composite_Snow_Free"].sel(time="2022-01-01")
            - r_annual["NearNadir_Composite_Snow_Free"].sel(time="2019-01-01")
        )
        / r_annual["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.Positron)

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=20);

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.

Contents