After the storm, which neighborhoods lost power — and who lives there?

At night, a city glows: streetlights, signs, windows. A satellite overhead measures how bright each square kilometre is. That number is called radiance (how much light is coming up toward space). When a hurricane knocks out the grid, the city goes dark and the brightness number plunges. Compare a clear night before the storm to a night after, and you can watch the blackout from orbit. Then comes the human question: roughly how many people were sitting in the dark.

Verified locally. Over Tampa Bay (27.6–28.1 °N), brightness averaged 19.18 nW/cm²/sr on the clear night of 8 Oct 2024, the eve of Hurricane Milton’s landfall. By the night of 11 Oct it had fallen to 9.35 nW/cm²/sr, a 51% drop. The lights then crept back (13.8 on the 12th, 15.3 on the 13th) as crews restored power. That recovery curve is the grid coming back online, written in light.

Which data to use, and how

Use this dataset: Black Marble (VNP46A2) (short name VNP46A2), NASA’s nightly nighttime-lights map. It comes already cleaned up: it fills cloud gaps and corrects for moonlight, so each pixel is a fair, comparable measure of how much light the ground threw off that night. It covers 2012 → today at roughly 500 m per pixel.

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

1. Pick your area: a small latitude/longitude box around the city or coastline — and the dates: a clear night before landfall and one or more nights after.
2. One number per night. Average the brightness inside your box for each night. A clean “before” night gives you the baseline; the “after” nights show the drop.
3. Map the drop. Because the grid is ~500 m, you can compare each pixel before vs. after and see exactly which neighborhoods lost the most light.
4. Add the people. Lay a free population grid, WorldPop (wpic1km), over the pixels that went dark and add up the population underneath. To name the places, clip your box against geoBoundaries county outlines (the ADM2 level) so you can say “Pinellas and Hillsborough,” not “a dark blob.”

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 see each step work.

What you can find out

• How much darker a city got after the storm: the before-vs-after drop in average brightness over your box (the verified Tampa Bay example was a 51% drop).
• Which neighborhoods went dark: at ~500 m resolution you can map the per-pixel drop and see which districts lost the most light.
• How fast the grid recovered: track the nightly average and watch the brightness climb back toward its pre-storm baseline as power is restored.
• Roughly how many people were in the dark zone: overlay WorldPop population on the pixels that lost light and sum the people underneath.
• Which counties to name: clip your box against geoBoundaries so your story has real place names.
• A before/after story for any storm since 2012: VIIRS has flown nightly since then, so the same recipe works for past hurricanes, ice storms, or wildfire evacuations.

What it can’t tell you

• A confirmed power outage. A darker pixel is consistent with a blackout, but clouds, smoke, storm surge, and people evacuating also dim a city. The satellite sees less light; it doesn’t know why. Cross-check utility outage maps before claiming the grid failed.
• Who, exactly, lost power. WorldPop is a modeled population estimate, not a customer list. It gives you the order of magnitude of people affected — thousands vs. millions — not addresses or households.
• Why one dark patch is dark. Because the product blends nights together to fill clouds, a single dark pixel can just mean missing data, not a real outage. Always check the quality flags.
• Damage in dollars or lives. Radiance is light, not cost or casualties. Impact assessment needs ground reporting, insurance data, and damage surveys.
• Houses dark vs. streetlights off. The satellite sees streetlights, signs, and window-spill all mixed together; it can’t separate them.
• Who was hit hardest. WorldPop gives head counts, not income, age, or vulnerability. To ask who was least able to recover, you’d pair it with census data.

Gotchas to watch for

• Pick a genuinely clear “before” night. Your baseline is everything. If the night before was hazy or moonlit, your “drop” might be the weather, not the storm. Black Marble already corrects for moonlight, but a clean, cloud-free baseline still gives the most trustworthy comparison.
• The gap-filling can hide missing data. To paper over clouds, the product borrows brightness from nearby nights (it’s “gap-filled”). That usually helps, but it means a dark pixel can be missing data rather than a real outage. Each pixel comes with quality information. Keep the high-confidence pixels and treat the rest with suspicion.
• Read the real fill value and scale from the file, don’t assume. Satellite files mark “no data” with a special placeholder number (the _FillValue) and sometimes store brightness in compressed units you multiply back up (a scale_factor). Forget to drop the fill value and it poisons your average; ignore the scale and your numbers are off. On the granules checked here the fill was -999.9 and the scale was 1.0, but read them from the file every time.
• WorldPop and the lights are on different grids. WorldPop is ~1 km; the lights are ~500 m. Before you can sum population under the dark pixels, line the two grids up (reproject one onto the other). Otherwise you’re adding apples to a different map.

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 Black Marble, open the nightly tiles, average brightness per night, map the per-pixel drop, overlay WorldPop population to count who went dark, and plot it. It needs a free Earthdata Login.

Verified locally. VNP46A2 is a daily HDF5 tile on a 2400×2400 grid. The brightness lives at HDFEOS/GRIDS/VIIRS_Grid_DNB_2d/Data Fields/Gap_Filled_DNB_BRDF-Corrected_NTL, with lat/lon arrays in the same group so you can subset to your bounding box before averaging. Read the real _FillValue and scale_factor from the dataset attributes — on the granules checked here the fill was -999.9 and the scale was 1.0, with units nW/cm²/sr.

import os, re, warnings, collections, io
warnings.filterwarnings("ignore")
import earthaccess, h5py, requests, numpy as np
import rasterio
from rasterio.warp import reproject, Resampling
import matplotlib.pyplot as plt

# load Earthdata creds from .env without `source` (passwords can break the shell)
for line in open(".env"):
    m = re.match(r'\s*(?:export\s+)?([A-Z0-9_]+)\s*=\s*(.*)\s*$', line)
    if m: os.environ.setdefault(m.group(1), m.group(2).strip().strip('"').strip("'"))
earthaccess.login(strategy="environment")   # free Earthdata Login

W, S, E, N = -82.8, 27.6, -82.2, 28.1       # Tampa Bay (Hurricane Milton)
GRP = 'HDFEOS/GRIDS/VIIRS_Grid_DNB_2d/Data Fields'
NTL = f'{GRP}/Gap_Filled_DNB_BRDF-Corrected_NTL'

# 1. Open one VNP46A2 tile and return the brightness grid + lat/lon for the AOI.
def aoi_grid(granule):
    h = h5py.File(earthaccess.open([granule])[0], 'r')
    ds = h[NTL]
    lat, lon = h[f'{GRP}/lat'][:], h[f'{GRP}/lon'][:]
    fill  = float(np.array(ds.attrs.get('_FillValue', -999.9)).ravel()[0])
    scale = float(np.array(ds.attrs.get('scale_factor', 1.0)).ravel()[0])
    iy = np.where((lat >= S) & (lat <= N))[0]
    ix = np.where((lon >= W) & (lon <= E))[0]
    sub = ds[iy.min():iy.max()+1, ix.min():ix.max()+1].astype('float64')
    sub[sub == fill] = np.nan                 # drop the no-data placeholder
    return sub * scale, lat[iy], lon[ix]      # real brightness in nW/cm2/sr

# 2. Search Black Marble and group granules by night.
res = earthaccess.search_data(short_name="VNP46A2",
        temporal=("2024-10-05", "2024-10-13"), bounding_box=(W, S, E, N))
by_date = collections.OrderedDict()
for g in res:
    d = g['umm']['TemporalExtent']['RangeDateTime']['BeginningDateTime'][:10]
    by_date.setdefault(d, []).append(g)

# 3. One brightness number per night -> the before/after drop.
grids, nightly = {}, collections.OrderedDict()
for d, gs in by_date.items():
    grids[d], glat, glon = aoi_grid(gs[0])
    nightly[d] = np.nanmean(grids[d])
    print(d, round(nightly[d], 3), "nW/cm2/sr")
before, after = "2024-10-08", "2024-10-11"   # clear eve of landfall, post-Milton
drop = 1 - nightly[after] / nightly[before]
print(f"VERDICT: {drop*100:.0f}% drop in mean radiance after Hurricane Milton "
      f"({nightly[before]:.1f} -> {nightly[after]:.1f} nW/cm2/sr)")

# 4. Per-pixel drop, then count the people underneath the dark pixels.
delta = grids[before] - grids[after]
dark  = (delta > 0.5 * grids[before]) & np.isfinite(delta)   # lost >50% of its light

# WorldPop wpic1km (free, ~1 km) -> reproject onto the ~500 m lights grid, then sum.
url = ("https://data.worldpop.org/GIS/Population/Global_2000_2020_Constrained/"
       "2020/BSGM/USA/usa_ppp_2020_UNadj_constrained.tif")
with rasterio.open(url) as src:
    win = src.window(W, S, E, N)
    pop = src.read(1, window=win); pop[pop < 0] = 0
    pop_t = src.window_transform(win)
lights_t = rasterio.transform.from_bounds(glon.min(), glat.min(),
            glon.max(), glat.max(), dark.shape[1], dark.shape[0])
pop_on_lights = np.zeros(dark.shape, "float32")
reproject(pop, pop_on_lights, src_transform=pop_t, src_crs="EPSG:4326",
          dst_transform=lights_t, dst_crs="EPSG:4326", resampling=Resampling.sum)
in_dark = int(pop_on_lights[dark].sum())
print(f"VERDICT: ~{in_dark:,} people lived under pixels that lost over half their light")

# 5. Plot the recovery curve and the per-pixel drop map.
fig, (a1, a2) = plt.subplots(1, 2, figsize=(11, 4))
a1.plot(list(nightly), list(nightly.values()), "o-"); a1.tick_params(axis="x", rotation=45)
a1.axhline(nightly[before], ls="--", c="k", label="pre-storm baseline")
a1.set_ylabel("mean radiance (nW/cm2/sr)"); a1.legend(); a1.set_title("grid coming back online")
im = a2.imshow(delta, cmap="inferno"); a2.set_title("light lost 10/08 -> 10/11 (per pixel)")
fig.colorbar(im, ax=a2, label="Δ radiance"); plt.tight_layout(); plt.show()

# Before you trust it: a dark pixel is consistent with a blackout but also with cloud,
# evacuation, or gap-filled missing data — cross-check a utility outage map (see gotchas).

Where the data comes from

The lights come from NASA’s Black Marble product, which turns the VIIRS day/night band into a clean nightly map. You get it through one free Earthdata Login. The two extras are free and separate: WorldPop supplies the ~1 km population grid (who went dark), and geoBoundaries supplies the county outlines (so you can name the places). Stitch the three together and you go from “a dark blob” to “roughly N people in Pinellas and Hillsborough lost light, and the grid was ~half recovered by the 13th.”

Sources

• Black Marble VNP46A2 user guide: https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A2/
• Black Marble product suite: https://blackmarble.gsfc.nasa.gov/
• WorldPop population grids: https://www.worldpop.org/
• geoBoundaries administrative boundaries: https://www.geoboundaries.org/