Refactor harvest pipeline: fix parallelism deadlocks and add output w… - #3
Open
evasinha-pnnl wants to merge 25 commits into
Open
Refactor harvest pipeline: fix parallelism deadlocks and add output w…#3evasinha-pnnl wants to merge 25 commits into
evasinha-pnnl wants to merge 25 commits into
Conversation
…riting - harvest.py: Switch from mp.Pool (fork, no initializer) to fork-based Pool with an initializer that sets plain numpy globals before any task runs. Workers are pure-compute functions that return regridded arrays; the main process writes results back to lt_year_data via imap_unordered — no proxy objects in worker processes, no manager deadlock. - harvest.py: Pre-read LUH2 and HYDE source data once in run() and share via pool initializer globals (copy-on-write after fork), eliminating 2520 redundant file reads per year. Build cellid_to_idx mapping to correctly index lt_year_data by positional row rather than raw HEALPix cell ID. - harvest.py: Switch from 10x10 to 5x5 degree spatial chunks (360 → ~1440 non-ocean chunks) and sort largest-first to reduce load imbalance; observed 10x10 chunk times ranged from 3s to 543s, causing long idle tails. - land_type.py: Allocate LtData with correct n_cells from HEALPix parquet mesh (was defaulting to n_cells=1). Add try/except ImportError guards around unimplemented modules (crop, urban, lake, ice, normalize_cell, veg_assoc, consistency). Write per-year NetCDF output after each year completes. - landgen_io.py: Add write_lt_year_data_to_netcdf() to write harvest_frac and grazing_frac with zlib compression to a per-year NetCDF file. Move calc_ll_limits() here from land_type.py as a shared utility, removing the circular import between harvest.py and land_type.py. Fix xr.Dataset coordinate construction syntax. - shared_data.py: Add getter/setter methods to LtData for manager proxy access (set_harvest_frac, set_grazing_frac, get_harvest_frac, etc.). - topography.py: Replace sys.exit(1) with a warning when active=False and output file is missing; remove unused sys import. - landcover.py: Stub out landcover_process body; comment out unimplemented module imports to allow the rest of the pipeline to run. - submit_landgen.sh: Add --cpus-per-task=128, use SLURM_CPUS_PER_TASK, remove --partition (not needed with --qos on Perlmutter), set MKL_NUM_THREADS=1 to prevent CPU oversubscription inside workers.
Trying to generalize some code. Starting setting up the hdf io for modis. (cherry picked from commit 9b96593)
evasinha-pnnl
force-pushed
the
evasinha/lnd/landgen
branch
from
May 7, 2026 15:54
6fc7afa to
8203393
Compare
evasinha-pnnl
force-pushed
the
evasinha/lnd/landgen
branch
from
May 7, 2026 21:08
9dd6ab4 to
0a8cdd2
Compare
evasinha-pnnl
force-pushed
the
evasinha/lnd/landgen
branch
from
May 7, 2026 21:15
0a8cdd2 to
3a48d92
Compare
harvest.run() previously recomputed spatial chunk boundaries internally using a 5x5-degree lat-lon grid built from the global parquet, creating a second chunk decomposition that was inconsistent with the tight-vertex bounding boxes computed by set_decomp_cell_idx_ll_limits() in landgen.py. Changes: - harvest.run() signature: move decomp_indices/decomp_ll_limits before manager/grid_manager to match the call convention in land_type.py. - Remove the internal chunk-building loop (5x5-degree lat-lon mask on the parquet) and the now-unused 'from . import land_type' import. - Use the pre-computed decomp_indices/decomp_ll_limits passed in by the caller, which carry tight per-chunk vertex bounding boxes — exactly the ll_limits regrid_to_landgen_grid needs so the source raster slice fully covers every polygon in the chunk. - land_type.py: update harvest.run() call to pass decomp_indices and decomp_ll_limits. - landgen_io.set_decomp_cell_idx_ll_limits(): store HEALPix cellid values (int64) in decomp_indices instead of NetCDF row indices; update docstring accordingly.
Two bugs caused cell values to be written to the wrong geographic locations in lt_year_data and the output NetCDF. Bug 1 — cellid_to_idx used sorted order instead of parquet row order: cellid_to_idx was built by sorting the global parquet cellids before enumerating them. lt_year_data arrays are allocated in parquet row order, and out_grid_data.cell_id stores cellids in parquet row order. Using sorted order made lt_year_data[k] hold the value for the k-th *sorted* cellid, while the NetCDF coordinate cell_id[k] was the k-th *parquet-order* cellid — scrambling every cell to the wrong location. HEALPix ring ordering groups same-latitude cells into consecutive cellids, so the mismatch produced the visible horizontal stripes. Fix: remove np.sort(); build cellid_to_idx directly from parquet order. Bug 2 — HYDE3.5 grazing values not converted to fractions: HYDE3.5 grazing variables are in km² per source grid cell. After area-weighted regridding the result is still in km², but was stored directly in grazing_frac without conversion, producing values >> 1. Fix: add HEALPIX_CELL_AREA_KM2 = 162.5086 km² constant and divide the regridded result by it to obtain a dimensionless fraction, then clip to [0, 1] to suppress floating-point rounding artefacts. Cleanup — rename misleading parameter in LtData.set_*() methods: set_harvest_frac(), set_harvest_mass(), set_grazing_frac(), set_pct_pft() had a parameter named cell_ids but received positional row indices. Renamed to row_indices to match actual usage.
…geotiff HYDE3.5 and LUH2 both store latitude north-to-south (lat[0] ≈ +90, decreasing). The previous code computed the GeoTIFF south/north bounds from chunk_lat[0] and chunk_lat[-1], assuming S-to-N order. With N-to-S data this swapped south and north, producing a spatially inverted GeoTIFF. uraster then assigned every regridded value to the wrong latitude — causing a systematic ~10° northward shift in all harvest and grazing output fields. Fix: use chunk_lat.min()/chunk_lat.max() (and same for lon) so the bounds are always geographically correct regardless of the source array ordering.
Collaborator
Author
|
landgen_vs_orig_data.pdf |
This just provides a current snapshot of my branch to others. It is not functioning yet. (cherry picked from commit 791c356)
The code runs successfully for the landcover submodule. For each chunk MODIS data are read, converted, regridded, writting to the output file, and plotted. This is still demonstration code to get a working structure that can be easily expanded. (cherry picked from commit 1062844)
The code runs successfully for the landcover submodule. For each chunk MODIS data are read, converted, regridded, writting to the output file, and plotted. This is still demonstration code to get a working structure that can be easily expanded.
Use the MOD44W data to separate ocean pixels. But need to include deep inland water and overlay vector coastline to include the mediterranean, black, and azov. Manually remove the capspian because it is not connected to ocean waterways and is included in the lakes data.
Collaborator
Author
|
@aldivi I have cherry-picked the latest commits from your branch |
… and grazing data
- Workers now return chunk LtData objects instead of array tuples
- Main process uses single copy_from() call to merge results (replaces nested loops with individual setters)
- Fixed grazing key consistency: uses stems ('pasture') throughout instead of mixing with filenames
- Updated docstrings to reflect new return types
- Simplifies code: 1 line per chunk vs 7+ lines, more maintainable and extensible
Unify management.py workflow with landcover.py approach: - Write mesh once per chunk instead of once per variable (10× → 1×) - Use write_mesh_to_geojson() + regrid_to_mesh() instead of monolithic regrid_to_landgen_grid() - Replace DataFrame-based mesh handling (global_mesh_df) with GridData object (out_grid_data) for consistency across modules - Add proper temp directory cleanup with try/finally + shutil.rmtree() - regrid_to_landgen_grid() and write_chunk_mesh_to_geojson() are now obsolete and and therefore removed from landgen_io.py. Performance: Reduces file I/O overhead by reusing mesh.geojson for all harvest and grazing variables in each spatial chunk.
Dropped - load_global_mesh_parquet - build_cellid_to_idx_map - get_cell_area_km2 - get_chunk_cell_ids All replaced by direct GridData access. Updated land_type.py to use out_grid_data.num_cells instead of reading parquet file.
evasinha-pnnl
force-pushed
the
evasinha/lnd/landgen
branch
from
June 30, 2026 21:07
ee0c605 to
351c203
Compare
Pre-compute row_indices once in run() using vectorized numpy operations instead of recalculating per chunk with dict lookups: - Use np.searchsorted() for O(log n) vectorized cellid→row_index conversion - Check if cell_id array is already sorted before sorting (typical case) - Remove cellid_to_idx dict from worker globals (no longer needed) - Workers receive pre-computed row_indices directly, eliminating per-chunk dictionary lookup overhead Performance: Reduces per-chunk overhead from O(n) dict lookups to one-time O(n log n) sorted array preparation shared across all chunks.
Simplify multiprocessing approach to match landcover.py: - Remove fork context + initializer complexity (global worker state) - Workers now read their own data via read_netcdf_ll with ll_limits - Use simple pool.starmap() instead of fork + imap_unordered - Remove unused locks (no longer needed without shared state) Benefits: simpler code, more portable (works on all platforms), easier to debug. Trade-off: workers read data subsets N times instead of once, but code maintainability gains outweigh the I/O cost.
Remove leftover code from refactoring: - Remove shapely.wkb import from landgen_io.py (unused after removing write_chunk_mesh_to_geojson) - Remove 'manager' parameter from management.py run() (no longer needed with starmap pattern) - Remove unused multiprocessing and shared_data module imports from management.py and land_type.py - Update management.run() call site in land_type.py to match new signature Fix missing output variables in land_type.py: - Add harvest_frac and grazing_frac to varnames list (were computed but not written to NetCDF) - Add harvest_frac and grazing_frac to timevars list (for annual file concatenation) - Add documentation comments explaining output variables
Per original design, decomp_indices already contains 0-based row indices (not arbitrary cellid values). Remove unnecessary searchsorted conversion logic added during earlier refactoring. Add explicit assumption documentation: this approach requires mesh cellid values to be sequential (0, 1, 2, ..., n-1) matching their NetCDF row positions. Non-sequential cellids will cause array indexing errors. Matches landcover.py pattern which also uses decomp_indices directly.
evasinha-pnnl
force-pushed
the
evasinha/lnd/landgen
branch
from
July 15, 2026 03:57
b75d47f to
70e538b
Compare
Owner
|
@evasinha-pnnl interestingly enough, the change to the grazing names in your last commit failed in my run. copilot suggested a fix that changes it back to how it was before. I am rerunning it now. |
This is a rebase to the remote so that these updates can be pushed. Pushing these updates to this branch before moving all the code to the new repo, for safety. Some efficiencies have been added. The plot function now has a scale-limits option. The grid file is read only once now. A few other things have been cleaned up. Some grid variables have been renamed to be more clear. This has been tested with both landcover and management processing on. This is ready to be converted to a main branch in a dedicated repo. On branch evasinha/lnd/landgen Changes to be committed: modified: components/elm/tools/landgen/README.md modified: components/elm/tools/landgen/config_template.json modified: components/elm/tools/landgen/load_landgen_env.sh modified: components/elm/tools/landgen/src/landgen/__main__.py modified: components/elm/tools/landgen/src/landgen/land_type.py modified: components/elm/tools/landgen/src/landgen/landcover.py modified: components/elm/tools/landgen/src/landgen/landgen.py modified: components/elm/tools/landgen/src/landgen/landgen_io.py modified: components/elm/tools/landgen/src/landgen/management.py modified: components/elm/tools/landgen/src/landgen/plot_landgen.py modified: components/elm/tools/landgen/src/landgen/shared_data.py modified: components/elm/tools/landgen/submit_landgen.sh
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
…riting
harvest.py: Switch from mp.Pool (fork, no initializer) to fork-based Pool with an initializer that sets plain numpy globals before any task runs. Workers are pure-compute functions that return regridded arrays; the main process writes results back to lt_year_data via imap_unordered — no proxy objects in worker processes, no manager deadlock.
harvest.py: Pre-read LUH2 and HYDE source data once in run() and share via pool initializer globals (copy-on-write after fork), eliminating 2520 redundant file reads per year. Build cellid_to_idx mapping to correctly index lt_year_data by positional row rather than raw HEALPix cell ID.
harvest.py: Switch from 10x10 to 5x5 degree spatial chunks (360 → ~1440 non-ocean chunks) and sort largest-first to reduce load imbalance; observed 10x10 chunk times ranged from 3s to 543s, causing long idle tails.
land_type.py: Allocate LtData with correct n_cells from HEALPix parquet mesh (was defaulting to n_cells=1). Add try/except ImportError guards around unimplemented modules (crop, urban, lake, ice, normalize_cell, veg_assoc, consistency). Write per-year NetCDF output after each year completes.
landgen_io.py: Add write_lt_year_data_to_netcdf() to write harvest_frac and grazing_frac with zlib compression to a per-year NetCDF file. Move calc_ll_limits() here from land_type.py as a shared utility, removing the circular import between harvest.py and land_type.py. Fix xr.Dataset coordinate construction syntax.
shared_data.py: Add getter/setter methods to LtData for manager proxy access (set_harvest_frac, set_grazing_frac, get_harvest_frac, etc.).
topography.py: Replace sys.exit(1) with a warning when active=False and output file is missing; remove unused sys import.
landcover.py: Stub out landcover_process body; comment out unimplemented module imports to allow the rest of the pipeline to run.
submit_landgen.sh: Add --cpus-per-task=128, use SLURM_CPUS_PER_TASK, remove --partition (not needed with --qos on Perlmutter), set MKL_NUM_THREADS=1 to prevent CPU oversubscription inside workers.