Ready.

# @title 1 · Choose the district and build a common 30 m grid
DISTRICT = "In Ouzzal (Hoggar) — Algeria"  # @param ["Mahd adh Dhahab — Saudi Arabia", "Bou Azzer — Morocco", "Sohar–Lasail — Oman", "In Ouzzal (Hoggar) — Algeria"]
D = {d["label"]: d for d in fetch_json("districts.json")}[DISTRICT]
p = f"districts/{D['id']}"
aster, a_tr, a_crs = fetch_raster(f"{p}/aster_indices.tif")      # grid of reference: ASTER 30 m, UTM
H, W = aster.shape[1:]
def onto(rel, band, resampling=Resampling.average):
    arr, tr, crs = fetch_raster(rel)
    out = np.full((H, W), np.nan, np.float32)
    reproject(arr[band].astype(np.float32), out, src_transform=tr, src_crs=crs, src_nodata=np.nan,
              dst_transform=a_tr, dst_crs=a_crs, dst_nodata=np.nan, resampling=resampling)
    return out
aloh = aster[0].astype(np.float32)
ferric = onto(f"{p}/s2_indices.tif", 0)
items = list(Client.open("https://earth-search.aws.element84.com/v1").search(collections=["cop-dem-glo-30"], bbox=D["bbox"]).items())
from odc.geo.geobox import GeoBox
gbox = GeoBox((H, W), a_tr, a_crs)
dem = load(items, bands=["data"], geobox=gbox, resampling="bilinear")["data"].isel(time=0).values.astype(np.float32)
known = fetch_vector(f"{p}/known_sites.geojson").to_crs(a_crs)
mined = features.rasterize([(g, 1) for g in known.geometry], out_shape=(H, W), transform=a_tr, fill=0).astype(bool)
valid = np.isfinite(aloh) & np.isfinite(ferric) & np.isfinite(dem)
print(f"{D['name']}: grid {W} × {H} at 30 m · {len(known)} mapped mining areas ({mined.sum() * 0.0009:.1f} km²)")

In Ouzzal (Amesmessa–Tirek), Hoggar: grid 2087 × 3328 at 30 m · 48 mapped mining areas (7.1 km²)

# @title 2 · Structures from the elevation model
def hillshade(z, az, alt=30, cell=30):
    gy, gx = np.gradient(z, cell)
    slope, aspect = np.arctan(np.hypot(gx, gy)), np.arctan2(-gx, gy)
    a, e = np.radians(360 - az + 90), np.radians(alt)
    return np.clip(np.sin(e) * np.cos(slope) + np.cos(e) * np.sin(slope) * np.cos(a - aspect), 0, 1)
z = np.nan_to_num(dem, nan=np.nanmedian(dem))
lines = []
for az in (0, 45, 90, 135):
    lines += probabilistic_hough_line(canny(hillshade(z, az), sigma=2.5, low_threshold=0.05, high_threshold=0.15),
                                      threshold=15, line_length=20, line_gap=3, rng=0)
L = np.array(lines, float)
strike = (np.degrees(np.arctan2(L[:, 1, 0] - L[:, 0, 0], L[:, 0, 1] - L[:, 1, 1])) + 360) % 180
def raster_lines(sel):
    m = np.zeros((H, W), bool)
    for (x0, y0), (x1, y1) in L[sel].astype(int):
        n = max(abs(x1 - x0), abs(y1 - y0)) + 1
        xs, ys = np.linspace(x0, x1, n).astype(int), np.linspace(y0, y1, n).astype(int)
        ok = (xs >= 0) & (xs < W) & (ys >= 0) & (ys < H)
        m[ys[ok], xs[ok]] = True
    return m
all_lin = raster_lines(np.ones(len(L), bool))
ns = (strike < 20) | (strike > 160)                       # roughly north–south structures
ns_lin = raster_lines(ns)
d_all = ndimage.distance_transform_edt(~all_lin) * 0.03   # km
d_ns = ndimage.distance_transform_edt(~ns_lin) * 0.03
print(f"{len(L):,} lineament segments; {ns.sum():,} strike within 20° of north–south")

428 lineament segments; 279 strike within 20° of north–south

# @title 3 · Build the favourability map
def rising(x):
    v = x[valid]
    lo, hi = np.percentile(v, 50), np.percentile(v, 98)
    return np.clip((x - lo) / (hi - lo), 0, 1)
def near(d, km=2.0):
    return np.clip(1 - d / km, 0, 1)
WEIGHTS = {"Al-OH alteration": 1.0, "ferric iron": 0.6, "near any structure": 0.6, "near N–S structure": 1.0}   # edit
mu = {"Al-OH alteration": rising(aloh), "ferric iron": rising(ferric),
      "near any structure": near(d_all), "near N–S structure": near(d_ns)}
w = {k: v / max(WEIGHTS.values()) for k, v in WEIGHTS.items()}
layers = [np.clip(mu[k] * w[k], 1e-6, 1) for k in mu]
prod = np.prod(layers, axis=0)
alg_sum = 1 - np.prod([1 - l for l in layers], axis=0)
GAMMA = 0.9
fav = (alg_sum ** GAMMA) * (prod ** (1 - GAMMA))
fav[~valid] = np.nan
shade = hillshade(z, 315)
p80, p995 = np.nanpercentile(fav, [80, 99.5])
fig, ax = plt.subplots(1, 2, figsize=(15, 6.8))
for a in ax:
    a.imshow(shade, cmap="gray", vmin=0, vmax=1)
im = ax[0].imshow(np.where(fav >= p80, fav, np.nan), cmap="inferno", vmin=p80, vmax=p995, alpha=0.9)
ax[0].contour(mined.astype(float), levels=[0.5], colors="cyan", linewidths=0.9)
ax[0].set_title("Most favourable 20 %\n(cyan = mapped mining areas)"); ax[0].axis("off")
fig.colorbar(im, ax=ax[0], shrink=0.7, label="favourability")
top10 = fav >= np.nanpercentile(fav, 90)
ax[1].imshow(np.where(top10, 1, np.nan), cmap=ListedColormap(["#F59E0B"]), interpolation="nearest", alpha=0.9)
ax[1].contour(mined.astype(float), levels=[0.5], colors="#0E7C86", linewidths=0.9)
ax[1].set_title("Top 10 % most favourable ground\n(teal = mapped mining areas)"); ax[1].axis("off")
plt.tight_layout(); plt.show()

# @title 4 · Success-rate curve
trapezoid = getattr(np, "trapezoid", None) or getattr(np, "trapz")   # name differs across NumPy versions
f = fav[valid]; t = mined[valid]
if t.sum() < 30:
    print("Too few mapped mining areas in this district for a meaningful validation — the map is shown for illustration only.")
else:
    order = np.argsort(-f)
    captured = np.cumsum(t[order]) / t.sum()
    area = np.arange(1, len(f) + 1) / len(f)
    auc = trapezoid(captured, area)
    k10 = captured[int(0.1 * len(f))]
    fig = go.Figure()
    fig.add_trace(go.Scatter(x=100 * area[::200], y=100 * captured[::200], name="favourability map", line=dict(color="#B45309", width=3)))
    fig.add_trace(go.Scatter(x=[0, 100], y=[0, 100], name="random", line=dict(color="#9CA3AF", dash="dash")))
    fig.update_layout(title=f"Success rate — AUC {auc:.2f}; the top 10 % of the ranking captures {100 * k10:.0f} % of mapped mining ground",
                      xaxis_title="% of district, from most favourable", yaxis_title="% of mapped mining ground captured",
                      height=440, **LAYOUT)
    fig.show()
    # Which evidence carries the result? Drop one layer at a time and re-score.
    rows = []
    for k in mu:
        keep = [np.clip(mu[j] * w[j], 1e-6, 1) for j in mu if j != k]
        fk = ((1 - np.prod([1 - l for l in keep], axis=0)) ** GAMMA) * (np.prod(keep, axis=0) ** (1 - GAMMA))
        ff = fk[valid]; o = np.argsort(-ff)
        rows.append({"left out": k, "AUC without it": round(float(trapezoid(np.cumsum(t[o]) / t.sum(), area)), 3)})
    display(pd.DataFrame(rows).set_index("left out"))
    key = min(rows, key=lambda r: r["AUC without it"])
    print(f"Most important evidence here: {key['left out']} — without it the AUC falls from {auc:.2f} to {key['AUC without it']:.2f}.")
    hurt = [r["left out"] for r in rows if r["AUC without it"] > auc + 0.01]
    if hurt:
        print("Not helping in this district: " + ", ".join(hurt) + " — the AUC rises when it is left out, so consider "
              "lowering its weight in WEIGHTS (and ask why: different deposit style, cover, or weathering).")

Most important evidence here: Al-OH alteration — without it the AUC falls from 0.73 to 0.61.
Not helping in this district: ferric iron — the AUC rises when it is left out, so consider lowering its weight in WEIGHTS (and ask why: different deposit style, cover, or weathering).

Evidence	Why it matters for orogenic gold	Data
Al-OH alteration (sericite)	fluids altered the rock around gold-bearing veins	ASTER 2000–2008
Ferric iron	weathered sulphides leave iron-oxide stains	Sentinel-2
Proximity to structures	shear zones channel the fluids	Copernicus DEM lineaments
Structural direction	the main shear zones here strike roughly north–south	Copernicus DEM lineaments

From evidence to targets — a transparent prospectivity map¶

Take-home notebook · Workshop 2 · IC-TWEP 2026 · QuestFeed¶

Fuzzy evidence and combination¶

Honest validation: a success-rate curve¶

What the map does — and does not — tell you¶

	AUC without it
left out
Al-OH alteration	0.615
ferric iron	0.788
near any structure	0.711
near N–S structure	0.683