The stages of project complexity

Why any of this exists

The short answer: research code fails at the seams between runs, not inside them. A script that computes the right number once is easy. What is hard is everything that happens afterwards — running it again with one input changed, running it four hundred times over scenarios, running it on someone else’s machine, running it on data too large to hold in memory, and still being able to say two years later exactly which code and which inputs produced figure 3.

Almost every research codebase answers those by copying: copy the script and edit one line, copy the block and change the index, comment out the slow part with a boolean. That works until it doesn’t, and when it stops working it does so silently — a stale intermediate file, a copied block where one index was never updated, a hardcoded path that meant something different on the machine the results came from.

The devstack’s answer is a single move applied over and over: make the thing that varies data rather than code. Scenarios become rows in a CSV instead of copied blocks. Machine-specific values become a parameters CSV instead of hardcoded paths. What has already been computed becomes files on disk instead of a commented-out boolean, so a re-run skips it automatically. What runs, and in what order, becomes a task tree instead of top-to-bottom statements — which is also what lets independent branches run in parallel. Hazelbean, and ProjectFlow inside it, are the machinery that makes that move cheap enough to be the default.

None of this is adopted up front. Each piece exists because a specific pain demanded it, and the rest of this page is that history: the questions that got progressively harder, and the tool each one forced. If your problem is still a one-script problem, the honest advice is to write the one script — the point of the stages is to let you recognize the moment you have outgrown it.

The narrative

As an introduction and/or motivation for using the software in the EE Devstack, I would like to talk through the process that I went through as a PhD student, postdoc and staff researcher to answer progressively harder questions. I break these out into 6 stages below. Solving these challenges is what led to the creation of Hazelbean and many other software solutions.

Stage 1, the old way: Simple question answered well

Here is an example script that you might write as an Earth-economy researcher. Suppose your adviser asks you “what is the total caloric yield on earth per hectare?” You might write a script like this:

import os
import numpy as np
import gdal

yield_per_hectare_raster_path = os.path.join('data', 'yield_per_cell.tif')
yield_per_hectare_raster = gdal.Open(yield_per_hectare_raster_path)
yield_per_hectare_array = yield_per_hectare_raster.ReadAsArray()

sum_of_yield = np.sum(yield_per_hectare_array)

print('The total caloric yield on earth per hectare is: ' + str(sum_of_yield))

Stage 2, the old way: Many similar questions. Creates a very long list.

This is where most research code goes to die, in my experience. Suppose your advisor now asks okay do this for the a bunch of different datasets on yield. The classic coder response is to make a longer script!

import os
import numpy as np
import gdal

yield_per_hectare_raster_path_1 = os.path.join('data', 'yield_per_cell_1.tif')
yield_per_hectare_raster_1 = gdal.Open(yield_per_hectare_raster_path_1)
yield_per_hectare_array_1 = yield_per_hectare_raster_1.ReadAsArray()

sum_of_yield_1 = np.sum(yield_per_hectare_array_1)

print('The total caloric yield on earth per hectare for dataset 1 is: ' + str(sum_of_yield_1))

yield_per_hectare_raster_path_2 = os.path.join('data', 'yield_per_cell_2.tif')
yield_per_hectare_raster_2 = gdal.Open(yield_per_hectare_raster_path_2)
yield_per_hectare_array_2 = yield_per_hectare_raster_2.ReadAsArray()

sum_of_yield_2 = np.sum(yield_per_hectare_array_2)

print('The total caloric yield on earth per hectare for dataset 2 is: ' + str(sum_of_yield_2))

yield_per_hectare_raster_path_3 = os.path.join('data', 'yield_per_cell_3.tif')
yield_per_hectare_raster_3 = gdal.Open(yield_per_hectare_raster_path_3)
yield_per_hectare_array_3 = yield_per_hectare_raster_3.ReadAsArray()

sum_of_yield_3 = np.sum(yield_per_hectare_array_2)

print('The total caloric yield on earth per hectare for dataset 3 is: ' + str(sum_of_yield_3))

yield_per_hectare_raster_path_4 = os.path.join('data', 'yield_per_cell_4.tif')
yield_per_hectare_raster_4 = gdal.Open(yield_per_hectare_raster_path_4)
yield_per_hectare_array_4 = yield_per_hectare_raster_4.ReadAsArray()

sum_of_yield_4 = np.sum(yield_per_hectare_array_4)

print('The total caloric yield on earth per hectare for dataset 4 is: ' + str(sum_of_yield_4))

This style of coding works, but will quickly cause you to lose your sanity. Who can find the reason the above code will cause your article to be retracted? Also, what if each of those summations takes a long time and you want to make a small change? You have to rerun the whole thing. This is bad.

Stage 3, the old way: Starting to deal with generalization, reusing code and shortening scripts.

The coding approach in stage 2 becomes intractable when there are lots of layers to consider. It’s also a pain to have to repeat code to do some common tasks, like loading the raster to a dataset and then to an array. This stage starts to apply good coding practices, such as defining helper functions like raster_to_array() below. Code is also made much shorter and more elegant by using loops. It minimizes the number of code statements, reduces bugs and scales better to long lists of input files.

import os
import numpy as np
import gdal

# NOTE 1: Helper function defined
def raster_to_array(raster_input_path):
    ds = gdal.Open(raster_input_path)
    print("Reading " + raster_input_path +'. This might take a while!')
    array = ds.ReadAsArray()
    return array

# NOTE 2: Inputs put into an iterable
input_paths = [
    'yield_per_cell_1.tif',
    'yield_per_cell_2.tif',
    'yield_per_cell_3.tif',
    'yield_per_cell_4.tif',
]

# NOTE 3: Calculation happens in loops, recording results to an output object
summations = []
for raster_path in input_paths:
    array = raster_to_array(raster_path)
    summations.append(np.sum(array))

print('Sums of layers: ' + str(summations))

Stage 4, the old way: Starting to deal with performance and generalization.

Below is a real-life script I created in around 2017 to calculate something for Johnson et al. 2016. Unlike the other levels, do not even attempt to run this, but just appreciate how awful it is. Please skim past it quickly to save me the personal embarassment! Instead, I provide a better example of code below that does things better, in the Earth-Economy Devstack way,

BAD CODE; DONT RUN.

import logging
import os
import csv
import math, time, random
from osgeo import gdal, gdalconst
import numpy as np

# NOTE 1: I started to pull in Cython (Python code compiled to C for speed) because my code was getting slow
import pyximport
pyximport.install(setup_args={"script_args":["--compiler=mingw32"],"include_dirs":numpy.get_include()}, reload_support=True)

# NOTE 2: I wrote my own Python Library (geoecon_utils), which went through several more 
# iterations (Numdal, Lol!), until it got finalized as hazelbean
import geoecon_utils.geoecon_utils as gu
import geoecon_utils.geoecon_cython_utils as gcu

# NOTE 3: Logging becomes important to manage information input-output used by the developer
log_id = gu.pretty_time()
LOGGER = logging.getLogger('ag_tradeoffs')
LOGGER.setLevel(logging.WARN) # warn includes the final output of the whole model, carbon saved.
file_handler = logging.FileHandler('logs/ag_tradeoffs_log_' + log_id + '.log')
LOGGER.addHandler(file_handler)

# NOTE 4: Defining inputs and outputs is now based on a workspace, from which everything
# else is defined with relative paths. Scales better with inputs and to other users.
workspace = 'E:/bulk_data/reusch_and_gibbs_carbon/data/datasets/c_1km/'
c_1km_file = 'c_1km.tif'
c_1km_uri = workspace + c_1km_file
ha_per_cell_5m_file = 'ha_per_cell_5m.tif'
ha_per_cell_5m_uri = workspace + ha_per_cell_5m_file

# NOTE 5: Here's an example of using custom libraries to 
ha_per_cell_5m = gu.as_array(ha_per_cell_5m_uri)

# NOTE 6: Here we start to deal with conditional running of code that skips outputs if they have already been
# created. This is often the first (and often most eficatious) optimization of code to run fast.
do_30s_resample = False
if do_30s_resample:
    # Define desired resample details. In this case, I am converting to 10x resolution of 5 min data (30 sec)
    desired_geotrans = (-180.0, 0.008333333333333, 0.0, 90.0, 0.0, -0.008333333333333)
    c_30s_unscaled_uri = workspace + 'c_30s_unscaled_' + gu.pretty_time() + '.tif'
    gu.aggregate_geotiff(c_1km_uri, c_30s_unscaled_uri, desired_geotrans)

# NOTE 7: An here, we see an incredibly slow approach that seems intuitive but is wrong
# because it is 1000x slower than correct vectorized calculations (which are provided by numpy)
array = raster_to_array(c_30s_unscaled_uri)
for row in range(array.shape[0]):
    for col in range(array.shape[1]):
        if array[row, col] < 0:
            array[row, col] = -9999 # Set negative values to the no-data-value

Good version; do run

# For reference, the correct way would have been as below. We will introduce hazelbean utils (like hb.as_array() soon.
array = hb.as_array(input_path)
array = np.where(array < 0, -9999, array)

Stage 5, the old way: Dealing with larger-than-memory data

Depending on the size of the array, the numpy where command used above would fail with a MemoryError or something similar. Below you will the first way that I dealt with this (BAD CODE) and then the correct way. In either case, it is almost always true that the most likely solution to larger-than-memory situations is to apply your algorithm in chunks. We’ll do this in both cases.

BAD CODE; DON’T RUN

In this code, I created a thing I made up, a “Tile Reference” which I implemented with geoecon_utils.Tr(). This returned a set of tiles, defined by their row, column, x_width and y_width. We used these to load just subsets of the array and do operations on the smaller thing.


literal_aggregation_to_5m_cell = False
if literal_aggregation_to_5m_cell:
    factor =  10
    shape = (2160, 4320)
    c_30s_tr = gu.Tr(c_30s_uri, chunkshape = shape)
    cell_sum_5m = np.zeros((c_30s_tr.num_rows / factor, c_30s_tr.num_cols / factor))
    cell_sum_5m_uri = os.path.split(ha_per_cell_5m_uri)[0] + '/c_5m_literal_aggregation_' + gu.pretty_time() + '.tif'
    for tile in c_30s_tr.tr_frame:
        tile_array = c_30s_tr.tile_to_array(tile)
        print 'Aggregating tile', tile
        for row in range(c_30s_tr.chunkshape[0] / factor):
            for col in range(c_30s_tr.chunkshape[1] / factor):
                cell_sum = np.sum(tile_array[row * factor : (row + 1) * factor, col * factor: (col + 1) * factor])
                cell_sum_5m[tile[0] / factor + row, tile[1] / factor + col] = cell_sum
    cell_sum_5m /= 100 #because 30s is in c per ha and i want c per 5min gridcell
    print gu.desc(cell_sum_5m)
    gu.save_array_as_geotiff(cell_sum_5m, cell_sum_5m_uri, ha_per_cell_5m_uri)

Better code, do run

The better way is to use a function that builds in the tiling functionality. Eventually this will expand to multi-computer approaches, but for now, we’ll just use the local but parallelized hb.raster_calculator()

def op(carbon_per_ha, ha_per_cell):
    return carbon_per_ha * ha_per_cell

hb.raster_calculator(
    [carbon_per_ha_path, ha_per_cell_path],
    op,
    carbon_per_cell_path)

Stage 5, redone with ProjectFlow, shows this call doing real work inside a run file.

Stage 6, the old way: Systematic file management.

The final level of complexity we will discuss (before just using the Earth Economy Devstack approach) arises when the number of files that must be managed becomes a challenge both for performance reasons and the challenges of managing complexity.

One example where this comes up is when a computation requires writing tiles of output. In many big-data applications and in most of the very large datasets that are available online, the data are themselves stored in tiles. On the one hand, this is nice because it automatically suggests a chunk-by-chunk parallelization strategy. On the other hand, it quickly becomes challenging when, for instance, you want to look at an area of interest (AOI) that spans multiple tiles. There are a plethora of software solutions to deal with this, such as GDAL’s Virtual Raster (VRT) file type, but many of these have limitations.

When the computation in question requires many complex steps which might be contingent on other intermediate products from adjacent tiles, even some of the most cutting-edge solutions that implement complex tiling architecture (like DASK with rioxarray) will not be sufficient. This was the challenge that arose when doing downscaling with the SEALS model, especially when the algorithm had to be trained on such tiles millions of times. The optimized algorithm in this complex data and computation dependency-tree situation required a new tool, which would also have to address all of the above challenges in project complexity.

This led to ProjectFlow, one of the key tools within the Earth Economy Devstack and a part of Hazelbean.

Three axes, not one ladder

The six stages above are a narrative: the pains, in the order I hit them, that produced the devstack. They are not the spec, and their numbers are not the spec’s numbers. When the spec says level, it means a position on one specific axis, defined here.

A run file’s complexity is three independent questions. A project sits somewhere on each, and moving along one axis does not require moving along another. Only the first is numbered, because only the first is a ladder climbed in order.

Configuration — how a run varies

level	what it adds	template
1	one task; constants inline	`run_template_1_minimal.py`
2	a task tree; constants inline	`run_template_2_script_with_tree.py`
3	+ a scenarios CSV: the rows of work become data	`run_template_3_canonical.py`
4	+ a parameters CSV: machine-specific values leave the code	`run_template_4_data_driven.py`

Code layout — where the code lives

position	what it means
single file	tasks defined in the run file
split	tasks in `<project>_tasks.py`; science helpers in `<project>_functions.py`; science-unaware helpers in `<project>_utils.py`; tree builders in `<project>_initialize_project.py`
library package	the above installed and imported by many projects

<project>_utils.py is not a filing category but a queue: science-unaware, and used by more than one project, means a function is a candidate for promotion into hazelbean. <project>_functions.py is terminal — it holds the science, and stays with the model. Task modules are plural and topic-named once one gets large (seals_process_coarse_timeseries.py, seals_visualization_tasks.py).

Ownership — who owns the code you run

position	what it means
self-contained	you own every file the run touches
devstack developer	you import a devstack library and could edit it
downstream user	the library is a read-only dependency in your own repo

The axes really are independent: global_invest’s service runners are configuration level 2 with a library-package layout. Templates 1–4 hold code layout and ownership fixed so that configuration is the only thing changing between them — that is what makes them readable as a diff. run_template_seals_example.py then moves the other two, staying at configuration level 3 while its layout becomes a library package.

Where to start: earning complexity rather than adopting it

Every template below is described on the Run Templates page, which is the place to start if you want to copy one.

Each tool in the six stages was adopted only when the previous stage’s pain demanded it — and that lesson applies to ProjectFlow itself. If your problem is still a stage-1 problem, do not start at configuration level 4. The spec scales down.

The smallest run file — a function becomes a task when you add it. The minimal ProjectFlow script is one task function plus a guarded three-line main block: p = hb.ProjectFlow(), p.add_task(your_function), p.execute(). Nothing is hidden — the add_task call is the moment a plain function becomes a task (and calling execute() with nothing added fails with a message telling you exactly that). See examples/run_templates/run_template_1_minimal.py — about 30 lines, and in exchange for roughly the effort of the stage-1 script you get an organized project directory, one folder per task, skip-if-already-computed re-runs, and logging.

The canonical form. When you want variants, graduate to configuration level 3: build_task_tree(p), run_project(p), and a __main__ guard that builds and configures the ProjectFlow. The rule that replaces a signature full of defaults: run_project(p) sets what no variant ever changes; the caller sets what a variant might. When a constant starts varying, it moves one line up, out of run_project and into the caller — no signature to edit and no default to keep in sync. See run_template_3_canonical.py, its annotated twin in examples/run_templates_annotated/ for the reasoning, and the conventions page for the full rules. For a project of your own that imports a devstack library as a read-only dependency, copy run_template_downstream_user.py into your repo as its run_<project>.py; for stock SEALS, copy run_seals_standard.py. The full-scale reference implementation is ngfs_pnas (a private repository — ask for access).

When to move along an axis. Configuration: add a scenarios CSV when you first have two variants of the same run (the CSV, not the code, is what varies); add a parameters.csv key the first time a machine-specific value — an ssh host, a credentials path — tries to enter your code; add a thin run_<project>_test.py wrapper when the full run gets long enough to want a pared smoke test. Code layout: split when a task module gets big enough to navigate, or when a second run file needs the same task. Ownership: you cross into downstream-user territory the moment the library you import is one you should not be editing.

Every one of these is the stage-2 lesson in a new costume: the moment you are tempted to copy a block — or a whole file — and edit one line, the spec has a mechanism that varies the data instead of forking the code.

Starting to walk the steps with ProjectFlow

Stage 1, redone with ProjectFlow

Old

import os
import numpy as np
import gdal

yield_per_hectare_raster_path = os.path.join('data', 'yield_per_cell.tif')
yield_per_hectare_raster = gdal.Open(yield_per_hectare_raster_path)
yield_per_hectare_array = yield_per_hectare_raster.ReadAsArray()

sum_of_yield = np.sum(yield_per_hectare_array)

print('The total caloric yield on earth per hectare is: ' + str(sum_of_yield))

New: P object organizes attributes

import os
import numpy as np
import gdal
import hazelbean as hb

# At this stage, p is honestly just a namespace — the code below works exactly as before.
# The payoff comes in the later stages, when everything already hanging off p gets you
# organized project directories, skip-if-already-computed re-runs, and auto-downloading
# data paths without restructuring your script. This line is the hook it all hangs on.
p = hb.ProjectFlow()

p.yield_per_hectare_raster_path = os.path.join('crops', 'yield_per_cell.tif') # Note for the future, we're defining a ref_path here but let's ignore that for now and assume it just goes in the p.input_dir
yield_per_hectare_raster = gdal.Open(p.yield_per_hectare_raster_path)
yield_per_hectare_array = yield_per_hectare_raster.ReadAsArray()

p.sum_of_yield = np.sum(yield_per_hectare_array)

print('The total caloric yield on earth per hectare is: ' + str(p.sum_of_yield))

One thing did already happen behind the scenes: hb.ProjectFlow() chose a sensible project directory for you. If your script lives inside a git repo (it probably does), ProjectFlow refuses to write outputs into the repo — it walks up to just outside it and uses <repo_parent>/projects/<name>, where <name> comes from your run file’s name. It logs the choice, and passing project_dir explicitly always overrides it. In the standard devstack layout this lands exactly where the full spec would put it, so nothing needs to change later as you climb the stages.

Stage 2, redone with ProjectFlow

Old

import os
import numpy as np
import gdal

yield_per_hectare_raster_path_1 = os.path.join('data', 'yield_per_cell_1.tif')
yield_per_hectare_raster_1 = gdal.Open(yield_per_hectare_raster_path_1)
yield_per_hectare_array_1 = yield_per_hectare_raster_1.ReadAsArray()

sum_of_yield_1 = np.sum(yield_per_hectare_array_1)

print('The total caloric yield on earth per hectare for dataset 1 is: ' + str(sum_of_yield_1))

yield_per_hectare_raster_path_2 = os.path.join('data', 'yield_per_cell_2.tif')
yield_per_hectare_raster_2 = gdal.Open(yield_per_hectare_raster_path_2)
yield_per_hectare_array_2 = yield_per_hectare_raster_2.ReadAsArray()

sum_of_yield_2 = np.sum(yield_per_hectare_array_2)

print('The total caloric yield on earth per hectare for dataset 2 is: ' + str(sum_of_yield_2))

yield_per_hectare_raster_path_3 = os.path.join('data', 'yield_per_cell_3.tif')
yield_per_hectare_raster_3 = gdal.Open(yield_per_hectare_raster_path_3)
yield_per_hectare_array_3 = yield_per_hectare_raster_3.ReadAsArray()

sum_of_yield_3 = np.sum(yield_per_hectare_array_2)

print('The total caloric yield on earth per hectare for dataset 3 is: ' + str(sum_of_yield_3))

yield_per_hectare_raster_path_4 = os.path.join('data', 'yield_per_cell_4.tif')
yield_per_hectare_raster_4 = gdal.Open(yield_per_hectare_raster_path_4)
yield_per_hectare_array_4 = yield_per_hectare_raster_4.ReadAsArray()

sum_of_yield_4 = np.sum(yield_per_hectare_array_4)

print('The total caloric yield on earth per hectare for dataset 4 is: ' + str(sum_of_yield_4))

New: Make it Iterate Over a list

import os
import numpy as np
import hazelbean as hb

p = hb.ProjectFlow()

# NOTE 1: No helper function needed — check hazelbean before writing one. The
# raster-to-array helper we wrote at stage 3 already exists as hb.as_array().
# (This is the devstack reuse rule.)

# NOTE 2: Inputs put into an iterable. The thing that varies is now DATA, not
# code — this list is the seed of what will later become a scenarios.csv.
p.input_paths = [
    os.path.join('data', 'yield_per_cell_1.tif'),
    os.path.join('data', 'yield_per_cell_2.tif'),
    os.path.join('data', 'yield_per_cell_3.tif'),
    os.path.join('data', 'yield_per_cell_4.tif'),
]

# NOTE 3: Calculation happens in a loop, recording results to an output object.
# This still re-reads every raster on every run — that pain is what the next
# stage's tasks (and their skip-if-already-computed re-runs) solve.
p.summations = []
for raster_path in p.input_paths:
    raster = hb.as_array(raster_path)
    sum_of_yield = np.sum(raster)
    p.summations.append(sum_of_yield)

print('Sums of layers: ' + str(p.summations))

Stage 3, redone with ProjectFlow

Old

import os
import numpy as np
import gdal

# NOTE 1: Helper function defined
def raster_to_array(raster_input_path):
    ds = gdal.Open(raster_input_path)
    print("Reading " + raster_input_path +'. This might take a while!')
    array = ds.ReadAsArray()
    return array

# NOTE 2: Inputs put into an iterable
input_paths = [
    'yield_per_cell_1.tif',
    'yield_per_cell_2.tif',
    'yield_per_cell_3.tif',
    'yield_per_cell_4.tif',
]

# NOTE 3: Calculation happens in loops, recording results to an output object
summations = []
for raster_path in input_paths:
    array = raster_to_array(raster_path)
    summations.append(np.sum(array))

print('Sums of layers: ' + str(summations))

New: The loop becomes a task

import os
import numpy as np
import hazelbean as hb

# NOTE 1: The loop moves inside a function, and the function becomes a TASK —
# explicitly, when we hand it to p.add_task() at the bottom. A task receives p
# and gets its own output directory, p.cur_dir, inside the project dir.
# Also note that we slightly deviate from python style: functions normally
# are verbs/command langugage (sum_yields). In EE spec, the name of a task
# becomes the folder that tasks' outputs are stored in, so we instead choose
# a noun.
def yield_summations(p):
    # Anything set BEFORE the p.run_this block is a project-level variable:
    # it is always assigned, even on runs where this task is skipped, so
    # later tasks can rely on it.
    p.input_paths = [
        os.path.join('data', 'yield_per_cell_1.tif'),
        os.path.join('data', 'yield_per_cell_2.tif'),
        os.path.join('data', 'yield_per_cell_3.tif'),
        os.path.join('data', 'yield_per_cell_4.tif'),
    ]

    if p.run_this:
        for raster_path in p.input_paths:
            # NOTE 2: One output file per input, each guarded by an existence
            # check. Re-runs recompute ONLY what is missing. This finally kills
            # the Level 2 pain: change one input, delete its one output file,
            # and re-run — the other three are skipped.
            sum_path = os.path.join(p.cur_dir, hb.file_root(raster_path) + '_sum.csv')
            if not hb.path_exists(sum_path):
                sum_of_yield = float(np.sum(hb.as_array(raster_path)))
                with open(sum_path, 'w') as f:
                    f.write(str(sum_of_yield) + '\n')
                hb.log('Summed ' + raster_path + ': ' + str(sum_of_yield))


if __name__ == '__main__':
    # NOTE 3: Ordinary Python hygiene, but the spec makes it mandatory:
    # importing a run file must never start a run.
    p = hb.ProjectFlow()
    p.add_task(yield_summations)
    p.execute()

Run it twice and the second run finishes instantly — every task reports its outputs already exist. The results are no longer trapped in a print that scrolls away: they live as files in <project_dir>/intermediate/yield_summations/, timestamped by the filesystem and inspectable after the fact. And notice what the existence check changed about how you work: “I updated dataset 3” is no longer a reason to rerun everything — it’s rm one file and re-run. This is the same mechanism that lets the full SEALS model resume a crashed global run without recomputing days of work; you now have it in a 30-line script.

Stage 4, redone with ProjectFlow

The Old code here is the 2017 confessional from stage 4 above — worth re-skimming its seven NOTEs, because every one of them is an ancestral form of something ProjectFlow now does properly. The New version answers them NOTE by NOTE.

Old

See the “BAD CODE; DONT RUN” block at stage 4 above (kept in one place to spare me reliving it twice). Its seven NOTEs, compressed: (1) Cython for speed, (2) a homemade utility library, (3) six lines of logging boilerplate, (4) a workspace hardcoded to E:/bulk_data/..., (5) custom raster loaders, (6) do_30s_resample = False booleans guarding expensive blocks, (7) a million-iteration double loop for what numpy does in one line.

New: Paths that find themselves, and switches that live in the tree

import os
import numpy as np
import hazelbean as hb


def build_task_tree(p):
    # With more than one task, the add_task calls move out of __main__ and into
    # this named function: build_task_tree(p) is the one place that answers
    # "what is this project's pipeline?" Every run file defines it — here it is
    # three lines; in a real project it grows; in the full spec it can delegate
    # to shared library builders.
    p.resampled_carbon_task = p.add_task(resampled_carbon, run=0)
    p.cleaned_carbon_task = p.add_task(cleaned_carbon)


def resampled_carbon(p):
    # Get path is given a ref_path, relative to the searched roots (project input/,
    # base_data, then the cloud). get_path finds it wherever it lives — or
    # downloads it
    p.c_1km_path = p.get_path(os.path.join('carbon', 'c_1km.tif'))
    p.ha_per_cell_30s_path = p.get_path(os.path.join('pyramids', 'ha_per_cell_30sec.tif'))

    # Output paths are declared ABOVE the run_this gate: a disabled task still
    # publishes where its products live (from the last run it was enabled), so
    # downstream tasks keep working. That is why this task can sit at run=0.
    p.c_30s_path = os.path.join(p.cur_dir, 'c_30s.tif')

    if p.run_this:
        if not hb.path_exists(p.c_30s_path):
            hb.resample_to_match(p.c_1km_path, p.ha_per_cell_30s_path, p.c_30s_path,
                                 resample_method='average')


def cleaned_carbon(p):
    p.cleaned_carbon_path = os.path.join(p.cur_dir, 'c_30s_cleaned.tif')
    if p.run_this:
        if not hb.path_exists(p.cleaned_carbon_path):
            # vectorized numpy, one line.
            array = hb.as_array(p.c_30s_path)
            array = np.where(array < 0, -9999, array)
            hb.save_array_as_geotiff(array, p.cleaned_carbon_path, p.c_30s_path)
            hb.log('Cleaned carbon written to ' + p.cleaned_carbon_path)


if __name__ == '__main__':
    p = hb.ProjectFlow()
    build_task_tree(p)
    p.execute()

Stage 5, redone with ProjectFlow

At global 30-second resolution a raster no longer fits in memory, and the stage 5 answer was a homemade “Tile Reference” class with manual chunk bookkeeping. The modern answer is one library call.

Old

See the “BAD CODE; DON’T RUN” block at stage 5 above: gu.Tr() hand-slices the raster into tiles, nested loops walk each tile’s rows and columns, and the chunk arithmetic (tile[0] / factor + row…) is exactly the kind of code where an off-by-one silently corrupts a publication figure.

New: One call streams the chunks

The Stage 4 tasks are unchanged; this rung adds one task to the tree, so only the new pieces are shown.

def build_task_tree(p):
    p.resampled_carbon_task = p.add_task(resampled_carbon, run=0)
    p.cleaned_carbon_task = p.add_task(cleaned_carbon)
    p.carbon_per_cell_task = p.add_task(carbon_per_cell)


def carbon_per_cell(p):
    # hb.raster_calculator streams the computation chunk by chunk
    # across all cores (dask under the hood). Pass the input paths and the
    # function to apply; the full rasters are never loaded at once. The op
    # converts carbon per hectare to carbon per cell — elementwise, so chunking
    # cannot change the result.
    p.carbon_per_cell_path = os.path.join(p.cur_dir, 'carbon_per_cell.tif')
    if p.run_this:

        def op(carbon_per_ha, ha_per_cell):
            # The operation applied to each chunk: receives arrays (one per
            # input raster) and returns an array.
            return carbon_per_ha * ha_per_cell

        if not hb.path_exists(p.carbon_per_cell_path):
            hb.raster_calculator(
                [p.cleaned_carbon_path, p.ha_per_cell_30s_path],
                op,
                p.carbon_per_cell_path)
            hb.log('Carbon per cell written to ' + p.carbon_per_cell_path)


if __name__ == '__main__':
    p = hb.ProjectFlow()
    build_task_tree(p)
    p.execute()

Notice the inputs: p.cleaned_carbon_path and p.ha_per_cell_30s_path were published by the Stage 4 tasks (above their run_this gates), so this task chains onto them without knowing or caring whether they ran this time or last week. When chunk-by-chunk parallelism isn’t enough — when the algorithm itself must run independently per zone, thousands of times — the same idea scales up from chunks to tasks via ProjectFlow’s iterators (p.add_iterator, run_in_parallel=1), which is how SEALS trains and allocates over millions of tiles. That machinery is beyond this walkthrough; the point of this rung is that you should never write chunk bookkeeping by hand again.

Stage 6, redone with ProjectFlow

For the final rung we return to the yield project from Stages 1–3, because Stage 2 left a promissory note there: “this list is the seed of what will later become a scenarios.csv.” Time to pay it.

Old

The Old is the Stage 3 script itself — nothing is wrong with its code. The pain is now operational: your adviser asks for a fifth dataset, or the same four for a different year, and the varying thing is hardcoded inside a task, so varying it means editing code. You want a quick smoke version versus the real run, and the Level 2 instinct — copy the file, edit one line — is exactly the retraction machine. And a collaborator wants to run it on their machine, where your paths and credentials don’t exist.

New: The run is driven by data — scenarios, parameters, and a thin test wrapper

The canonical run file separates three things that were tangled together: scenarios (the rows of work this run iterates over — what varies), parameters (values constant across all scenarios — including machine-specific ones), and code (which no longer changes when either of those do). Both CSVs live in a tracked input_template/ directory next to the run file; constructing the ProjectFlow copies anything missing into the project’s untracked input/ working copy on first run, and the working copy is never overwritten — so per-machine values survive re-runs.

input_template/yield_scenarios.csv — one row per unit of work:

scenario_label	yield_path	year
cell_1	data/yield_per_cell_1.tif	2017
cell_2	data/yield_per_cell_2.tif	2017
cell_3	data/yield_per_cell_3.tif	2017
cell_4	data/yield_per_cell_4.tif	2017

input_template/yield_parameters.csv — vertical key,value, constant per run. Machine-specific keys ship blank in the tracked template; each machine fills its own untracked input/ copy:

key	value
ndv	-9999
data_credentials_path

run_yield.py:

"""Sum caloric yield for every scenario in the scenarios CSV."""
import os
import numpy as np
import pandas as pd
import hazelbean as hb


def build_task_tree(p):
    p.yield_summations_task = p.add_task(yield_summations)


def yield_summations(p):
    if p.run_this:
        for _, row in p.scenarios_df.iterrows():
            # One output file per scenario ROW, each behind an existence check —
            # the same skip logic as Stage 3, now driven by the CSV. A finished
            # scenario is skipped on re-run; a finished TASK (all rows done)
            # costs only these existence checks. Delete one file to redo one
            # scenario; add a row to the CSV to extend the run.
            sum_path = os.path.join(p.cur_dir, row['scenario_label'] + '_sum.csv')
            if not hb.path_exists(sum_path):
                array = hb.as_array(row['yield_path'])
                total = float(array[array != p.ndv].sum())  # p.ndv came from parameters.csv
                with open(sum_path, 'w') as f:
                    f.write('scenario_label,year,total\n')
                    f.write(row['scenario_label'] + ',' + str(row['year']) + ',' + str(total) + '\n')
                hb.log('Summed ' + row['scenario_label'] + ': ' + str(total))


def run_project(scenario_definitions_filename='yield_scenarios.csv',
                project_name='yield',
                run_mode='check',
                tasks_to_skip=None):
    """The full run and the test run differ ONLY by which scenarios CSV drives
    the task tree. A stable project_name (run_mode='check') means repeated runs
    resume in place, skipping anything already computed. Returns p."""
    # One call does the whole directory setup: it places the project, creates it,
    # and copies input_template/ -> input/ (skip-existing).
    p = hb.ProjectFlow(project_name=project_name, run_mode=run_mode)

    # Parameters: constant across scenarios. Hydrated onto p, so tasks just use
    # p.ndv, p.data_credentials_path, etc. Blank values read as None.
    parameters_df = pd.read_csv(os.path.join(p.input_dir, 'yield_parameters.csv'))
    for _, row in parameters_df.iterrows():
        setattr(p, row['key'], None if pd.isna(row['value']) else row['value'])
    p.ndv = float(p.ndv)

    # Scenarios: the rows of work. Which CSV loads is the run_project argument —
    # this is the ONLY thing the test wrapper changes.
    p.scenarios_df = pd.read_csv(os.path.join(p.input_dir, scenario_definitions_filename))

    build_task_tree(p)
    p.skip_tasks(tasks_to_skip)

    p.L = hb.get_logger(p.project_name)
    p.execute()
    return p


if __name__ == '__main__':
    run_project()

And the smoke test, run_yield_test.py, is the whole answer to the Level 2 copy-the-file instinct — the variant is a CSV plus a few lines, never a fork:

"""Pared test of the yield project: same task tree, a one-row scenarios CSV, stable dir."""
from run_yield import run_project

if __name__ == '__main__':
    run_project(scenario_definitions_filename='yield_scenarios_test.csv',
                project_name='yield_test')

Run run_yield.py twice: the second run loads the CSVs, checks four existence checks, and exits — everything finished is skipped, which is what makes a stable project dir the default (the test project resumes in place the same way). Add a fifth row to yield_scenarios.csv and re-run: only the new scenario computes. Change a machine value (a credentials path, an ssh host) in input/yield_parameters.csv and no code changes anywhere. That is the whole canon in miniature: code defines the pipeline once; scenarios say what to run over; parameters say what is constant; the filesystem remembers what is done.

This is where the walkthrough ends and the spec begins. The full rules — variant runs via tasks_to_skip, shared library builders, ref_paths in scenario CSVs, the run_<project>_test.py conventions — live on the conventions page; a copy-me seed for a project that imports a devstack library is run_template_downstream_user.py; and the full-scale reference implementation, running a two-pass global economy-environment pipeline with exactly this anatomy, is ngfs_pnas (a private repository — ask for access). What was deliberately not shown here: the libraries’ richer scenario hydration (rows assigned onto p and iterated through the tree), parallel iterators, and the machine-config registry — each is the same idea you just saw, grown up.