#!/usr/bin/env python3 # -*- coding: utf-8 -*- """Figure A3: dfs_in and dfs_out from the released seq_dfs and the upg of MERIT Hydro, on the 8 x 8 example of Figure 4. No input data. Output: ../Figures/FigA4_DFS_Conversion_v11.png and .pdf. Run from this folder; ../Data holds the inputs and the figure is written to ../Figures. """ import os import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt mpl.rcParams.update({ "font.family": "Helvetica", "font.size": 9, "axes.titlesize": 10, "mathtext.fontset": "cm", "axes.linewidth": 0.6, }) HERE = os.path.dirname(os.path.abspath(__file__)) OUT_DIR = os.path.join(HERE, "..", "Figures") # ---- the 8 x 8 example of Fig. 4 ------------------------------------------------------------------------------------- NR, NC = 8, 8 SIZE = NR * NC DIR = np.array([ [ 1, 2, 4, 8, 8, 8, 16, 16], [128, 2, 4, 8, 8, 16, 32, 32], [ 2, 2, 4, 8, 8, 8, 8, 8], [ 4, 2, 4, 8, 8, 8, 8, 8], [ 2, 4, 8, 8, 16, 16, 8, 8], [ 2, 4, 8, 8, 4, 8, 8, 8], [ 0, 0, 0, 8, 2, 4, 8, 8], [247, 247, 0, 0, 0, 0, 0, 0]]) MSK = np.array([ [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 0, 1, 1, 1, 1], [0, 1, 0, 0, 1, 1, 1, 0], [0, 0, 0, 0, 0, 1, 0, 0]]) BSN = np.array([ [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 2], [1, 1, 1, 1, 1, 1, 2, 2], [1, 1, 1, 0, 2, 2, 2, 2], [0, 1, 0, 0, 2, 2, 2, 0], [0, 0, 0, 0, 0, 2, 0, 0]]) D8 = {1: (0, 1), 2: (1, 1), 4: (1, 0), 8: (1, -1), 16: (0, -1), 32: (-1, -1), 64: (-1, 0), 128: (-1, 1)} BLUE, GREEN, RED = "#2563EB", "#16A34A", "#E53935" SHADE = "#fff1b8" T1 = 26 def build(): down = np.full(SIZE, -1) for r in range(NR): for c in range(NC): if MSK[r, c]: d = int(DIR[r, c]) down[r * NC + c] = r * NC + c if d == 0 else (r + D8[d][0]) * NC + c + D8[d][1] children = [[] for _ in range(SIZE)] outlets = [] for i in range(SIZE): if down[i] == i: outlets.append(i) elif down[i] >= 0: children[down[i]].append(i) def upstream(i): out = [i] for k in children[i]: out += upstream(k) return out upg = np.array([len(upstream(i)) if down[i] >= 0 else 0 for i in range(SIZE)]) for i in range(SIZE): children[i].sort(key=lambda k: -upg[k]) din = np.full(SIZE, -1) t = 0 for o in outlets: # basins in the order of their outlet pixel, as FlowTopo stack = [o] while stack: n = stack.pop() din[n] = t t += 1 stack.extend(reversed(children[n])) seq = np.full(SIZE, -1) for o in outlets: # FlowTopo: the block of the basin, reversed members = upstream(o) a, b = din[members].min(), din[members].max() for i in members: seq[i] = a + b - din[i] # the conversion of the appendix, from seq_dfs and upg only din2 = np.full(SIZE, -1) blocks = {} for o in outlets: b = seq[o] a = b - upg[o] + 1 blocks[o] = (a, b) for i in upstream(o): assert a <= seq[i] <= b din2[i] = a + b - seq[i] assert np.array_equal(din2, din), "dfs_in from seq_dfs and upg differs from the depth-first numbers" dout = np.where(din >= 0, din + upg - 1, -1) for i in range(SIZE): if din[i] >= 0: up = upstream(i) assert sorted(din[up]) == list(range(din[i], din[i] + upg[i])), i assert dout[i] == din[up].max(), i assert (seq[T1], upg[T1], din[T1], dout[T1]) == (37, 20, 3, 22), (seq[T1], upg[T1], din[T1], dout[T1]) return outlets, blocks, upg, seq, din, dout, set(upstream(T1)) def colour(r, c): return BLUE if BSN[r, c] == 1 else GREEN def grid(ax, title, note): ax.set_xlim(-0.5, NC - 0.5) ax.set_ylim(NR - 0.5, -0.5) ax.set_aspect("equal") ax.set_xticks(np.arange(-0.5, NC, 1)) ax.set_yticks(np.arange(-0.5, NR, 1)) ax.set_xticklabels([]) ax.set_yticklabels([]) ax.grid(True, color="0.4", lw=0.35, alpha=0.45, zorder=1) ax.tick_params(length=0) ax.set_title(title, loc="left", fontweight="bold", pad=3) ax.text(0.5, -0.04, note, fontsize=7.2, va="top", ha="center", transform=ax.transAxes, linespacing=1.35, bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="0.6", lw=0.4), clip_on=False) def panel(ax, values, block, outlets, title, note, framed=False): for i in range(SIZE): r, c = divmod(i, NC) if not MSK[r, c]: continue if i in block: ax.add_patch(plt.Rectangle((c - 0.5, r - 0.5), 1, 1, facecolor=SHADE, edgecolor="none", zorder=0.5)) kw = dict(ha="center", va="center", fontsize=8.5, fontweight="bold", color=colour(r, c), zorder=5) if framed: kw.update(fontsize=7.5, bbox=dict(boxstyle="round,pad=0.04", fc="white", ec="black", lw=0.35)) ax.text(c, r, str(int(values[i])), **kw) if i in outlets: ax.plot(c + 0.33, r + 0.30, marker="v", ms=4.2, mfc=RED, mec="black", mew=0.4, zorder=6) r, c = divmod(T1, NC) ax.add_patch(plt.Rectangle((c - 0.46, r - 0.46), 0.92, 0.92, fill=False, edgecolor="black", lw=0.9, zorder=6)) ax.text(c - 0.5, r - 0.5, "T$_1$", ha="center", va="center", fontsize=5.5, zorder=7, bbox=dict(boxstyle="round,pad=0.12", fc="white", ec="black", lw=0.5)) grid(ax, title, note) def main(): outlets, blocks, upg, seq, din, dout, block = build() (oa, (aa, ba)), (ob, (ab, bb)) = sorted(blocks.items(), key=lambda kv: kv[1][0]) fig, axes = plt.subplots(2, 2, figsize=(6.4, 7.6)) panel(axes[0, 0], seq, block, outlets, r"(a) $\mathit{seq\_dfs}$ (FlowTopo)", "numbered as $\\mathit{dfs\\_in}$, reversed within each basin:\n" f"the outlet holds the basin's largest number $b$; T$_1$: {seq[T1]}") panel(axes[0, 1], upg, block, outlets, r"(b) $\mathit{upg}$ (MERIT Hydro)", f"upstream pixel count, the pixel included;\nT$_1$: {upg[T1]}") panel(axes[1, 0], din, block, outlets, r"(c) $\mathit{dfs\_in}$", r"$\mathit{dfs\_in} = a + b - \mathit{seq\_dfs}$, $b = \mathit{seq\_dfs}$(outlet), $a = b - \mathit{upg}$(outlet) + 1" f"\nbasin A [{aa}, {ba}], basin B [{ab}, {bb}]; T$_1$: {aa} + {ba} $-$ {seq[T1]} = {din[T1]}") panel(axes[1, 1], dout, block, outlets, r"(d) $\mathit{dfs\_out}$", r"$\mathit{dfs\_out} = \mathit{dfs\_in} + \mathit{upg} - 1$, the largest $\mathit{dfs\_in}$ upstream" f"\nT$_1$: {din[T1]} + {upg[T1]} $-$ 1 = {dout[T1]}", framed=True) fig.subplots_adjust(left=0.02, right=0.98, top=0.96, bottom=0.08, wspace=0.10, hspace=0.32) for ext in ("png", "pdf"): out = os.path.join(OUT_DIR, f"FigA4_DFS_Conversion_v11.{ext}") if os.path.exists(out): raise SystemExit(f"{out} is there already") fig.savefig(out, dpi=300, bbox_inches="tight", facecolor="white") print(f"Saved: {out}") print(f" basin A block [{aa}, {ba}] (outlet i{oa}), basin B [{ab}, {bb}] (outlet i{ob}); T1 = i{T1}: seq_dfs " f"{seq[T1]}, upg {upg[T1]}, dfs_in {din[T1]}, dfs_out {dout[T1]}; seq_dfs of T1's upstream " f"{min(seq[list(block)])}..{max(seq[list(block)])}") plt.close(fig) if __name__ == "__main__": main()