Ready.

# @title 1 · Choose a district and load the elevation model
DISTRICT = "Mahd adh Dhahab — Saudi Arabia"  # @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]
bbox = D["bbox"]
lon0 = (bbox[0] + bbox[2]) / 2
crs = f"EPSG:{(32600 if bbox[1] >= 0 else 32700) + int((lon0 + 180) // 6) + 1}"
items = list(Client.open("https://earth-search.aws.element84.com/v1").search(collections=["cop-dem-glo-30"], bbox=bbox).items())
dem_ds = load(items, bands=["data"], bbox=bbox, crs=crs, resolution=30, resampling="bilinear")
dem_tr = dem_ds.odc.geobox.transform
dem = dem_ds["data"].isel(time=0).values.astype("float32")
dem[dem < -100] = np.nan
print(f"{D['name']}: elevation {np.nanmin(dem):.0f}–{np.nanmax(dem):.0f} m, grid {dem.shape[1]} × {dem.shape[0]} at 30 m")

Mahd adh Dhahab district, Arabian Shield: elevation 896–1425 m, grid 1722 × 1499 at 30 m

# @title 2 · Multi-directional hillshade
def hillshade(z, azimuth, altitude=30, cell=30):
    gy, gx = np.gradient(z, cell)
    slope = np.arctan(np.hypot(gx, gy))
    aspect = np.arctan2(-gx, gy)
    az, alt = np.radians(360 - azimuth + 90), np.radians(altitude)
    return np.clip(np.sin(alt) * np.cos(slope) + np.cos(alt) * np.sin(slope) * np.cos(az - aspect), 0, 1)

z = np.nan_to_num(dem, nan=np.nanmedian(dem))
shades = {az: hillshade(z, az) for az in (0, 45, 90, 135)}
multi = np.mean(list(shades.values()), axis=0)
fig, ax = plt.subplots(1, 2, figsize=(15, 6.5))
ax[0].imshow(dem, cmap="terrain"); ax[0].set_title("Copernicus DEM (m)")
ax[1].imshow(multi, cmap="gray"); ax[1].set_title("Multi-directional hillshade")
for a in ax:
    a.axis("off")
plt.tight_layout(); plt.show()

# @title 3 · Lineaments and their directions
segments = []
for az, hs in shades.items():
    edges = canny(hs, sigma=2.5, low_threshold=0.05, high_threshold=0.15)
    segments += probabilistic_hough_line(edges, threshold=15, line_length=20, line_gap=3, rng=0)
segments = np.array(segments, dtype=float)                    # (n, 2, 2): [(x0, y0), (x1, y1)] in pixels
dx, dy = segments[:, 1, 0] - segments[:, 0, 0], segments[:, 0, 1] - segments[:, 1, 1]
length_km = np.hypot(dx, dy) * 0.03
azimuth = (np.degrees(np.arctan2(dx, dy)) + 360) % 180       # 0 = N-S, 90 = E-W
print(f"{len(segments):,} lineament segments, total length {length_km.sum():,.0f} km")

fig = plt.figure(figsize=(15, 6.5))
a1 = fig.add_subplot(1, 2, 1)
a1.imshow(multi, cmap="gray")
for (x0, y0), (x1, y1) in segments:
    a1.plot([x0, x1], [y0, y1], color="#E4572E", lw=0.8)
a1.set_title("Automatically extracted lineaments"); a1.axis("off")
a2 = fig.add_subplot(1, 2, 2, projection="polar")
bins = np.radians(np.arange(0, 181, 10))
hist, _ = np.histogram(np.radians(azimuth), bins=bins, weights=length_km)
for h, b in zip(hist, bins[:-1]):
    for off in (0, np.pi):                                     # a rose diagram is symmetric
        a2.bar(b + off + np.radians(5), h, width=np.radians(10), color="#B45309", alpha=0.85, edgecolor="white")
a2.set_theta_zero_location("N"); a2.set_theta_direction(-1); a2.set_yticklabels([])
a2.set_title("Rose diagram (length-weighted)")
plt.tight_layout(); plt.show()
main = int(np.argmax(hist)) * 10
print(f"Dominant direction: {main}–{main + 10}° from north")

928 lineament segments, total length 724 km

Dominant direction: 60–70° from north

# @title 4 · Lineament density
mask = np.zeros(dem.shape, bool)
for (x0, y0), (x1, y1) in segments.astype(int):
    n = int(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 < mask.shape[1]) & (ys >= 0) & (ys < mask.shape[0])
    mask[ys[ok], xs[ok]] = True
# each lineament pixel stands for ~0.03 km of line; a Gaussian window (σ ≈ 0.75 km) gives km of lineament per km²
density = ndimage.gaussian_filter(mask.astype(float), sigma=25) / 0.03
fig, ax = plt.subplots(figsize=(9, 6.5))
im = ax.imshow(density, cmap="magma"); ax.set_title("Lineament density (km per km², smoothed over ~1.5 km)"); ax.axis("off")
fig.colorbar(im, ax=ax, shrink=0.75); plt.show()

# @title 5 · Alteration versus distance to lineaments
emit, e_tr, e_crs = fetch_raster(f"districts/{D['id']}/emit_minerals.tif")
cls = np.zeros(dem.shape, np.uint16)
reproject(emit[0], cls, src_transform=e_tr, src_crs=e_crs, dst_transform=dem_tr, dst_crs=crs, resampling=Resampling.nearest)
dist_km = ndimage.distance_transform_edt(~mask) * 0.03
altered = np.isin(cls, [4, 5, 6])          # kaolinite group, smectite clays, white mica / illite
known = cls > 0
near, far = known & (dist_km <= 0.5), known & (dist_km > 1.5)
rate_near, rate_far = altered[near].mean(), altered[far].mean()
print(f"Clay / white-mica share of identified pixels: {100 * rate_near:.1f} % within 500 m of lineaments, "
      f"{100 * rate_far:.1f} % beyond 1.5 km → enrichment ratio {rate_near / rate_far:.2f}")
bins_km = np.arange(0, 3.01, 0.25)
shares = [altered[known & (dist_km >= a) & (dist_km < b)].mean() * 100 for a, b in zip(bins_km[:-1], bins_km[1:])]
fig = go.Figure(go.Bar(x=[f"{a:.2f}–{b:.2f}" for a, b in zip(bins_km[:-1], bins_km[1:])], y=shares, marker_color="#B45309"))
fig.update_layout(title="Share of clay / white-mica pixels by distance to the nearest lineament", xaxis_title="km",
                  yaxis_title="% of identified pixels", height=400, **LAYOUT)
fig.show()

Clay / white-mica share of identified pixels: 33.6 % within 500 m of lineaments, 22.6 % beyond 1.5 km → enrichment ratio 1.49

Structures from elevation data — lineaments, faults and alteration¶

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

Lighting the landscape from several directions¶

Extracting lineaments automatically¶

Testing a hypothesis: is alteration found near structures?¶

Reading the result critically¶