Source code: Sentinel-1 ship detection in Python: a CA-CFAR detector off the Port of Rotterdam

"""Ship detection in Sentinel-1 SAR imagery with a CA-CFAR detector, off the Port of Rotterdam.

Searches the Microsoft Planetary Computer for Sentinel-1 RTC scenes from one orbit, streams VV and
VH backscatter for the Maasmond approaches and anchorages onto a 10 m grid, masks land with ESA
WorldCover, runs a cell-averaging constant false alarm rate (CA-CFAR) detector written with
NumPy and SciPy, clusters the detections into targets, separates fixed structures (wind turbines,
platforms) from vessels by checking which targets persist across six dates, flags possible
azimuth ambiguities and writes the vessels to GeoJSON.
Run with:  python ship_detection_cfar.py
"""
import json
import platform
from importlib.metadata import version
from pathlib import Path

import geopandas as gpd
import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import planetary_computer
import pystac_client
import rasterio
import requests
from pyproj import Transformer
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.vrt import WarpedVRT
from scipy import ndimage, stats
from scipy.sparse import coo_matrix
from scipy.sparse.csgraph import connected_components
from scipy.spatial import cKDTree
from shapely.geometry import Point, box, shape
from shapely.ops import unary_union

HERE = Path(__file__).parent
OUT, DATA = HERE / "output", HERE / "data"
for d in (OUT, DATA):
    d.mkdir(exist_ok=True)
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
pd.set_option("display.width", 130)

print("=== versions ===")
print(f"python {platform.python_version()} | pystac-client {pystac_client.__version__} | "
      f"planetary-computer {version('planetary-computer')} | rasterio {rasterio.__version__} | "
      f"numpy {np.__version__} | scipy {version('scipy')} | geopandas {gpd.__version__}")


def section(title: str) -> None:
    print(f"\n=== {title} ===")


STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"
BBOX = (3.55, 51.85, 4.25, 52.30)   # lon/lat: Maasmond, the Rotterdam anchorages and Hollandse Kust Zuid
ORBIT = 110                         # descending; covers the whole box every 6 days (Sentinel-1C and 1D)
WINDOW = "2026-09-01/2026-10-16"
MAIN_DATE = "2026-10-11"
RES, CRS = 10, "EPSG:32631"         # 10 m pixels in UTM zone 31N, the native grid of these scenes
PFA = 1e-6                          # design probability of false alarm per pixel
GUARD, BACKGROUND = 41, 101         # CFAR window sizes in pixels (410 m guard, 1,010 m background)
MIN_PIXELS = 4                      # smallest target kept: 4 pixels = 400 m²
MERGE = 250                         # blobs closer than this (m) are one target
FIXED_M = 100                       # within this distance on another pass = the same object
AMB_KM = (4.0, 6.5)                 # azimuth ambiguity offset band for IW mode

to_utm = Transformer.from_crs("EPSG:4326", CRS, always_xy=True)
to_lonlat = Transformer.from_crs(CRS, "EPSG:4326", always_xy=True)
xs, ys = to_utm.transform([BBOX[0], BBOX[2], BBOX[0], BBOX[2]], [BBOX[1], BBOX[1], BBOX[3], BBOX[3]])
x0, x1 = np.floor(min(xs) / RES) * RES, np.ceil(max(xs) / RES) * RES
y0, y1 = np.floor(min(ys) / RES) * RES, np.ceil(max(ys) / RES) * RES
WIDTH, HEIGHT = int((x1 - x0) / RES), int((y1 - y0) / RES)
TRANSFORM = from_origin(x0, y1, RES, RES)


def warp(href: str, resampling=Resampling.nearest, nodata=np.nan, dtype="float32") -> np.ndarray:
    """Stream a remote COG onto the 10 m analysis grid."""
    with rasterio.open(href) as src, WarpedVRT(src, crs=CRS, transform=TRANSFORM, width=WIDTH, height=HEIGHT,
                                               resampling=resampling, src_nodata=src.nodata, nodata=nodata,
                                               dtype=dtype) as vrt:
        return vrt.read(1)


def to_db(power):
    with np.errstate(divide="ignore", invalid="ignore"):
        return 10 * np.log10(np.where(power > 0, power, np.nan))


# --- Step 1: find the scenes ---------------------------------------------------------------
catalog = pystac_client.Client.open(STAC_URL, modifier=planetary_computer.sign_inplace)
items = catalog.search(collections=["sentinel-1-rtc"], bbox=BBOX, datetime=WINDOW,
                       query={"sat:relative_orbit": {"eq": ORBIT}}).item_collection()
aoi = box(*BBOX)
by_date = {}
for it in items:
    by_date.setdefault(it.datetime.strftime("%Y-%m-%d"), []).append(it)
cover = {d: unary_union([shape(i.geometry) for i in its]).intersection(aoi).area / aoi.area for d, its in by_date.items()}
DATES = sorted(d for d, c in cover.items() if c > 0.99)
section("STAC search")
print(f"{len(items)} scene slices from relative orbit {ORBIT}; dates covering the whole box: {', '.join(DATES)}")
main = by_date[MAIN_DATE][0]
p = main.properties
print(f"{MAIN_DATE}: {', '.join(i.id for i in by_date[MAIN_DATE])}")
print(f"platform {p['platform']} | {p['sat:orbit_state']} | polarisations {p['sar:polarizations']} | "
      f"ENL {p['sar:looks_equivalent_number']} | {p['proj:code']}")
ENL = p["sar:looks_equivalent_number"]


def read_date(date: str, pol: str) -> np.ndarray:
    cache = DATA / f"{date}_{pol}.npy"
    if cache.exists():
        return np.load(cache)
    img = np.full((HEIGHT, WIDTH), np.nan, "float32")
    for it in by_date[date]:
        part = warp(it.assets[pol].href)
        img = np.where(np.isnan(img), part, img)
    np.save(cache, img)
    return img


vv, vh = read_date(MAIN_DATE, "vv"), read_date(MAIN_DATE, "vh")

# --- Step 2: land mask from ESA WorldCover ----------------------------------------------------
wc = catalog.search(collections=["esa-worldcover"], bbox=BBOX,
                    query={"esa_worldcover:product_version": {"eq": "2.0.0"}}).item_collection()
landcover = np.zeros((HEIGHT, WIDTH), "uint8")
for it in wc:
    part = warp(it.assets["map"].href, nodata=0, dtype="uint8")
    landcover = np.where(landcover == 0, part, landcover)
land = (landcover != 0) & (landcover != 80)          # 80 = permanent water, 0 = no data (open sea)
land = ndimage.binary_dilation(land, iterations=25)  # 250 m buffer: quays, breakwaters, moored ships
lab, n = ndimage.label(~land & np.isfinite(vv))
sizes = ndimage.sum(np.ones_like(lab), lab, index=np.arange(1, n + 1))
sea = lab == 1 + int(np.argmax(sizes))               # keep the open sea, drop harbour basins and lakes
PIXEL_KM2 = RES * RES / 1e6
section("land mask")
print(f"grid {WIDTH} x {HEIGHT} pixels at {RES} m ({WIDTH * RES / 1000:.1f} x {HEIGHT * RES / 1000:.1f} km)")
print(f"WorldCover land {((landcover != 0) & (landcover != 80)).mean():.1%}; "
      f"open sea after a 250 m coastal buffer: {sea.sum() * PIXEL_KM2:,.0f} km² ({sea.mean():.1%})")


# --- Step 3: CA-CFAR --------------------------------------------------------------------------
def ca_cfar(img: np.ndarray, sea: np.ndarray, pfa: float = PFA, enl: float = ENL):
    """Cell-averaging CFAR. The clutter mean comes from a ring of sea pixels between the guard
    and background windows; the multiplier comes from a gamma model of multilook speckle."""
    x = np.where(sea, np.nan_to_num(img), 0).astype("float64")
    m = sea.astype("float64")
    box_sum = lambda a, k: ndimage.uniform_filter(a, k, mode="constant") * k * k
    n_ring = box_sum(m, BACKGROUND) - box_sum(m, GUARD)
    with np.errstate(divide="ignore", invalid="ignore"):
        clutter = (box_sum(x, BACKGROUND) - box_sum(x, GUARD)) / n_ring
    factor = stats.gamma.isf(pfa, a=enl, scale=1 / enl)   # P(I > factor * mean) = pfa for gamma clutter
    ok = sea & (n_ring > 0.5 * (BACKGROUND ** 2 - GUARD ** 2))
    return ok & (img > factor * clutter), clutter, factor


det_vv, clutter_vv, factor = ca_cfar(vv, sea)
det_vh, clutter_vh, _ = ca_cfar(vh, sea)
section("CA-CFAR")
print(f"Pfa {PFA:g}, ENL {ENL}: threshold = {factor:.2f} x local clutter mean ({10 * np.log10(factor):.1f} dB above it)")
print(f"background ring {GUARD * RES} m to {BACKGROUND * RES} m wide; "
      f"expected false-alarm pixels over pure gamma speckle: {PFA * sea.sum():.0f}")
print(f"detected pixels: VV {det_vv.sum():,}, VH {det_vh.sum():,}")
for name, img, det in (("VV", vv, det_vv), ("VH", vh, det_vh)):
    lab, n = ndimage.label(det)
    lone = np.isin(lab, 1 + np.flatnonzero(ndimage.sum(det, lab, np.arange(1, n + 1)) < MIN_PIXELS))
    print(f"{name}: sea clutter median {np.nanmedian(to_db(img[sea & ~det])):.1f} dB; isolated speckle hits "
          f"(clusters < {MIN_PIXELS} px) {lone.sum():,} = {lone.sum() / sea.sum():.1e} of sea pixels")


# --- Step 4: cluster detections into targets -------------------------------------------------------
def targets(det_a, det_b, img_a, img_b, clut_a, clut_b) -> pd.DataFrame:
    """Label connected detections in either polarisation, then merge blobs closer than MERGE metres."""
    det = det_a | det_b
    lab, n = ndimage.label(det, structure=np.ones((3, 3)))
    idx = np.arange(1, n + 1)
    rows, cols = np.indices(det.shape)
    ratio = np.where(det, np.nan_to_num(img_a / clut_a) + np.nan_to_num(img_b / clut_b), 0)
    tot = ndimage.sum(ratio, lab, idx)
    sl = ndimage.find_objects(lab)
    blobs = pd.DataFrame({
        "row": ndimage.sum(ratio * rows, lab, idx) / tot, "col": ndimage.sum(ratio * cols, lab, idx) / tot,
        "pixels": ndimage.sum(det, lab, idx), "in_a": ndimage.maximum(det_a, lab, idx).astype(bool),
        "in_b": ndimage.maximum(det_b, lab, idx).astype(bool),
        "extent_m": [RES * np.hypot(s[0].stop - s[0].start, s[1].stop - s[1].start) for s in sl],
        "tcr_a": to_db(ndimage.maximum(np.where(det, img_a / clut_a, 0), lab, idx)),
        "tcr_b": to_db(ndimage.maximum(np.where(det, img_b / clut_b, 0), lab, idx)),
        "peak_a": to_db(ndimage.maximum(np.where(det, img_a, 0), lab, idx)),
        "peak_b": to_db(ndimage.maximum(np.where(det, img_b, 0), lab, idx))})
    links = cKDTree(blobs[["row", "col"]] * RES).query_pairs(MERGE, output_type="ndarray")
    graph = coo_matrix((np.ones(len(links)), (links[:, 0], links[:, 1])), shape=(n, n))
    blobs["target"] = connected_components(graph, directed=False)[1]
    blobs["score"] = blobs[["tcr_a", "tcr_b"]].max(axis=1)
    first = blobs.sort_values("score", ascending=False).groupby("target")
    df = first.first()
    df["pixels"] = first.pixels.sum()
    df["blobs"] = first.size()
    df["in_a"], df["in_b"] = first.in_a.any(), first.in_b.any()
    df = df[df.pixels >= MIN_PIXELS].reset_index(drop=True)
    df["x"], df["y"] = x0 + (df.col + 0.5) * RES, y1 - (df.row + 0.5) * RES
    return df.rename(columns={"in_a": "in_vv", "in_b": "in_vh", "tcr_a": "tcr_vv_db", "tcr_b": "tcr_vh_db",
                              "peak_a": "peak_vv_db", "peak_b": "peak_vh_db"}).drop(columns=["score"])


t = targets(det_vv, det_vh, vv, vh, clutter_vv, clutter_vh)
section("targets")
n_blobs = ndimage.label(det_vv | det_vh, structure=np.ones((3, 3)))[1]
print(f"{n_blobs:,} connected blobs in VV or VH -> {len(t)} targets after merging blobs within {MERGE} m "
      f"and dropping those under {MIN_PIXELS} pixels")
print(f"seen in both {(t.in_vv & t.in_vh).sum()}, VV only {(t.in_vv & ~t.in_vh).sum()}, "
      f"VH only {(~t.in_vv & t.in_vh).sum()}")
both = t[t.in_vv & t.in_vh]
print(f"median target-to-clutter ratio of targets seen in both: VV {both.tcr_vv_db.median():.1f} dB, "
      f"VH {both.tcr_vh_db.median():.1f} dB")


# --- Step 5: sidelobes and azimuth ambiguities -------------------------------------------------------
def arm_angle(img: np.ndarray, r: float, c: float) -> float:
    """Bearing (degrees clockwise from grid north) of the brightest straight arm through pixel (r, c)."""
    db = np.nan_to_num(to_db(img), nan=-40)
    s = np.arange(150, 2500, RES) / RES                  # sample 150 m to 2.5 km either side, in pixels
    angles = np.arange(0, 180, 0.25)
    score = [ndimage.map_coordinates(db, [np.r_[r - s * np.cos(a), r + s * np.cos(a)],
                                          np.r_[c + s * np.sin(a), c - s * np.sin(a)]], order=1).mean()
             for a in np.radians(angles)]
    return float(angles[int(np.argmax(score))])


bright = t.loc[t.tcr_vh_db.idxmax()]
a1 = arm_angle(vh, bright.row, bright.col)
AZ = min((a1, (a1 + 90) % 180), key=lambda a: min(a, 180 - a))   # the arm closest to north-south
u_az = np.array([np.sin(np.radians(AZ)), np.cos(np.radians(AZ))])
u_rg = np.array([u_az[1], -u_az[0]])
section("sidelobes and ambiguities")
print(f"brightest target: VH {bright.tcr_vh_db:.1f} dB above clutter at row {bright.row:.0f}, col {bright.col:.0f}")
print(f"its sidelobe cross gives azimuth {AZ:.1f}° and range {(AZ + 90) % 180:.1f}° (grid bearings)")
print(f"azimuth ambiguity offset lambda*R*PRF/(2V) for IW: {AMB_KM[0]}-{AMB_KM[1]} km")


def flag_artefacts(df: pd.DataFrame) -> pd.DataFrame:
    xy = df[["x", "y"]].to_numpy()
    pk = df[["peak_vv_db", "peak_vh_db"]].to_numpy()
    side, amb = np.zeros(len(df), bool), np.zeros(len(df), bool)
    for k in range(len(df)):
        d = xy - xy[k]
        along, across = np.abs(d @ u_az), np.abs(d @ u_rg)
        brighter = np.any(pk > pk[k] + 10, axis=1)          # 10 dB brighter in VV or VH
        on_arm = ((along < 3000) & (across < 60)) | ((across < 3000) & (along < 60))
        side[k] = np.any(brighter & on_arm)
        amb[k] = np.any(brighter & (along > AMB_KM[0] * 1000) & (along < AMB_KM[1] * 1000) & (across < 300))
    return df.assign(sidelobe=side, ambiguity=amb & ~side)


t = flag_artefacts(t)
print(f"targets on the sidelobe cross of one >= 10 dB brighter: {t.sidelobe.sum()}")
print(f"targets at the ambiguity offset of one >= 10 dB brighter: {t.ambiguity.sum()}")

# --- Step 6: fixed structures persist, ships move ----------------------------------------------------
per_date = {MAIN_DATE: t}
for d in DATES:
    if d not in per_date:
        a, b = read_date(d, "vv"), read_date(d, "vh")
        ok = sea & np.isfinite(a) & np.isfinite(b)
        (da, ca, _), (db_, cb, _) = ca_cfar(a, ok), ca_cfar(b, ok)
        per_date[d] = flag_artefacts(targets(da, db_, a, b, ca, cb))
others = [d for d in DATES if d != MAIN_DATE]
t["dates_seen"] = sum((cKDTree(per_date[d][["x", "y"]]).query(t[["x", "y"]], distance_upper_bound=FIXED_M)[0]
                       < FIXED_M).astype(int) for d in others)
t["fixed"] = t.dates_seen >= len(others) - 1        # there on (almost) every other pass
t["class"] = np.select([t.sidelobe, t.ambiguity, t.fixed], ["sidelobe", "possible ambiguity", "fixed structure"],
                       "vessel")
section("persistence across dates")
print(f"{len(others)} other passes; main-date targets also found within {FIXED_M} m on n other passes: "
      f"{t.dates_seen.value_counts().sort_index().to_dict()}")
print(t["class"].value_counts().to_string())

turbines = None
try:      # optional sanity check: offshore wind turbines mapped in OpenStreetMap inside the box
    cache = DATA / "osm_turbines.json"
    if not cache.exists():
        q = (f"[out:json][timeout:60];node[\"generator:source\"=\"wind\"]"
             f"({BBOX[1]},{BBOX[0]},{BBOX[3]},{BBOX[2]});out;")
        r = requests.post("https://overpass-api.de/api/interpreter", data={"data": q}, timeout=90,
                          headers={"User-Agent": "TechnicalGIS tutorial (https://technicalgis.com)"})
        r.raise_for_status()
        cache.write_text(r.text, encoding="utf-8")
    nodes = json.loads(cache.read_text(encoding="utf-8"))["elements"]
    tx, ty = to_utm.transform([e["lon"] for e in nodes], [e["lat"] for e in nodes])
    tc, tr = (np.array(tx) - x0) / RES, (y1 - np.array(ty)) / RES
    keep = (tc >= 0) & (tc < WIDTH) & (tr >= 0) & (tr < HEIGHT)
    keep[keep] = sea[tr[keep].astype(int), tc[keep].astype(int)]
    turbines = pd.DataFrame({"x": np.array(tx)[keep], "y": np.array(ty)[keep]})
    dist = cKDTree(t[["x", "y"]]).query(turbines[["x", "y"]])[0]
    fx = t[t["class"] == "fixed structure"]
    dist_fx = cKDTree(turbines[["x", "y"]]).query(fx[["x", "y"]])[0]
    section("OpenStreetMap check")
    print(f"{len(turbines)} OSM wind turbines over open sea in the box; {np.mean(dist < 100):.0%} have a "
          f"target within 100 m (median offset {np.median(dist):.0f} m)")
    print(f"fixed structures within 100 m of an OSM turbine: {np.sum(dist_fx < 100)} of {len(fx)}")
except Exception as exc:
    print(f"OpenStreetMap check skipped: {exc}")

# --- Step 7: outputs ---------------------------------------------------------------------------------
lon, lat = to_lonlat.transform(t.x.to_numpy(), t.y.to_numpy())
cols = ["class", "pixels", "blobs", "extent_m", "tcr_vv_db", "tcr_vh_db", "in_vv", "in_vh", "dates_seen"]
gdf = gpd.GeoDataFrame(t[cols].round(1), geometry=gpd.points_from_xy(np.round(lon, 5), np.round(lat, 5)),
                       crs="EPSG:4326")
gdf.to_file(OUT / f"ship_detections_{MAIN_DATE}.geojson", driver="GeoJSON")
ships = gdf[gdf["class"] == "vessel"]
section(f"result {MAIN_DATE}")
print(f"{len(ships)} vessels written to output/ship_detections_{MAIN_DATE}.geojson (all {len(gdf)} targets kept, "
      f"with a class column)")
print(f"vessel extent median {ships.extent_m.median():.0f} m, largest {ships.extent_m.max():.0f} m; "
      f"seen in VV {ships.in_vv.mean():.0%}, in VH {ships.in_vh.mean():.0%}, in both {(ships.in_vv & ships.in_vh).mean():.0%}")
print(ships.sort_values("tcr_vv_db", ascending=False).head(5)
      .assign(lon=ships.geometry.x, lat=ships.geometry.y).drop(columns="geometry").to_string(index=False))

fixed_tree = cKDTree(t.loc[t.fixed & ~t.sidelobe, ["x", "y"]])
rows_out = []
for d, td in per_date.items():
    near_fixed = fixed_tree.query(td[["x", "y"]], distance_upper_bound=FIXED_M)[0] < FIXED_M
    art = td.sidelobe | td.ambiguity
    rows_out.append({"date": d, "platform": by_date[d][0].properties["platform"], "targets": len(td),
                     "fixed": int((near_fixed & ~art).sum()), "artefacts": int(art.sum()),
                     "vessels": int((~near_fixed & ~art).sum())})
counts = pd.DataFrame(rows_out).sort_values("date")
section("vessel counts per pass")
print(counts.to_string(index=False))
print(f"mean {counts.vessels.mean():.0f} vessels per pass, range {counts.vessels.min()}-{counts.vessels.max()}")

# --- Step 8: figures -----------------------------------------------------------------------------------
EXTENT = [0, WIDTH * RES / 1000, 0, HEIGHT * RES / 1000]           # km from the lower-left corner
km = lambda df: ((df.col + 0.5) * RES / 1000, HEIGHT * RES / 1000 - (df.row + 0.5) * RES / 1000)
STYLE = {"vessel": ("o", "#e8590c", 22), "fixed structure": ("s", "#1d6fa3", 10),
         "sidelobe": ("x", "#868e96", 14), "possible ambiguity": ("^", "#9c36b5", 22)}


def save(fig, name):
    fig.savefig(OUT / name, facecolor="white")
    plt.close(fig)
    print(f"saved output/{name}")


section("figures")
fig, axs = plt.subplots(1, 2, figsize=(11, 5.6), sharey=True)
for ax, img, name, lim in ((axs[0], vv, "VV", (-22, 0)), (axs[1], vh, "VH", (-32, -10))):
    ax.imshow(to_db(img), cmap="gray", vmin=lim[0], vmax=lim[1], extent=EXTENT)
    ax.contour(np.flipud(sea), levels=[0.5], colors="#e8590c", linewidths=0.7, extent=EXTENT, origin="lower")
    ax.set_title(f"Sentinel-1 {name} backscatter (dB), {MAIN_DATE}")
    ax.set_xlabel("km east")
axs[0].set_ylabel("km north")
axs[1].annotate("Hollandse Kust Zuid\nwind farm", (17, 46), color="white", fontsize=8)
axs[1].annotate("Maasvlakte /\nPort of Rotterdam", (37, 2), color="black", fontsize=8)
fig.suptitle("Orange line: open-sea mask (ESA WorldCover water, 250 m coastal buffer)", y=0.04, fontsize=8.5)
save(fig, "01-scene-vv-vh-land-mask.png")

ship = t[t["class"] == "vessel"].sort_values("tcr_vv_db").iloc[len(t[t["class"] == "vessel"]) // 2]
r0, c0 = int(ship.row), int(ship.col)
h = 70
fig, axs = plt.subplots(1, 2, figsize=(11, 4.6), gridspec_kw={"width_ratios": [1, 1.25]})
axs[0].imshow(to_db(vv[r0 - h:r0 + h, c0 - h:c0 + h]), cmap="gray", vmin=-22, vmax=5,
              extent=[-h * RES, h * RES, -h * RES, h * RES])
for size, color, label in ((GUARD, "#1d6fa3", "guard window"), (BACKGROUND, "#e8590c", "background window")):
    half = size * RES / 2
    axs[0].add_patch(plt.Rectangle((-half, -half), 2 * half, 2 * half, fill=False, color=color, lw=1.4, label=label))
axs[0].plot(0, 0, "+", color="yellow", ms=10, label="cell under test")
axs[0].legend(loc="lower left", fontsize=7.5)
axs[0].set_title(f"CA-CFAR windows around a typical vessel (VV)")
axs[0].set_xlabel("metres")
z = (vv / clutter_vv)[sea]
z = z[np.isfinite(z)]
bins = np.logspace(-1.5, 2.5, 110)
axs[1].hist(z, bins=bins, density=True, color="#adb5bd", label="sea pixels, VV (I / local clutter mean)")
xx = np.logspace(-1.5, 2.5, 400)
axs[1].plot(xx, stats.gamma.pdf(xx, a=ENL, scale=1 / ENL), color="#1d6fa3", lw=1.5, label=f"gamma model, ENL {ENL}")
axs[1].axvline(factor, color="#e8590c", lw=1.5, label=f"threshold {factor:.2f} (Pfa {PFA:g})")
axs[1].set_xscale("log")
axs[1].set_yscale("log")
axs[1].set_ylim(1e-8, 5)
axs[1].set_xlabel("intensity / local clutter mean")
axs[1].set_ylabel("probability density")
axs[1].set_title("Speckle model and threshold")
axs[1].legend(fontsize=7.5, loc="lower left")
save(fig, "02-cfar-window-and-threshold.png")

sel = t[t["class"] == "vessel"]
h = 250                                              # 5 x 5 km window on the busiest patch of sea
busy, re, ce = np.histogram2d(sel.row, sel.col, bins=[np.arange(0, HEIGHT, 2 * h), np.arange(0, WIDTH, 2 * h)])
i, j = np.unravel_index(np.argmax(busy), busy.shape)
r0, c0 = int(re[i] + h), int(ce[j] + h)
fig, axs = plt.subplots(1, 2, figsize=(11, 5.4), sharey=True)
win = (slice(r0 - h, r0 + h), slice(c0 - h, c0 + h))
ext = [(c0 - h) * RES / 1000, (c0 + h) * RES / 1000, HEIGHT * RES / 1000 - (r0 + h) * RES / 1000,
       HEIGHT * RES / 1000 - (r0 - h) * RES / 1000]
for ax, img, det, name, lim in ((axs[0], vv, det_vv, "VV", (-22, 0)), (axs[1], vh, det_vh, "VH", (-32, -10))):
    ax.imshow(to_db(img[win]), cmap="gray", vmin=lim[0], vmax=lim[1], extent=ext)
    ax.imshow(np.ma.masked_where(~det[win], det[win]), cmap="autumn", alpha=0.9, extent=ext, interpolation="nearest")
    ax.set_title(f"{name}: {det[win].sum():,} detected pixels in this 5 x 5 km window")
    ax.set_xlabel("km east")
axs[0].set_ylabel("km north")
save(fig, "03-vv-vs-vh-anchorage.png")

fig, ax = plt.subplots(figsize=(8.6, 8.4))
ax.imshow(to_db(vv), cmap="gray", vmin=-24, vmax=2, extent=EXTENT, alpha=0.85)
for cls, (m, color, s) in STYLE.items():
    sub = t[t["class"] == cls]
    ax.scatter(*km(sub), marker=m, s=s, facecolor="none" if m in "os^" else color, edgecolor=color,
               color=color if m == "x" else None, lw=1, label=f"{cls} ({len(sub)})")
if turbines is not None:
    ax.scatter((turbines.x - x0) / 1000, (turbines.y - y0) / 1000, s=2,
               color="#74c0fc", label=f"OSM wind turbines ({len(turbines)})")
ax.legend(loc="lower left", fontsize=8)
ax.set_title(f"CA-CFAR targets off Rotterdam, Sentinel-1 VV + VH, {MAIN_DATE}")
ax.set_xlabel("km east")
ax.set_ylabel("km north")
save(fig, "04-detections-map.png")

r0, c0, h = int(bright.row), int(bright.col), 260
win = (slice(max(r0 - h, 0), r0 + h), slice(c0 - h, min(c0 + h, WIDTH)))
fig, ax = plt.subplots(figsize=(7.4, 6.4))
ax.imshow(to_db(vh[win]), cmap="gray", vmin=-32, vmax=-5, extent=[win[1].start * RES / 1000, win[1].stop * RES / 1000,
          HEIGHT * RES / 1000 - win[0].stop * RES / 1000, HEIGHT * RES / 1000 - win[0].start * RES / 1000])
near = t[(t.row.between(win[0].start, win[0].stop)) & (t.col.between(win[1].start, win[1].stop))]
for cls, (m, color, s) in STYLE.items():
    sub = near[near["class"] == cls]
    ax.scatter(*km(sub), marker=m, s=s * 2, facecolor="none" if m in "os^" else color, edgecolor=color,
               color=color if m == "x" else None, lw=1.2, label=f"{cls} ({len(sub)})")
bx, by = km(bright.to_frame().T)
for ang, name in ((AZ, "azimuth"), (AZ + 90, "range")):
    dx, dy = 2.4 * np.sin(np.radians(ang)), 2.4 * np.cos(np.radians(ang))
    ax.plot([float(bx.iloc[0]) - dx, float(bx.iloc[0]) + dx], [float(by.iloc[0]) - dy, float(by.iloc[0]) + dy],
            ":", color="#ffd43b", lw=0.8, alpha=0.6)
    ax.annotate(name, (float(bx.iloc[0]) + dx, float(by.iloc[0]) + dy), color="#ffd43b", fontsize=8)
ax.legend(loc="lower left", fontsize=7.5)
ax.set_title("Brightest target (VH): sidelobe cross and nearby turbines")
ax.set_xlabel("km east")
ax.set_ylabel("km north")
save(fig, "05-sidelobes-and-turbines.png")

fig, axs = plt.subplots(1, 2, figsize=(11, 3.8))
x = np.arange(len(counts))
axs[0].bar(x, counts.vessels, color="#e8590c", label="vessels")
axs[0].bar(x, counts.fixed, bottom=counts.vessels, color="#1d6fa3", label="fixed structures")
axs[0].bar(x, counts.artefacts, bottom=counts.vessels + counts.fixed, color="#adb5bd", label="sidelobes / ambiguities")
axs[0].set_xticks(x, [f"{d[5:]}\n{p[-2:].upper()}" for d, p in zip(counts.date, counts.platform)])
axs[0].set_ylabel("targets")
axs[0].set_title("Targets per pass, relative orbit 110")
axs[0].legend(fontsize=7.5)
v = t[t["class"] == "vessel"]
axs[1].hist(v.extent_m, bins=np.arange(0, 560, 20), color="#e8590c", label="vessels")
axs[1].hist(t[t["class"] == "fixed structure"].extent_m, bins=np.arange(0, 560, 20), color="#1d6fa3", alpha=0.6,
            label="fixed structures")
axs[1].set_xlabel("brightest blob extent (bounding-box diagonal, m)")
axs[1].set_ylabel("count")
axs[1].set_title(f"Target sizes, {MAIN_DATE}")
axs[1].legend(fontsize=7.5)
save(fig, "06-counts-and-sizes.png")