Night light analysis.

aoi_bbox = (11.016088, 51.062939, 15.334464, 53.700949)
aoi_geom = [box(*aoi_bbox).__geo_interface__]

with rasterio.open(SRC_PATH) as src:
  clipped, clipped_transform = mask(src, aoi_geom, crop=True)
  clipped_meta = src.meta.copy()
  clipped_meta.update({
      "height": clipped.shape[1],
      "width": clipped.shape[2],
      "transform": clipped_transform,
  })

with rasterio.open(AOI_CLIPPED) as src:
    data = src.read(1)
    bounds = src.bounds

fig, ax = plt.subplots(figsize=(10, 8))
im = ax.imshow(
    np.log1p(np.clip(data, 0, None)),
    extent=(bounds.left, bounds.right, bounds.bottom, bounds.top),
    cmap="inferno",
)

vals = aoi_data[~np.isnan(aoi_data)] if aoi_nodata is None else aoi_data[aoi_data != aoi_nodata]
# We can now get vals.min(), .max(), .mean(), .std()...
# percentiles with `float(np.percentile(vals, 50))` for P50, etc.

# Radiance is power-law-distributed: use log-spaced bins on log-x so the
# city-scale tail is actually visible alongside the dim-pixel mass.
log_bins = np.logspace(np.log10(vals.min()), np.log10(vals.max()), 80)
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].hist(vals, bins=log_bins, color="steelblue")
axes[0].set_xscale("log")

axes[1].hist(vals, bins=log_bins, color="steelblue")
axes[1].set_xscale("log")
axes[1].set_yscale("log")
plt.show()

(Simplified the code above, but you get the idea.)

from rasterio.features import shapes as rio_shapes

def count_polygons(mask_arr, transform):
    n = 0
    for _, val in rio_shapes(mask_arr, transform=transform):
        if val == 1:
            n += 1
    return n

def otsu_threshold(x):
    # Standard 1-D Otsu on a 256-bin histogram
    hist, edges = np.histogram(x, bins=256)
    p = hist / hist.sum()
    omega = np.cumsum(p)
    mu = np.cumsum(p * (edges[:-1] + np.diff(edges) / 2))
    mu_t = mu[-1]
    denom = omega * (1 - omega)
    sigma_b2 = np.where(denom > 0, (mu_t * omega - mu) ** 2 / denom, 0)
    k = np.argmax(sigma_b2)
    return float((edges[k] + edges[k + 1]) / 2)

otsu = otsu_threshold(vals)
print(f"Otsu threshold: {otsu:.3f}")

with rasterio.open(AOI_CLIPPED) as src:
    aoi_transform = src.transform

thresholds = [NOISE_FLOOR, 1.0, 2.0, 5.0, 10.0, stats["P95"], stats["P99"], otsu]
thresholds = sorted(set(round(t, 3) for t in thresholds))

rows = []
for t in thresholds:
    mask_arr = (
        (aoi_data > t)
        & (aoi_data > NOISE_FLOOR)
        & (aoi_data != aoi_nodata)
    ).astype(np.uint8)
    lit_px = int(mask_arr.sum())
    n_polys = count_polygons(mask_arr, aoi_transform)
    rows.append((t, lit_px, lit_px / vals.size * 100, n_polys))

print(f"{'thresh':>8} {'lit_px':>10} {'lit_%':>8} {'n_polys':>8}")
for t, lp, pct, n in rows:
    print(f"{t:>8.3f} {lp:>10,} {pct:>7.2f}% {n:>8}")

# we then use matplotlib to plot `rows`

def polygonize_night_lights(raster_path, threshold=1.0, output_path=LIT_AREAS, noise_floor=0.0):
  """Convert lit pixels (radiance > max(threshold, noise_floor)) to polygons.

  threshold and noise_floor are in nW/cm^2/sr.
  """
  features = []
  with rasterio.open(raster_path) as src:
      data = src.read(1)
      nodata = src.nodata
      cut = max(threshold, noise_floor)
      mask_arr = np.where(
          (data > cut) & (data != nodata),
          1,
          0,
      ).astype(np.uint8)
      for geom, val in shapes(mask_arr, transform=src.transform):
          if val == 1:
              features.append({"geometry": shape(geom), "value": val})
      crs = src.crs

  gdf = gpd.GeoDataFrame(features, crs=crs)
  gdf = gdf.dissolve()
  gdf.to_file(output_path, driver="GPKG")
  print(f"Saved {len(gdf)} polygons to {output_path}")
  return gdf

gdf = polygonize_night_lights(AOI_CLIPPED, threshold=4.0)
gdf.head()

# Explode the dissolved multipolygon back into individual polygons.
polys = gdf.explode(index_parts=False).reset_index(drop=True)

# Reproject to ETRS89 LAEA Europe (EPSG:3035) for equal-area calculation.
polys_eq = polys.to_crs("EPSG:3035")
polys["area_km2"] = polys_eq.geometry.area / 1e6

n = len(polys)
areas = polys["area_km2"].to_numpy()
print(f"# polygons: {n}")
print(f"min:    {areas.min():.4f} km²")
print(f"max:    {areas.max():.2f} km²")
print(f"mean:   {areas.mean():.2f} km²")
print(f"median: {np.median(areas):.4f} km²")
print(f"total:  {areas.sum():.1f} km² lit")
for q in [50, 75, 90, 95, 99]:
    print(f"P{q}:    {np.percentile(areas, q):.3f} km²")

print("\nTop 10 largest polygons:")
print(polys.nlargest(10, "area_km2")[["area_km2"]].to_string())

from datetime import datetime, timezone
import math

ACQ_TIME = datetime(2025, 12, 22, 22, 17, 30, tzinfo=timezone.utc)
_ref_new_moon = datetime(2000, 1, 6, 18, 14, tzinfo=timezone.utc)
_synodic = 29.530588853
_age = ((ACQ_TIME - _ref_new_moon).total_seconds() / 86400.0) % _synodic
_illum = (1 - math.cos(2 * math.pi * _age / _synodic)) / 2
print(f"Moon age:    {_age:.2f} days post-new")
print(f"Illuminated: {_illum * 100:.1f}%  (lower → less moonlight scatter, more emitter-only signal)")

Out:

Moon age:    2.85 days post-new
Illuminated: 8.9%  (lower → less moonlight scatter, more emitter-only signal)

import json
import urllib.request
from rasterio.windows import from_bounds
from rasterio.vrt import WarpedVRT
from rasterio.enums import Resampling
from rasterio.features import geometry_mask
from scipy.ndimage import gaussian_filter
from shapely.geometry import shape as shp_shape

STAC_ITEM_URL = (
    "https://stac.dataspace.copernicus.eu/v1/collections/"
    "sentinel-2-night-time-acquisitions/items/"
    "S2A_OPER_PRD_R058_V20251222T221730_20251222T221745"
)
stac_item = json.loads(urllib.request.urlopen(STAC_ITEM_URL).read())

east_polygon = {
    "type": "Polygon",
    "coordinates": stac_item["geometry"]["coordinates"][1],
}
east_shape = shp_shape(east_polygon)
EAST_BBOX = east_shape.bounds  # (minx, miny, maxx, maxy)
print(f"East strip bbox: ({EAST_BBOX[0]:.4f}, {EAST_BBOX[1]:.4f}, {EAST_BBOX[2]:.4f}, {EAST_BBOX[3]:.4f})")

DECIM = 10  # 10m → 100m pixels
S2_REF_URL = band_vsi("B03")

def read_decim_on_ref(band, decim=DECIM):
    with rasterio.open(S2_REF_URL) as ref, rasterio.open(band_vsi(band)) as src:
        with WarpedVRT(
            src,
            crs=ref.crs,
            transform=ref.transform,
            width=ref.width,
            height=ref.height,
            resampling=Resampling.bilinear,
        ) as vrt:
            win = from_bounds(*EAST_BBOX, transform=vrt.transform)
            out_shape = (int(win.height // decim), int(win.width // decim))
            d = vrt.read(1, window=win, out_shape=out_shape,
                         out_dtype=np.float32, resampling=Resampling.average)
            nd = vrt.nodata
            win_transform = vrt.window_transform(win)
            scale_x = win.width / out_shape[1]
            scale_y = win.height / out_shape[0]
            decim_transform = win_transform * win_transform.scale(scale_x, scale_y)
    valid = (d != nd) & np.isfinite(d) & (np.abs(d) < 1000)
    return np.where(valid, d, 0.0), decim_transform

print("Reading decimated Sentinel-2 bands over east strip...")
b12_e, s2_grid_transform = read_decim_on_ref("B12")
b08_e, _ = read_decim_on_ref("B08")
b03_e, _ = read_decim_on_ref("B03")
print(f"Sentinel-2 decimated shape: {b03_e.shape}  (~100m pixels)")

def local_column_destripe(arr, mask, height=51):
    """Subtract a moving column-median computed over `height` rows.

    The Sen2VM detector-overlap stripes run along the orbit ground track
    (~15° off vertical here). A tall narrow window captures the local
    per-column stripe offset without needing to rotate into orbit coordinates.
    """
    from numpy.lib.stride_tricks import sliding_window_view
    arr_m = np.where(mask, arr, np.nan).astype(np.float32)
    H = arr_m.shape[0]
    half = height // 2
    bg = np.full_like(arr_m, np.nan)
    if H >= height:
        for j in range(arr_m.shape[1]):
            padded = np.concatenate([np.full(half, np.nan), arr_m[:, j], np.full(half, np.nan)])
            windows = sliding_window_view(padded, height)
            with np.errstate(all="ignore"):
                bg[:, j] = np.nanmedian(windows, axis=1)
    bg = np.where(np.isnan(bg), 0.0, bg)
    return np.where(mask, arr - bg, 0.0)

footprint_mask = geometry_mask(
    [east_polygon],
    out_shape=b03_e.shape,
    transform=s2_grid_transform,
    invert=True,  # True inside the polygon
)
print(f"Polygon covers {100 * footprint_mask.sum() / footprint_mask.size:.1f}% of bbox")

print("Destriping (Sen2VM detector-overlap noise)...")
b12_e = local_column_destripe(b12_e, footprint_mask, height=51)
b08_e = local_column_destripe(b08_e, footprint_mask, height=51)
b03_e = local_column_destripe(b03_e, footprint_mask, height=51)

lit_e = (
    np.where(b12_e > 0.020, b12_e, 0.0) / 0.1
    + np.where(b08_e > 0.060, b08_e, 0.0) / 0.5
    + np.where(b03_e > 0.100, b03_e, 0.0) / 0.5
)
lit_e = np.where(footprint_mask, lit_e, 0.0)
n_lit = int((lit_e > 0).sum())
print(f"Lit pixels in east strip: {n_lit:,} / {int(footprint_mask.sum()):,} ({100 * n_lit / max(footprint_mask.sum(), 1):.3f}% of swath)")

sigma_px_e = 750 / 100 / 2.355
s2_blurred_e = gaussian_filter(lit_e, sigma=sigma_px_e)
s2_blurred_e = np.where(footprint_mask, s2_blurred_e, np.nan)