Data and code for "Global Analysis of Hack's Law" (version 1.1)

This package holds the basin tables, the programs that compute them from the input products, and the analysis and figure scripts of Global Analysis of Hack's Law: Decomposing the Geometric Origins of h > 0.5 Across Scales and Sampling Designs (Water Resources Research, manuscript 2026WR044573, revised version of October 2026). Version 1.1 accompanies the revised manuscript. It replaces version 1.0 (https://doi.org/10.4211/hs.d29b097e15524e13afc61cd229819b96), which the submitted manuscript cites; CHANGES.md lists the differences. Text S1–S2, Tables S1–S2 and Figures S1–S4 named in this package are those of the material for the reviewers of the manuscript.

Files

Unpack the four zip files into one folder:

SourceCode/Python/ the programs that compute the tables from the input products (SourceCode/Python/README.md) Analysis/scripts/ the analysis scripts; run_analyses.sh runs them Analysis/out/ where the analysis scripts write (empty in the package) Analysis/audit/ the record of the check of every along-stem sample against the flow directions Materials/Scripts/ the figure scripts and their shared helpers Materials/Data/ the 90 m tables (Materials/Data/README_tables.txt); results/ holds the analysis outputs that the figure scripts read Materials/Data30m/ the 1 arc-second tables of South and North America Materials/Data30m_match/ for each 1 arc-second basin, the 90 m basin it overlaps most, with the intersection over union Materials/DataThreshold/ the channel-head experiment Materials/DataLakeLines/ the junction samples excluded as flow lines on lake surfaces Materials/Data_version1/ the basin attributes of version 1.0, read where the two versions are compared

The tables

The inventory holds every basin of at least 1 km² on MERIT Hydro (3 arc-seconds, about 90 m): 340,991 basins, each the drainage area of one terminal pixel of the flow-direction field, with its outlet type. The analyses use the 331,299 basins whose outlet lies on the coast (outlet_dir == 0 in basin_attributes_fine_only.csv; outlet_type == "coastal" in the 1 arc-second tables).

The symbols of the paper are columns of the tables: the length of the longest flow path L is lfp_length_km, the straight source-to-outlet distance D is source_outlet_dist_km, the tortuosity Ω = L/D is lfp_sinuosity, the basin length L_b is basin_axial_length_km (the 90 m tables hold the same values once more as basin_length_km) and the elongation ratio R_e is elongation_ratio. basin_transverse_length_km is the width of the basin across L_b; the paper does not analyze it, and Figure 10 uses it only to keep the basins whose width is positive and finite. The perimeter P (basin_perimeter_km) is the length of a Moore-neighborhood walk around the boundary pixels of the basin, with the weights of Vossepoel and Smeulders (1982, https://doi.org/10.1016/0146-664X(82)90057-0): 0.980 of the pixel edge for a step to an edge neighbor, 1.406/√2 of the pixel diagonal for a step to a corner neighbor, and 0.091 of the mean pixel edge taken off wherever the walk changes from one kind of step to the other, each step with the edge lengths of the pixels at its latitude. The shape indices, among them the Gravelius compactness G_c = P/(2√(πA)) (gravelius_compactness), take A as basin_area_km2, the upstream area at the outlet pixel. Distances and lengths are measured on the WGS84 ellipsoid.

Materials/Data/README_tables.txt describes every table. The tables of Materials/Data/ are put together from the global tables that the programs write, without changing any value: basin_attributes_fine_only.csv joins, basin by basin, the tables of the flow path, the planform and the terrain (relief and hypsometry); basin_attributes_120th_only.csv those of aridity, fault distance, Köppen–Geiger and rainfall erosivity; basin_attributes_400th_only.csv those of landform and lithology and of soil; basin_attributes_020th_only.csv is the table of peak ground acceleration. The tables of a group hold the same basins in the same order under the same five key columns (region_id, global_basin_id, local_basin_id, outlet_dir, basin_area_km2); a column that two of them write is kept once; the flow-path table's lfp_length_vs_km and lfp_sinuosity_vs are left out; the planform table's basin_length_km is named basin_axial_length_km; and ai_frac_arid_incl_hyper is added as ai_frac_hyper_arid + ai_frac_arid (−9999 where either is missing). basin_perbasin_hackfit_fine.csv, alongstem_global.csv and tributaryjunction_global.csv are the global tables as the programs write them; partition_regions.csv and basin_blocks.csv come from the region table of the partition (hacklaw_blocks.py describes the blocks); the three tables of Materials/DataThreshold/ are the 90 region tables of the channel-head program one after the other, in numeric region order.

Figures

The figure scripts are run from Materials/Scripts/ and write to Materials/Figures/:

cd Materials/Scripts && python3 Fig08_Geometric_Decomposition.py

On the authors' machine (macOS, matplotlib 3.10.8) the scripts draw the figures of the paper pixel for pixel: Figure 7 with matplotlib's default backend on macOS, the others with the Agg backend (MPLBACKEND=Agg); with other fonts or another backend the margins of a figure can differ by a few pixels. The literature values of Figure 1 are written in its script, each with its source. Figure 7 computes the 95% bands of its sliding windows on its first run (a block bootstrap, 2,000 draws) and keeps them next to the figure in Fig07_window_intervals.csv (delete it to recompute the bands after changing the tables); it also prints the window maxima quoted in Section 4.1. Fig07_30m_Exponent_vs_Area.py and Fig08_30m_Geometric_Decomposition.py are the 1 arc-second counterparts of Figures 7 and 8; FigS_30m_Supplement.py draws Figures S1 and S2 with the panels of the second, and paper_numbers.py uses the windows of the first. Data_Appendix.py made the example basins of Figures A1, A3 and A5 (Materials/Data/cases/) from the coarse basin raster in Materials/Data/Global/ and, for the reaches of Figure A5, from per-region Strahler reach files of the 90 m chain that are not part of the package. rf_robustness_checks.py runs the robustness checks of the random forest (gradient boosting on the same folds) and writes Materials/Data/rf_robustness_results.json. The shared helpers are coastal_filter.py (coastal basins only), figure_helpers.py (log axes and their tick labels), hacklaw_blocks.py (the blocks of the bootstrap and of the cross-validation) and hacklaw_lake_lines.py (drops the junction samples on lake surfaces).

Analyses and numbers

bash Analysis/scripts/run_analyses.sh lane_a # intervals, decomposition, channel-head fits, 1 arc-second comparison bash Analysis/scripts/run_analyses.sh lane_b # random forests (several hours) bash Analysis/scripts/run_analyses.sh lane_c # block bootstrap of the random-forest importances (500 draws) bash Analysis/scripts/run_analyses.sh numbers # the numbers of the paper not read off a table (paper_numbers.py)

The runner can be started from any directory; set PYTHON to choose the interpreter (default python3). It runs every step on one thread, so that the random forests give the same result on every run, records what it ran in Analysis/out/progress_<lane>.txt, and stops if the along-stem samples do not match the record in Analysis/audit/. The numbers lane reads the outputs of the other three. The analysis outputs that the figure scripts read are included in Materials/Data/results/ and Materials/DataThreshold/fits/, so the figures can be drawn without running the analyses. To draw the figures from a new run, copy its outputs from Analysis/out/ into Materials/Data/results/, and those of Analysis/out/threshold/ into Materials/DataThreshold/fits/; then draw Figure 11 before Figure 12, because Figure 11 rewrites Materials/Data/rf_importance_results.pkl, from which Figure 12 takes its predictors.

Further steps of lane_a: threshold_basin_length.py (the exponent of L_b and the trend of R_e at every channel-head threshold, and the bounds on the independent-outlet exponent, Section 5.1 and Text S2), tortuosity_by_product_class.py (the median tortuosity of the two products by area class, and the share of the coastal basins above 100 km² by count and by area), paired_product_comparison.py (the two products on the same basins, with paired block-bootstrap intervals of the changes of h, h_D and β, Text S1), shape_term_checks.py (the shape term by area class at 90 m and at 1 arc-second, and the change of h_D and β between the tables of versions 1.0 and 1.1) and singer_convention_and_blocks.py (the fits in the convention of Singer et al., 2024, of Section 5.1: L in m and A in m² by nonlinear least squares over the area range of GDBM; and the intervals of Table 3 and of the decomposition under other blocks and other random seeds).

Scripts run on their own, from the package root (python3 Analysis/scripts/<name>.py); they write to Analysis/out/ or print:

coastal_alongstem.py loads the coastal basins and their along-stem samples, with D measured from each basin's source pixel, for alongpath_decomposition_l2.py and the two same_basins_*.py scripts.

validate_dir_sampling.py made the record in Analysis/audit/: it follows the flow directions from each basin's source pixel to the outlet and checks every along-stem sample and every per-basin fit against them. It reads the region tiles of the partition (flow direction, basin and upstream area, and the flow distances ldn and lup as Float32 rasters), which are not part of the package. The record is the summary of the check; the per-region replay files and the inventory extract it names (validation_details/) are not included.

The junction samples on lake surfaces are dropped inside every script that fits the pooled junction samples (they are all of orders 1 and 2, which the per-basin junction fits leave out); a fit made directly on Materials/Data/tributaryjunction_global.csv without hacklaw_lake_lines.py includes them.

Computing the tables

SourceCode/Python/README.md describes what each program computes, the inputs it reads and how to run it. In short:

cd SourceCode/Python/merit90m && python run_merit90m.py --regions all --programs all --merge
cd ../channel_head_threshold && python run_threshold.py --regions all --merge
cd ../hydrosheds30m && python run_hydrosheds30m.py --continents south-america,north-america \
    --programs lfp,alongstem,tributaryjunction,perbasin_hackfit,planform,basin_match,attributes

The programs read the products listed in DATA_SOURCES.md from a data root (--data-root or HACKLAW_DATA_ROOT). --out-root is the output folder of one experiment; without it, the programs write to merit90m/ or hydrosheds30m/ under HACKLAW_OUT_ROOT (default <data root>/Data_HackLaw/tables/). The 1 arc-second programs write their tables to <output folder>/<continent>/table/, which the package holds as Materials/Data30m/<continent>/global/table/, with the attribute tables beside them in Materials/Data30m/ and the basin matches in Materials/Data30m_match/<continent>/. --programs all also runs six attribute steps that the paper does not use (qacc_era5, qacc_grades, koppen, landcover, rgi2000, strainrate), most of whose inputs DATA_SOURCES.md does not describe; to leave them out, use --programs shape,aridity,soil,landform_lithology,faultdist,pga475,rfactor,hack (the merge skips the tables that were not made). On another machine the last digit of a few quantities can differ, because sine, cosine, logarithm and powers of ten come from that machine's math library.

Requirements

Python 3.10 or later. Analysis and figures: numpy, pandas, scipy, statsmodels, scikit-learn (1.3 or later), matplotlib, geopandas, shapely. Programs: numpy, numba, pandas, rasterio and GDAL's Python bindings; SourceCode/Python/tools/build_elv_raw_global.sh runs under zsh. The versions used are listed in DATA_SOURCES.md. The package was run on macOS. On Windows use WSL: the runner is a bash script, the scripts print UTF-8 text, and the programs use fork-based worker processes and hard links.

License

Code: MIT (LICENSE). Tables: CC BY 4.0. The input products keep their own licenses (DATA_SOURCES.md).

Use of AI tools

Section 3.6 of the paper states how AI tools were used in the study. Claude Code (Anthropic) and Codex (OpenAI) were also used to help organize this package (file names, comments and documentation) and to check that the scripts reproduce the figures and numbers of the paper; the authors reviewed every change.

Citation

Jiang, L., Wu, H., Chen, W., Qu, L., Yang, T., & Xia, Y. (2026). Data and code for "Global Analysis of Hack's Law: Decomposing the Geometric Origins of h > 0.5 Across Scales and Sampling Designs" (Version 1.1) [Dataset and Software]. HydroShare. https://doi.org/10.4211/hs.ce5b73513a474abcb57c11bd1a7f82c4