#!/usr/bin/env python # -*- coding: utf-8 -*- """Figure 3: (a) the 90 processing regions; (b) example rows of the Regional Attribute Table. Input: ../Data/rgn_global_200th_v11.gpkg, region_fine_global_v11.csv. Output: ../Figures/Fig03_GlobalRGN_RAT_v11.png and .pdf. Run from this folder; ../Data holds the inputs and the figure is written to ../Figures. """ from pathlib import Path import time import colorsys import numpy as np import pandas as pd import rasterio import geopandas as gpd import matplotlib as mpl import matplotlib.pyplot as plt import matplotlib.colors as mcolors from matplotlib.colors import ListedColormap from matplotlib.patches import Rectangle # ============================================================================ # Paths # ============================================================================ HERE = Path(__file__).resolve().parent VECTOR = HERE.parent / "Data/rgn_global_200th_v11.gpkg" CSV = HERE.parent / "Data/region_fine_global_v11.csv" OUT_DIR = HERE.parent / "Figures" # ============================================================================ # Region grid — CSV rows are at 3-arcsec resolution (1°=1200 cells) # ============================================================================ CELLS_PER_DEG = 1200 W_FULL = 360 * CELLS_PER_DEG # 432000 LON_MIN = -180.0 LAT_MIN = -60.0 LAT_MAX = 85.0 def col_to_lon(c: float) -> float: return LON_MIN + c / CELLS_PER_DEG def row_to_lat(r: float) -> float: return LAT_MAX - r / CELLS_PER_DEG # ============================================================================ # Golden-ratio HSV palette (mirrors Figure_Region.ipynb) # ============================================================================ def distinct_palette(n: int, s_cycle=(0.72, 0.62, 0.82), v_cycle=(0.92, 0.82, 0.72)): phi = (1 + 5 ** 0.5) / 2.0 g = 1.0 / phi out = [] for i in range(n): h = (i * g) % 1.0 s = s_cycle[i % len(s_cycle)] v = v_cycle[i % len(v_cycle)] r, gg, b = colorsys.hsv_to_rgb(h, s, v) out.append((r, gg, b)) return out # ============================================================================ # Panel (a) # ============================================================================ def draw_panel_a(ax, gdf_boundary: gpd.GeoDataFrame, df: pd.DataFrame): # --- raster basemap, coloured per region_id via palette --- t0 = time.time() uniq_ids = np.sort(gdf_boundary["region_id"].astype(int).unique()).astype(np.int64) palette = distinct_palette(max(len(uniq_ids), 100)) cmap = ListedColormap(palette[: len(uniq_ids)]) id2idx = {int(rid): i for i, rid in enumerate(uniq_ids)} gdf_fill = gdf_boundary.copy() gdf_fill["_idx"] = gdf_fill["region_id"].astype(int).map(id2idx) gdf_fill.plot(ax=ax, column="_idx", cmap=cmap, edgecolor="none", linewidth=0, rasterized=True) print(f" filled {len(gdf_fill)} regions: {time.time()-t0:.1f}s") # --- region boundary (black, thin) --- t0 = time.time() gdf_boundary.boundary.plot(ax=ax, linewidth=0.4, color="k", zorder=4, rasterized=True) # map as image in the PDF print(f" boundary overlay: {time.time()-t0:.1f}s") # --- region rectangles + labels from CSV --- def draw_rect(lonL, lonR, latB, latT, edgecolor, ls="--", lw=1.0, z=5): if lonR <= lonL or latT <= latB: return ax.add_patch(Rectangle( (lonL, latB), lonR - lonL, latT - latB, fill=False, ec=edgecolor, lw=lw, linestyle=ls, zorder=z, )) def draw_label(lonL, lonR, latB, latT, text, fs=14, color="black", box=False, where="center"): cx = (lonL + lonR) / 2.0 if where == "above": # just above the rectangle cy, va = latT + 0.6, "bottom" elif where == "below": # just below the rectangle cy, va = latB - 0.6, "top" else: cy, va = (latB + latT) / 2.0, "center" kw = {} if box: # keep the dashed outlines from crossing the digits kw["bbox"] = dict(boxstyle="round,pad=0.12", fc="white", ec="none", alpha=0.85) ax.text(cx, cy, text, ha="center", va=va, fontsize=fs, color=color, weight="bold", zorder=6, **kw) BLUE_DARK = "#1f4ed8" # antimerid label / rect GREEN_DARK = "#0a7a2a" # island label / rect SMALL_DEG = 2.0 # island windows this small (either side) get the label above them NEAR_LON, NEAR_LAT = 6.0, 4.0 # an island label this close to a continental label moves below its window def halves(R): """Window of a region in degrees, split at the dateline if it crosses it: list of (lonL, lonR, latB, latT); the first entry carries the label.""" rmin = int(R["region_row_min"]); rmax = int(R["region_row_max"]) cmin = int(R["region_col_min"]); cmax = int(R["region_col_max"]) latT, latB = row_to_lat(rmin), row_to_lat(rmax) if cmin >= 0 and cmax <= W_FULL: return [(col_to_lon(cmin), col_to_lon(cmax), latB, latT)] if cmin < 0: # negative start cmin_n, cmax_p = cmin, cmax else: # v1.1 form: end beyond 432000 cmin_n, cmax_p = cmin - W_FULL, cmax - W_FULL # right (Asian) half first: the label sits there only return [(col_to_lon(W_FULL + cmin_n), col_to_lon(W_FULL), latB, latT), (col_to_lon(0), col_to_lon(cmax_p), latB, latT)] # Continental (region_id < 10000): no rectangle, BLACK 14pt label # Antimeridian (10000–19999): BLUE dashed rectangle + BLUE 11pt label # Islands (>= 20000): GREEN dashed rectangle + GREEN 10pt label, on a white box # A region crossing the dateline is split there (v1.1 also has two # continental ones, 3513 and 55); its label sits on the right half. cont_xy = [] for _, R in df.iterrows(): rid = int(R["region_id"]) if rid < 10000: lonL, lonR, latB, latT = halves(R)[0] draw_label(lonL, lonR, latB, latT, str(rid), fs=14, color="black") cont_xy.append(((lonL + lonR) / 2.0, (latB + latT) / 2.0)) for _, R in df.iterrows(): rid = int(R["region_id"]) if rid < 10000: continue parts = halves(R) if rid < 20000: ec, fs, col, lw = BLUE_DARK, 11, BLUE_DARK, 1.4 else: ec, fs, col, lw = GREEN_DARK, 10, GREEN_DARK, 1.2 for lonL, lonR, latB, latT in parts: draw_rect(lonL, lonR, latB, latT, edgecolor=ec, lw=lw) lonL, lonR, latB, latT = parts[0] cx, cy = (lonL + lonR) / 2.0, (latB + latT) / 2.0 where = "center" if min(lonR - lonL, latT - latB) <= SMALL_DEG: where = "above" if any(abs(cx - x) < NEAR_LON and abs(cy - y) < NEAR_LAT for x, y in cont_xy): where = "below" draw_label(lonL, lonR, latB, latT, str(rid), fs=fs, color=col, box=True, where=where) # --- axis style (matches Fig 2 (a) and the user's ipynb) --- ax.set_xlabel("Longitude", fontsize=20, fontweight="bold") ax.set_ylabel("Latitude", fontsize=20, fontweight="bold") ax.set_xlim(-180, 180) ax.set_ylim(-60, 85) ax.set_xticks(np.arange(-180, 181, 60)) ax.set_yticks(np.arange(-60, 86, 30)) ax.set_xticks(np.arange(-180, 181, 30), minor=True) ax.set_yticks(np.arange(-60, 86, 15), minor=True) ax.grid(which="major", linestyle="--", linewidth=0.4, alpha=0.6) ax.grid(which="minor", linestyle=":", linewidth=0.25, alpha=0.4) ax.tick_params(axis="both", which="major", labelsize=18, length=6, width=1) ax.tick_params(axis="both", which="minor", labelsize=18, length=3, width=0.8) # ============================================================================ # Panel (b) — Regional Attribute Table # ============================================================================ ELL = "…" PICK_IDS = [11, 12, 43, 5618, 10001, 10002, 20101, 20901] DISPLAY_COLS = [ "region_id", "_ell1", "region_row_min", "region_row_max", "region_col_min", "region_col_max", "_ell2", "region_nrow", "region_ncol", "_ell3", "region_grid_count", "_ell4", ] # Relative column widths (sum = 1). Wider for region_cell_count which holds # 9-digit numbers (e.g., 381,369,152); ellipsis columns are slim. COL_WIDTHS = [ 0.09, # region_id up to "20,029" 0.04, # _ell1 0.10, # row_min up to "146,400" 0.10, # row_max up to "166,800" 0.09, # col_min includes -6,000 / -1,200 0.10, # col_max up to "363,600" 0.04, # _ell2 0.09, # nrow up to "60,000" 0.09, # ncol up to "73,200" 0.04, # _ell3 0.18, # region_grid_count up to "991,364,889" (region 5618) 0.04, # _ell4 ] assert abs(sum(COL_WIDTHS) - 1.0) < 1e-9 assert len(COL_WIDTHS) == len(DISPLAY_COLS) def fmt_int(v) -> str: if isinstance(v, str): return v return f"{int(round(v)):,}" def wrap_at_first_underscore(name: str) -> str: """'region_cell_count' → 'region\\ncell_count'; '_ellN' → ellipsis.""" if name.startswith("_ell"): return ELL if "_" not in name: return name a, b = name.split("_", 1) return f"{a}\n{b}" def fetch_rat_picks(df_full: pd.DataFrame) -> pd.DataFrame: sub = df_full[df_full["region_id"].isin(PICK_IDS)].copy() # Order by PICK_IDS list (manual ordering to match display sequence) sub["_ord"] = sub["region_id"].map({rid: i for i, rid in enumerate(PICK_IDS)}) return sub.sort_values("_ord").drop(columns="_ord").reset_index(drop=True) def draw_panel_b(ax, df_b: pd.DataFrame, all_ids: np.ndarray): """Insert an ellipsis row between picks ONLY when the source CSV actually contains region_ids strictly between them, and never as a trailing row if no rows follow the last pick in the CSV.""" all_ids = np.sort(np.asarray(all_ids)) body = [] for k in range(len(df_b)): cur_id = int(df_b.iloc[k]["region_id"]) if k > 0: prev_id = int(df_b.iloc[k-1]["region_id"]) has_between = bool(((all_ids > prev_id) & (all_ids < cur_id)).any()) if has_between: body.append({c: ELL for c in DISPLAY_COLS}) rec = {} for c in DISPLAY_COLS: rec[c] = ELL if c.startswith("_ell") else fmt_int(df_b.iloc[k][c]) body.append(rec) # Trailing ellipsis only if any real id follows the last pick. last_id = int(df_b.iloc[-1]["region_id"]) if (all_ids > last_id).any(): body.append({c: ELL for c in DISPLAY_COLS}) cell_text = [[r[c] for c in DISPLAY_COLS] for r in body] col_labels = [wrap_at_first_underscore(c) for c in DISPLAY_COLS] ax.axis("off") table = ax.table( cellText=cell_text, colLabels=col_labels, colWidths=COL_WIDTHS, bbox=[0.0, 0.0, 1.0, 1.0], cellLoc="center", colLoc="center", ) table.auto_set_font_size(False) # Per-cell styling — mirrors Fig 2 v4 (b) for (r, c), cell in table._cells.items(): cell.set_edgecolor("black") cell.set_linewidth(0.6) if r == 0: cell.set_facecolor("#1F4E79") cell.get_text().set_color("white") cell.get_text().set_weight("bold") cell.get_text().set_fontsize(16) cell._text.set_multialignment("center") else: cell.get_text().set_fontsize(16) cell.set_facecolor("#f7f7f7" if r % 2 == 0 else "white") # Double-height header row. n_rows = 1 + len(cell_text) body_h = 1.0 / (n_rows + 1) head_h = 2.0 * body_h for (r, c), cell in table._cells.items(): cell.set_height(head_h if r == 0 else body_h) # ============================================================================ # Main # ============================================================================ def main(): print("Loading inputs …") t0 = time.time() gdf = gpd.read_file(VECTOR) print(f" vector: {len(gdf):,} regions, {time.time()-t0:.1f}s") df_full = pd.read_csv(CSV, sep=r"\s+") df_full = df_full.rename(columns={ "row_min": "region_row_min", "row_max": "region_row_max", "col_min": "region_col_min", "col_max": "region_col_max", "nrow": "region_nrow", "ncol": "region_ncol", }) df_b = fetch_rat_picks(df_full) print(f" csv: {len(df_full)} rows; picks = {df_b['region_id'].tolist()}") # --- figure geometry (mirrors Fig 2 v4) --- fig_w_in = 18.0 bot_pad_in = 0.05 b_h_in = 4.60 gap_in = 1.00 a_h_in = 5.80 top_pad_in = 0.55 fig_h_in = bot_pad_in + b_h_in + gap_in + a_h_in + top_pad_in # 12.0 b_bot_in = bot_pad_in b_top_in = b_bot_in + b_h_in a_bot_in = b_top_in + gap_in a_top_in = a_bot_in + a_h_in fig = plt.figure(figsize=(fig_w_in, fig_h_in)) a_left, a_w = 0.06, 0.88 ax_a = fig.add_axes([ a_left, a_bot_in / fig_h_in, a_w, (a_top_in - a_bot_in) / fig_h_in, ]) b_left, b_w = 0.02, 0.96 ax_b = fig.add_axes([ b_left, b_bot_in / fig_h_in, b_w, (b_top_in - b_bot_in) / fig_h_in, ]) print("\nDraw panel (a) …") draw_panel_a(ax_a, gdf, df_full) print("\nDraw panel (b) …") draw_panel_b(ax_b, df_b, df_full["region_id"].values) # --- panel labels at fixed figure positions --- fig.text(a_left, (a_top_in + 0.10) / fig_h_in, "(a)", fontsize=30, fontweight="bold", va="bottom") b_label_y_in = b_top_in + gap_in * 0.30 fig.text(b_left, b_label_y_in / fig_h_in, "(b)", fontsize=30, fontweight="bold", va="bottom") OUT_DIR.mkdir(parents=True, exist_ok=True) for ext in ("png", "pdf"): p = OUT_DIR / f"Fig03_GlobalRGN_RAT_v11.{ext}" fig.savefig(p, dpi=300, bbox_inches="tight") print(f"\n wrote {p}") print(f" fig {fig_w_in:.1f} × {fig_h_in:.2f} in") plt.close(fig) if __name__ == "__main__": main()