Is light pollution / electrification changing in this region?

Nighttime lights are exactly what they sound like: a satellite looking straight down at Earth after dark, measuring how much light shines up from the ground. More light usually means more people, more buildings, more electricity. A newly lit village, a growing city, a factory running at night. When the lights go out, that’s often a power cut. Watch the same patch of ground year after year and you can see electrification spread, cities sprawl, and disasters knock the power out.

Which data to use, and how

Use this dataset: Black Marble VNP46A2 (short name VNP46A2), which covers 2012 → today at about ~500 m per pixel. It’s the easiest to work with because NASA has already cleaned it up for you. They’ve removed the glow of the moon and filled gaps left by clouds (this is the “BRDF-corrected, gap-filled” version). There’s also a rawer version, Black Marble VNP46A1 (VNP46A1, at-sensor radiance), but skip it. It hasn’t been cleaned, so beginners should stay on VNP46A2.

Then do four steps (the runnable code further down does all of this — this is just the idea):

1. Pick your region (a small latitude/longitude box around it) and the years you care about.
2. Decide what counts as “lit”. Each pixel has a brightness number; mark a pixel as lit if it’s above a threshold (e.g. brightness over 5 nW/cm²/sr, a standard cutoff).
3. One clean picture per year. For each year, take the middle (median) brightness of each pixel across all the nights. The median throws out odd bright nights (a passing flare, leftover moonlight) and gives you one stable map per year.
4. Count and compare. Count how many pixels are “lit” each year to get an electrification timeline, then subtract one year’s map from another (or a disaster month from the month before) to see exactly where lights appeared or vanished.

That’s the whole method. Below, “How a scientist answers this” names the exact tools, and the Run it cell does it live on synthetic data so you can watch each step work.

What you can find out

• Cities growing. The bright footprint spreading outward over the years.
• Rural electrification. A pixel that lights up where there was darkness before.
• Power outages after disasters. The classic Black Marble result was Puerto Rico going dark after Hurricane Maria in 2017. Comparing the disaster month to the month before draws a map of who lost power.
• Light pollution as a wildlife clue. Brighter skies are a rough proxy (stand-in) for how much artificial light is disturbing nocturnal animals and ecosystems.
• Industrial activity. Gas flaring (burning off waste gas at oil and gas sites) shows up as intense bright spots.

What it can’t tell you

• What kind of light it is. The satellite sees brightness, not its source. On its own it can’t tell a housing block from a factory from a gas flare. To separate them you need extra maps (land use, infrastructure layers).
• Outages faster than the satellite flies over. VIIRS passes over roughly once per night (~12-hour revisit), so a flicker that comes and goes within hours is invisible. Sub-hourly outage tracking needs commercial systems, not this.
• The color of the light. VIIRS measures one broad band of brightness (it’s “broadband panchromatic”), so it can’t break light into colors — questions about blue-light pollution specifically are out of reach here.

Gotchas to watch for

• Use the cleaned version, not the raw one. Analyze VNP46A2 (gap-filled, BRDF-corrected), not raw VNP46A1. The raw product still has moonlight and cloud gaps baked in, which will fool you.
• The moon is a lamp too. Moonlight reflecting off the ground can masquerade as city lights. VNP46A2 already removes most of it, but each pixel ships with a quality flag (QF). Check it and drop pixels the flag marks as bad.
• Snow and ice fake new lights. At high latitudes in winter, bright white snow can reflect light and suddenly look like a town that wasn’t there. The fix: compare like-for-like seasons (or mask out snowy months) so a snowfall isn’t mistaken for electrification.
• Gas flares blow out the picture. Flares are so intense they can saturate the detector (max it out) and dominate the brightness statistics over oilfields. If your region has them, expect them to drown out subtler changes nearby.
• Don’t naïvely glue the old record onto the new one. The pre-2012 lights record (DMSP-OLS, 1992–2013) came from a different, coarser sensor that doesn’t line up with VIIRS. Stacking them directly creates fake jumps. If you truly need pre-2012 history, use a published intercalibrated (harmonized) product instead.

The real-data code

The Run it cell above runs this method on synthetic data with no login. Below is the whole recipe against the real archive: search, open the daily lights, build one median map per year, count the lit pixels, difference the first and last year, and plot it. It needs a free Earthdata Login.

import earthaccess, numpy as np, xarray as xr
import matplotlib.pyplot as plt

earthaccess.login(strategy="netrc")

aoi = (75.5, 18.5, 76.5, 19.5)   # your region as (W, S, E, N) — here, rural Maharashtra
years = range(2013, 2024)
LIT = 5.0                        # "lit" threshold in nW/cm2/sr (a standard cutoff)
NTL = "Gap_Filled_DNB_BRDF-Corrected_NTL"   # the cleaned, gap-filled lights layer
QF  = "Mandatory_Quality_Flag"              # 0/1 = good, 2 = poor/gap — drop those

# 1. One clean brightness map per year: open every VNP46A2 granule for the year,
#    keep good-quality pixels, and take the per-pixel median across all nights.
#    The median rejects odd bright nights (leftover moonlight, a passing flare).
def yearly_median(y):
    grans = earthaccess.search_data(
        short_name="VNP46A2", bounding_box=aoi,
        temporal=(f"{y}-01-01", f"{y}-12-31"))
    nights = []
    for fh in earthaccess.open(grans):
        ds = xr.open_dataset(fh, group="HDFEOS/GRIDS/VNP_Grid_DNB/Data Fields")
        ntl = ds[NTL].where(ds[QF] <= 1)        # mask poor / gap-filled-uncertain pixels
        ntl = ntl * 0.1                          # apply scale factor -> nW/cm2/sr
        nights.append(ntl)
    if not nights:
        return None
    stack = xr.concat(nights, dim="night")
    return stack.median(dim="night")

maps = {y: m for y in years if (m := yearly_median(y)) is not None}

# 2. Electrification timeline: count "lit" pixels (> threshold) in each yearly map.
lit_pixels = np.array([int((maps[y] > LIT).sum()) for y in maps])
yrs = np.array(list(maps.keys()))
growth = 100.0 * (lit_pixels[-1] - lit_pixels[0]) / max(lit_pixels[0], 1)
trend = "brightening / electrifying" if growth > 0 else "going dark"
print(f"Lit area is {trend}: {lit_pixels[0]} lit px in {yrs[0]} -> "
      f"{lit_pixels[-1]} in {yrs[-1]} ({growth:+.0f}%)")

# 3. Where the lights appeared/vanished: difference the last and first yearly map.
diff = maps[yrs[-1]] - maps[yrs[0]]

# 4. Plot the electrification timeline and the change map side by side.
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4))
ax1.plot(yrs, lit_pixels, "-o"); ax1.set_xlabel("year")
ax1.set_ylabel(f"lit pixels (> {LIT} nW/cm2/sr)"); ax1.set_title("Electrification timeline")
im = diff.plot(ax=ax2, cmap="RdBu_r", center=0, add_colorbar=True)
ax2.set_title(f"Change in lights {yrs[0]}->{yrs[-1]} (red=brighter)")
plt.tight_layout(); plt.show()

# Before you trust it: replay the Puerto Rico / Hurricane Maria case (median Sep-2017 vs
# Sep-2016) as a known-answer check, and watch for snow or gas flares faking new lights.

Where the data comes from

It’s all NASA, through one free login (Earthdata Login). Black Marble is produced at NASA Goddard and served from the LAADS DAAC archive, so a single sign-in gets you everything. No juggling separate accounts.

Sources

• Black Marble: https://blackmarble.gsfc.nasa.gov/
• VNP46A2 product: https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A2/
• Black Marble after-Maria study: https://earthobservatory.nasa.gov/images/91047/puerto-rico-power-outage