From 5dd5ddd1b2f40241c477b9d843c883a02dc1401d Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 3 Dec 2025 12:55:38 -0700 Subject: [PATCH 01/61] Update --- revrt/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/revrt/__init__.py b/revrt/__init__.py index 71cfa24c..6dc698af 100644 --- a/revrt/__init__.py +++ b/revrt/__init__.py @@ -1,4 +1,4 @@ -"""Routing analysis library for the reV model""" +"""Routing analysis library""" import importlib.metadata From a30e6686fb89c4188a6332c89b5f2522240e78b7 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 3 Dec 2025 13:01:13 -0700 Subject: [PATCH 02/61] Update --- .github/ISSUE_TEMPLATE/bug_report.yaml | 2 +- .github/ISSUE_TEMPLATE/feature_request.yaml | 2 +- README.rst | 2 +- docs/source/conf.py | 4 ++-- docs/source/glossary.rst | 7 +------ pixi.lock | 2 +- pyproject.toml | 4 ++-- 7 files changed, 9 insertions(+), 14 deletions(-) diff --git a/.github/ISSUE_TEMPLATE/bug_report.yaml b/.github/ISSUE_TEMPLATE/bug_report.yaml index e36ebc81..1c1eea2e 100644 --- a/.github/ISSUE_TEMPLATE/bug_report.yaml +++ b/.github/ISSUE_TEMPLATE/bug_report.yaml @@ -77,4 +77,4 @@ body: attributes: label: Charge code placeholder: > - If you are at NREL and fixing this bug is urgent, please provide a charge code + If you are at NLR and fixing this bug is urgent, please provide a charge code diff --git a/.github/ISSUE_TEMPLATE/feature_request.yaml b/.github/ISSUE_TEMPLATE/feature_request.yaml index 5f305b31..9296bdd1 100644 --- a/.github/ISSUE_TEMPLATE/feature_request.yaml +++ b/.github/ISSUE_TEMPLATE/feature_request.yaml @@ -59,4 +59,4 @@ body: attributes: label: Charge code placeholder: > - If you are at NREL and implementing this feature is urgent, please provide a charge code + If you are at NLR and implementing this feature is urgent, please provide a charge code diff --git a/README.rst b/README.rst index 5d700493..ea48261a 100644 --- a/README.rst +++ b/README.rst @@ -20,7 +20,7 @@ Welcome to reV Routing (reVRt)! .. |Pixi| image:: https://img.shields.io/endpoint?url=https://raw.githubusercontent.com/prefix-dev/pixi/main/assets/badge/v0.json :target: https://pixi.sh -.. |SWR| image:: https://img.shields.io/badge/SWR--25--112_-blue?label=NREL +.. |SWR| image:: https://img.shields.io/badge/SWR--25--112_-blue?label=NLR :alt: Static Badge .. |Zenodo| image:: https://zenodo.org/badge/944738283.svg diff --git a/docs/source/conf.py b/docs/source/conf.py index bd7d82d8..7b29bce0 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -21,7 +21,7 @@ project = "reVRt" copyright = "2025, Alliance for Sustainable Energy, LLC" -author = "NREL: Guilherme Pimenta Castelao, Paul Pinchuk" +author = "NLR: Guilherme Pimenta Castelao, Paul Pinchuk" pkg = os.path.dirname(os.path.abspath(os.path.dirname(__file__))) pkg = os.path.dirname(pkg) @@ -135,7 +135,7 @@ html_context = { "display_github": True, - "github_user": "nrel", + "github_user": "nlr", "github_repo": "reVRt", "github_version": "main", "conf_py_path": "/docs/source/", diff --git a/docs/source/glossary.rst b/docs/source/glossary.rst index ecb48faf..44230e44 100644 --- a/docs/source/glossary.rst +++ b/docs/source/glossary.rst @@ -119,7 +119,7 @@ supporting technologies referenced throughout the reVRt documentation. ``revrt._rust``. reV - NREL's Renewable Energy Potential platform that models generation + NLR's Renewable Energy Potential platform that models generation profiles for renewable resources. reVRt complements reV by adding transmission routing capabilities that feed downstream analyses. @@ -128,11 +128,6 @@ supporting technologies referenced throughout the reVRt documentation. simulations, and exports supporting analytics for transmission planning studies. - reVX - NREL's valuation framework that builds on reV outputs to assess - market and financial metrics. reVRt outputs can flow into reVX to - inform siting and transmission decisions. - routing file HDF5 file produced by the cost preparation pipeline that stores coordinate grids, friction and barrier layers, metadata, and any diff --git a/pixi.lock b/pixi.lock index 1502115f..1aa74c88 100644 --- a/pixi.lock +++ b/pixi.lock @@ -16885,7 +16885,7 @@ packages: - pypi: ./ name: nrel-revrt version: 0.2.1 - sha256: b79b7d6fe54b011d9b70599a4d70c8fc316d10124946018122ff94c4af90eca5 + sha256: 283ea597b21e2a863365183e656f8e6f7ccf29eeaae1753dc008d9614aa2ade5 requires_dist: - affine>=2.4.0,<3 - dask>=2025.4.1,<2026 diff --git a/pyproject.toml b/pyproject.toml index 124076ef..35d5e15d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ build-backend = "maturin" [project] name = "NREL-reVRt" version = "0.2.1" -description = "National Renewable Energy Laboratory's (NREL's) Transmission for reVX tool" +description = "National Laboratory of the Rockies' (NLR's) Transmission for reVX tool" readme = "README.rst" authors = [ {name = "Guilherme Castelão", email = "gpimenta@nrel.gov"}, @@ -16,7 +16,7 @@ maintainers = [ {name = "Paul Pinchuk", email = "ppinchuk@nrel.gov"}, ] license = {text = "BSD-3-Clause"} -keywords = ["NREL", "routing", "transmission", "reVRt", "reV", "reVX"] +keywords = ["NLR", "routing", "transmission", "reVRt", "reV", "reVX"] requires-python = ">= 3.11" classifiers=[ "Development Status :: 2 - Pre-Alpha", From 515b2ac7ebd63f058be813e7a0edafa4c5750746 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 14:52:43 -0700 Subject: [PATCH 03/61] Update lockfile --- pixi.lock | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/pixi.lock b/pixi.lock index 1aa74c88..dc8a008e 100644 --- a/pixi.lock +++ b/pixi.lock @@ -2606,6 +2606,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/linux-64/ruff-0.14.7-h813ae00_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/ruff-lsp-0.0.62-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/rust-1.90.0-h53717f1_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/rust-src-1.90.0-unix_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-unknown-linux-gnu-1.90.0-h2c6d0dc_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/s2n-1.6.0-h8399546_1.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/scikit-image-0.25.2-py313h08cd8bf_2.conda @@ -3070,6 +3071,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-64/ruff-0.14.7-hd9f4cfa_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/ruff-lsp-0.0.62-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/rust-1.90.0-h34a2095_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/rust-src-1.90.0-unix_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-apple-darwin-1.90.0-h38e4360_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/scikit-image-0.25.2-py313h2f264a9_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/scikit-learn-1.7.2-py313hfce3ed8_0.conda @@ -3510,6 +3512,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/osx-arm64/ruff-0.14.7-h382de68_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/ruff-lsp-0.0.62-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/rust-1.90.0-h4ff7c5d_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/rust-src-1.90.0-unix_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.90.0-hf6ec828_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scikit-image-0.25.2-py313h7d16b84_2.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/scikit-learn-1.7.2-py313h6eea1c5_0.conda @@ -3950,6 +3953,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/win-64/ruff-0.14.7-h15e3a1f_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/ruff-lsp-0.0.62-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/rust-1.90.0-hf8d6059_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/rust-src-1.90.0-win_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-x86_64-pc-windows-msvc-1.90.0-h17fc481_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/scikit-image-0.25.2-py313hc90dcd4_2.conda - conda: https://conda.anaconda.org/conda-forge/win-64/scikit-learn-1.7.2-py313he28f1d7_0.conda @@ -16885,7 +16889,7 @@ packages: - pypi: ./ name: nrel-revrt version: 0.2.1 - sha256: 283ea597b21e2a863365183e656f8e6f7ccf29eeaae1753dc008d9614aa2ade5 + sha256: 8db001fdf27ec3ffa8b7a17346a3e9965d791d1cacdff0da92b63d0544a8a25a requires_dist: - affine>=2.4.0,<3 - dask>=2025.4.1,<2026 @@ -20045,6 +20049,30 @@ packages: purls: [] size: 251268829 timestamp: 1758352792592 +- conda: https://conda.anaconda.org/conda-forge/noarch/rust-src-1.90.0-unix_0.conda + sha256: 02f3d1a707b62b91ef2ae0f9a810e5e819f8b0559bfbda8e709545da23d568d3 + md5: 0cd189e71b846188cca9b789de366369 + depends: + - __unix + constrains: + - rust >=1.90.0,<1.90.1.0a0 + license: MIT + license_family: MIT + purls: [] + size: 4099159 + timestamp: 1758349704068 +- conda: https://conda.anaconda.org/conda-forge/noarch/rust-src-1.90.0-win_0.conda + sha256: 04aaf10bee07340c0cc9f05a6bcf90ac10f3797554fedf73ad1d5b73d4d8a9d8 + md5: 49b13cc66e9ab3f9d88985df52afcad9 + depends: + - __win + constrains: + - rust >=1.90.0,<1.90.1.0a0 + license: MIT + license_family: MIT + purls: [] + size: 4081987 + timestamp: 1758351259295 - conda: https://conda.anaconda.org/conda-forge/noarch/rust-std-aarch64-apple-darwin-1.90.0-hf6ec828_0.conda sha256: 3c39f49325d0335c24ded6c2944734bb1ca3b072255a0a20942e4a0970065ad2 md5: 051ab4f9607160ed71ad648cf0f222ab From a3044b0dace75c4f9c44fc1b49ca189fbb420dba Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 15:10:13 -0700 Subject: [PATCH 04/61] Add new exception type --- revrt/exceptions.py | 4 ++++ tests/python/unit/test_exceptions.py | 5 +++++ 2 files changed, 9 insertions(+) diff --git a/revrt/exceptions.py b/revrt/exceptions.py index 1921e1c8..0c5ffedd 100644 --- a/revrt/exceptions.py +++ b/revrt/exceptions.py @@ -33,6 +33,10 @@ class revrtFileNotFoundError(revrtError, FileNotFoundError): """revrt FileNotFoundError""" +class revrtInvalidStartCostError(revrtError, ValueError): + """revrt FileNotFoundError""" + + class revrtKeyError(revrtError, KeyError): """revrt KeyError""" diff --git a/tests/python/unit/test_exceptions.py b/tests/python/unit/test_exceptions.py index c77dbf98..bace3e9d 100644 --- a/tests/python/unit/test_exceptions.py +++ b/tests/python/unit/test_exceptions.py @@ -10,6 +10,7 @@ revrtConfigurationError, revrtFileExistsError, revrtFileNotFoundError, + revrtInvalidStartCostError, revrtKeyError, revrtNotImplementedError, revrtProfileCheckError, @@ -68,6 +69,10 @@ def test_exceptions_log_uncaught_error(assert_message_was_logged): revrtFileNotFoundError, [revrtError, revrtFileNotFoundError, FileNotFoundError], ), + ( + revrtInvalidStartCostError, + [revrtError, revrtInvalidStartCostError, ValueError], + ), ( revrtNotImplementedError, [revrtError, revrtNotImplementedError, NotImplementedError], From b0e4dc31a57c85dd7e513404f473515f3ed16e6b Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 15:12:03 -0700 Subject: [PATCH 05/61] Add another error type --- revrt/exceptions.py | 6 +++++- tests/python/unit/test_exceptions.py | 5 +++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/revrt/exceptions.py b/revrt/exceptions.py index 0c5ffedd..50b3f252 100644 --- a/revrt/exceptions.py +++ b/revrt/exceptions.py @@ -34,13 +34,17 @@ class revrtFileNotFoundError(revrtError, FileNotFoundError): class revrtInvalidStartCostError(revrtError, ValueError): - """revrt FileNotFoundError""" + """revrt InvalidStartCostError""" class revrtKeyError(revrtError, KeyError): """revrt KeyError""" +class revrtLeastCostPathNotFoundError(revrtError, RuntimeError): + """revrt LeastCostPathNotFoundError""" + + class revrtNotImplementedError(revrtError, NotImplementedError): """revrt NotImplementedError""" diff --git a/tests/python/unit/test_exceptions.py b/tests/python/unit/test_exceptions.py index bace3e9d..d00bc28f 100644 --- a/tests/python/unit/test_exceptions.py +++ b/tests/python/unit/test_exceptions.py @@ -12,6 +12,7 @@ revrtFileNotFoundError, revrtInvalidStartCostError, revrtKeyError, + revrtLeastCostPathNotFoundError, revrtNotImplementedError, revrtProfileCheckError, revrtRuntimeError, @@ -73,6 +74,10 @@ def test_exceptions_log_uncaught_error(assert_message_was_logged): revrtInvalidStartCostError, [revrtError, revrtInvalidStartCostError, ValueError], ), + ( + revrtLeastCostPathNotFoundError, + [revrtError, revrtLeastCostPathNotFoundError, RuntimeError], + ), ( revrtNotImplementedError, [revrtError, revrtNotImplementedError, NotImplementedError], From 60c1fca7cb24876c7317a1c71da5563284552992 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 15:14:41 -0700 Subject: [PATCH 06/61] Log class names --- revrt/exceptions.py | 4 +++- revrt/warn.py | 4 +++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/revrt/exceptions.py b/revrt/exceptions.py index 50b3f252..1a5ddc55 100644 --- a/revrt/exceptions.py +++ b/revrt/exceptions.py @@ -14,7 +14,9 @@ def __init__(self, *args, **kwargs): """Init exception and broadcast message to logger""" super().__init__(*args, **kwargs) if args: - logger.error(str(args[0]), stacklevel=2) + logger.error( + "<%s> %s", self.__class__.__name__, args[0], stacklevel=2 + ) class revrtAttributeError(revrtError, AttributeError): diff --git a/revrt/warn.py b/revrt/warn.py index 0d6755d3..1b2f642a 100644 --- a/revrt/warn.py +++ b/revrt/warn.py @@ -14,7 +14,9 @@ def __init__(self, *args, **kwargs): """Init exception and broadcast message to logger""" super().__init__(*args, **kwargs) if args: - logger.warning(str(args[0]), stacklevel=2) + logger.warning( + "<%s> %s", self.__class__.__name__, args[0], stacklevel=2 + ) class revrtDeprecationWarning(revrtWarning, DeprecationWarning): From c8809454d714318e7573875a2694d1583c3fdb4d Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 15:14:48 -0700 Subject: [PATCH 07/61] Update tests to check for name --- tests/python/unit/test_exceptions.py | 1 + tests/python/unit/test_warn.py | 1 + 2 files changed, 2 insertions(+) diff --git a/tests/python/unit/test_exceptions.py b/tests/python/unit/test_exceptions.py index d00bc28f..daf63f9c 100644 --- a/tests/python/unit/test_exceptions.py +++ b/tests/python/unit/test_exceptions.py @@ -101,6 +101,7 @@ def test_catching_error_by_type( raise raise_type(BASIC_ERROR_MESSAGE) assert BASIC_ERROR_MESSAGE in str(exc_info.value) + assert_message_was_logged(raise_type.__name__, "ERROR") assert_message_was_logged(BASIC_ERROR_MESSAGE, "ERROR") diff --git a/tests/python/unit/test_warn.py b/tests/python/unit/test_warn.py index fa2f0a55..c307c45a 100644 --- a/tests/python/unit/test_warn.py +++ b/tests/python/unit/test_warn.py @@ -18,6 +18,7 @@ def test_warnings_log_message(warning_class, assert_message_was_logged): """Test that a raised warning logs message, if any""" warn(BASIC_WARNING_MESSAGE, warning_class) + assert_message_was_logged(warning_class.__name__, "WARNING") assert_message_was_logged(BASIC_WARNING_MESSAGE, "WARNING") From 160190774dd3233d0cefdf039375e6271ac43e52 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 20:57:30 -0700 Subject: [PATCH 08/61] Update logic to look for shape dataset --- crates/revrt/src/dataset.rs | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/crates/revrt/src/dataset.rs b/crates/revrt/src/dataset.rs index 71f5c3aa..4281f255 100644 --- a/crates/revrt/src/dataset.rs +++ b/crates/revrt/src/dataset.rs @@ -80,9 +80,16 @@ impl Dataset { let name = entry.split('/').next().unwrap_or("").to_ascii_lowercase(); const EXCLUDES: [&str; 6] = ["latitude", "longitude", "band", "x", "y", "spatial_ref"]; - !name.ends_with(".json") && !EXCLUDES.iter().any(|needle| name.contains(needle)) + !name.ends_with(".json") && !EXCLUDES.iter().any(|needle| name == *needle) }) - .expect("no suitable variables found in source dataset"); + .unwrap_or_else(|| { + panic!( + "no non-coordinate variables found in source dataset: {:?}", + source + .list() + .unwrap_or_else(|err| panic!("failed to list dataset entries: {err}")) + ) + }); // Skip coordinate axes when selecting a representative variable for cost storage. let varname = first_entry.split('/').next().unwrap().to_string(); debug!("Using '{}' to determine shape of cost data", varname); From 74bd763766e0ea9e2b46ce98d031961d589cba3e Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 20:57:51 -0700 Subject: [PATCH 09/61] Add first few tests for routing --- .../routing/test_routing_point_to_many.py | 182 ++++++++++++++++++ 1 file changed, 182 insertions(+) create mode 100644 tests/python/unit/routing/test_routing_point_to_many.py diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py new file mode 100644 index 00000000..444945d0 --- /dev/null +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -0,0 +1,182 @@ +"""reVrt tests for routing one point to many endpoints""" + +from pathlib import Path + +import pytest +import numpy as np +import xarray as xr +from rasterio.transform import from_origin + +from revrt.utilities import LayeredFile +from revrt.routing.point_to_many import find_all_routes, RoutingScenario + + +@pytest.fixture(scope="module") +def sample_layered_data(tmp_path_factory): + """Sample layered data files to use across tests""" + data_dir = tmp_path_factory.mktemp("routing_data") + + layered_fp = data_dir / "test_layered.zarr" + layer_file = LayeredFile(layered_fp) + + height, width = (7, 8) + cell_size = 1.0 + x0, y0 = 0.0, float(height) + transform = from_origin(x0, y0, cell_size, cell_size) + x_coords = ( + x0 + np.arange(width, dtype=np.float32) * cell_size + cell_size / 2 + ) + y_coords = ( + y0 - np.arange(height, dtype=np.float32) * cell_size - cell_size / 2 + ) + + layer_1 = np.array( + [ + [ + [7, 7, 8, 0, 9, 9, 9, 0], + [8, 1, 2, 2, 9, 9, 9, 0], + [9, 1, 3, 3, 9, 1, 2, 3], + [9, 1, 2, 1, 9, 1, 9, 0], + [9, 9, 9, 1, 9, 1, 9, 0], + [9, 9, 9, 1, 1, 1, 9, 0], + [9, 9, 9, 9, 9, 9, 9, 0], + ] + ], + dtype=np.float32, + ) + + layer_2 = np.array( + [ + [ + [8, 7, 6, 5, 5, 6, 7, 9], + [7, 1, 1, 2, 3, 3, 2, 8], + [6, 2, 9, 6, 5, 2, 1, 7], + [7, 3, 8, 1, 2, 3, 2, 6], + [8, 4, 7, 2, 8, 4, 3, 5], + [9, 5, 6, 3, 4, 4, 3, 4], + [9, 6, 7, 4, 5, 5, 4, 3], + ] + ], + dtype=np.float32, + ) + + layer_3 = np.array( + [ + [ + [6, 6, 6, 6, 6, 7, 8, 9], + [5, 2, 2, 3, 4, 5, 6, 8], + [4, 3, 7, 7, 6, 4, 5, 7], + [5, 4, 6, 2, 3, 4, 4, 6], + [6, 5, 5, 3, 7, 5, 5, 5], + [7, 6, 6, 4, 5, 5, 4, 4], + [8, 7, 7, 5, 6, 5, 4, 3], + ] + ], + dtype=np.float32, + ) + + layer_4 = np.array( + [ + [ + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + ] + ], + dtype=np.float32, + ) + + layer_5 = np.array( + [ + [ + [0, 0, 0, 1, 1, 1, 1, 1], + [0, 0, 0, 1, 1, 1, 1, 1], + [0, 0, 0, 1, 1, 1, 1, 1], + [0, 1, 1, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 1, 0], + [0, 0, 0, 1, 0, 1, 1, 0], + [0, 0, 0, 1, 0, 0, 0, 0], + ] + ], + dtype=np.float32, + ) + + for ind, routing_layer in enumerate( + [ + layer_1, + layer_2, + layer_3, + layer_4, + layer_5, + ], + start=1, + ): + da = xr.DataArray( + routing_layer, + dims=("band", "y", "x"), + coords={"y": y_coords, "x": x_coords}, + ) + da = da.rio.write_crs("EPSG:4326") + da = da.rio.write_transform(transform) + + geotiff_fp = data_dir / f"layer_{ind}.tif" + da.rio.to_raster(geotiff_fp, driver="GTiff") + + layer_file.write_geotiff_to_file( + geotiff_fp, f"layer_{ind}", overwrite=True + ) + return layered_fp + + +def test_basic_single_route_layered_file_short_path(sample_layered_data): + """Test routing using a LayeredFile-generated cost surface""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + output = find_all_routes( + scenario, route_definitions=[((1, 1), [(1, 2)])], save_paths=False + ) + + assert len(output) == 1 + assert output.iloc[0]["cost"] == pytest.approx((1 + 2) / 2) + assert output.iloc[0]["length_km"] == 1 / 1000 + assert output.iloc[0]["cost"] == output.iloc[0]["optimized_objective"] + + +def test_basic_single_route_layered_file(sample_layered_data): + """Test routing using a LayeredFile-generated cost surface""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + output = find_all_routes( + scenario, + route_definitions=[((1, 1), [(2, 6)]), ((1, 2), [(2, 6)])], + save_paths=False, + ) + + assert len(output) == 2 + assert output.iloc[0]["cost"] == pytest.approx(11.192389) + assert output.iloc[0]["length_km"] == pytest.approx(0.0090710678) + assert np.isclose( + output.iloc[0]["cost"], output.iloc[0]["optimized_objective"] + ) + + assert output.iloc[1]["cost"] == pytest.approx(12.278174) + assert output.iloc[1]["length_km"] == pytest.approx(0.008656854) + assert np.isclose( + output.iloc[1]["cost"], output.iloc[1]["optimized_objective"] + ) + + +if __name__ == "__main__": + pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) From 6fec6e92e09f9e100d4a3121f690ffa01ddefe64 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 20:57:58 -0700 Subject: [PATCH 10/61] Routing WIP!! --- revrt/routing/point_to_many.py | 539 +++++++++++++++++++++++++++++++++ 1 file changed, 539 insertions(+) create mode 100644 revrt/routing/point_to_many.py diff --git a/revrt/routing/point_to_many.py b/revrt/routing/point_to_many.py new file mode 100644 index 00000000..9b41f851 --- /dev/null +++ b/revrt/routing/point_to_many.py @@ -0,0 +1,539 @@ +"""reVRt routing from a point to many points""" + +import json +import time +import logging +from warnings import warn + +import rasterio +import numpy as np +import pandas as pd +import xarray as xr +import dask.array as da +import geopandas as gpd +from shapely.geometry import Point +from shapely.geometry.linestring import LineString + +from revrt import find_paths + +# from revrt.exceptions import ( +# # revrtValueError, +# # revrtInvalidStartCostError, +# revrtLeastCostPathNotFoundError, +# ) +from revrt.warn import revrtWarning + +logger = logging.getLogger(__name__) +LCP_AGG_COST_LAYER_NAME = "lcp_agg_costs" +"""Special name reserved for internally-built cost layer""" + + +class RoutingScenario: + def __init__( + self, + cost_fpath, + cost_layers, + friction_layers=None, + tracked_layers=None, + cost_multiplier_layer=None, + cost_multiplier_scalar=1, + use_hard_barrier=True, + ): + self.cost_fpath = cost_fpath + self.cost_layers = cost_layers + self.friction_layers = friction_layers or [] + self.tracked_layers = tracked_layers or {} + self.cost_multiplier_layer = cost_multiplier_layer + self.cost_multiplier_scalar = cost_multiplier_scalar + self.use_hard_barrier = use_hard_barrier + + +class RoutingLayers: + """Class to build a routing layer from user input""" + + SOFT_BARRIER_MULTIPLIER = 100 + """Multiplier to apply to max cost to use for barriers + + This value is only used if ``use_hard_barrier=False``. + """ + + def __init__(self, routing_scenario, chunks="auto"): + """ + + Parameters + ---------- + cost_fpath : str + Full path to HDF5 file containing cost arrays. The cost data + layers should be named ``"tie_line_costs_{capacity}MW"``, + where ``capacity`` is an integer value representing the + capacity of the line (the integer values must matches at + least one of the integer capacity values represented by the + keys in the ``power_classes`` portion of the transmission + config). + cell_size : int, optional + Side length of each cell, in meters. Cells are assumed to be + square. By default, :obj:`CELL_SIZE`. + use_hard_barrier : bool, optional + Optional flag to treat any cost values of <= 0 as a hard + barrier (i.e. no paths can ever cross this). If ``False``, + cost values of <= 0 are set to a large value to simulate a + strong but permeable barrier. By default, ``True``. + """ + self.routing_scenario = routing_scenario + self._layer_fh = xr.open_dataset( + self.routing_scenario.cost_fpath, + chunks=chunks, + consolidated=False, + engine="zarr", + ) + # self.cost_layer_map = {} + self.li_cost_layer_map = {} + # self.friction_layer_map = {} + self.tracked_layers = [] + + self.transform = self._layer_fh.rio.transform() + self._full_shape = self._layer_fh.rio.shape + self.cost_crs = self._layer_fh.rio.crs + # self.cell_size = abs(self.transform.a) + self._layers = set(self._layer_fh.variables) + + self.cost = None + self.li_cost = None + + @property + def latitudes(self): + return self._layer_fh["latitude"] + + @property + def longitudes(self): + return self._layer_fh["longitude"] + + def _verify_layer_exists(self, layer_name): + """Verify that layer exists in cost file""" + if layer_name not in self._layers: + msg = ( + f"Did not find layer {layer_name!r} in cost " + f"file {str(self.routing_scenario.cost_fpath)!r}" + ) + logger.error(msg) + raise KeyError(msg) + + def close(self): + """Close the underlying xarray file handle""" + self._layer_fh.close() + + def build(self): + """Extract clipped cost arrays from exclusion files + + Parameters + ---------- + cost_layers : List[str] + List of layers in H5 that are summed to determine total + costs raster used for routing. Costs and distances for each + individual layer are also reported (e.g. wet and dry costs). + determining path using main cost layer. + friction_layers : List[dict] | None, optional + List of layers in H5 to be added to the cost raster to + influence routing but NOT reported in final cost. Should + have the same format as the `cost_layers` input. + By default, ``None``. + tracked_layers : dict, optional + Dictionary mapping layer names to strings, where the strings + are numpy methods that should be applied to the layer along + the LCP. For example, + ``tracked_layers={'layer_1': 'mean', 'layer_2': 'max}`` + would report the average of ``layer_1`` values along the + least cost path and the max of ``layer_2`` values along the + least cost path. Examples of numpy methods (non-exhaustive): + + - mean + - max + - min + - mode + - median + - std + + By default, ``None``, which does not track any extra layers. + cost_multiplier_layer : str, optional + Name of layer containing final cost layer spatial + multipliers. By default, ``None``. + cost_multiplier_scalar : int | float, optional + Final cost layer multiplier. By default, ``1``. + """ + + logger.debug( + # TODO: Turn this input a repr of routing scenario and use + # that here + "Building cost layer with the following inputs:" + "\n\t- cost_layers: %r" + "\n\t- friction_layers: %r" + "\n\t- cost_multiplier_layer: %r" + "\n\t- cost_multiplier_scalar: %r", + self.routing_scenario.cost_layers, + self.routing_scenario.friction_layers, + self.routing_scenario.cost_multiplier_layer, + self.routing_scenario.cost_multiplier_scalar, + ) + + self._build_cost_layer() + self._build_mcp_cost() + self._build_tracked_layers() + + return self + + def _build_cost_layer(self): + """Build out the main cost layer""" + self.cost = da.zeros(self._full_shape, dtype=np.float32) + self.li_cost = da.zeros(self._full_shape, dtype=np.float32) + for layer_info in self.routing_scenario.cost_layers: + layer_name = layer_info["layer_name"] + is_li = layer_info.get("is_invariant", False) + cost = self._extract_and_scale_layer(layer_info) + + if is_li: + self.li_cost += cost + self.li_cost_layer_map[layer_name] = cost + else: + self.cost += cost + + if layer_info.get("include_in_report", True): + self.tracked_layers.append( + CharacterizedLayer( + layer_name, cost, is_length_invariant=is_li + ) + ) + + if mll := self.routing_scenario.cost_multiplier_layer: + self._verify_layer_exists(mll) + self.cost *= self._layer_fh[mll].isel(band=0).astype(np.float32) + + self.cost *= self.routing_scenario.cost_multiplier_scalar + + def _build_mcp_cost(self): + """Build out the routing array""" + self.mcp_cost = self.cost.copy() + self.mcp_cost += self.li_cost + + # for li_cost_layer in self.li_cost_layer_map.values(): + # # normalize by cell size for routing only + # self.mcp_cost += li_cost_layer / self.cell_size + + friction_costs = da.zeros(self._full_shape, dtype=np.float32) + for layer_info in self.routing_scenario.friction_layers: + layer_name = layer_info["layer_name"] + friction_layer = self._extract_and_scale_layer( + layer_info, allow_cl=True + ) + # self.friction_layer_map[layer_name] = friction_layer + if layer_info.get("include_in_report", False): + self.tracked_layers.append( + CharacterizedLayer(layer_name, friction_layer) + ) + + friction_costs += friction_layer + + # Must happen at end of loop so that "lcp_agg_cost" + # remains constant + self.mcp_cost += friction_costs + + max_val = np.max(self.mcp_cost) * self.SOFT_BARRIER_MULTIPLIER + self.mcp_cost = da.where( + self.mcp_cost <= 0, + -1 if self.routing_scenario.use_hard_barrier else max_val, + self.mcp_cost, + ) + # logger.debug( + # "MCP cost min: %.2f, max: %.2f, median: %.2f", + # np.min(self.mcp_cost), + # np.max(self.mcp_cost), + # np.median(self.mcp_cost), + # ) + + def _extract_and_scale_layer(self, layer_info, allow_cl=False): + """Extract layer based on name and scale according to input""" + cost = self._extract_layer(layer_info["layer_name"], allow_cl=allow_cl) + + multiplier_layer_name = layer_info.get("multiplier_layer") + if multiplier_layer_name: + cost *= self._extract_layer( + multiplier_layer_name, allow_cl=allow_cl + ) + + cost *= layer_info.get("multiplier_scalar", 1) + return cost + + def _extract_layer(self, layer_name, allow_cl=False): + """Extract layer based on name""" + if allow_cl and layer_name == LCP_AGG_COST_LAYER_NAME: + return self.mcp_cost.copy() + + self._verify_layer_exists(layer_name) + return self._layer_fh[layer_name].isel(band=0).astype(np.float32) + + def _build_tracked_layers(self): + """Build out a dictionary of tracked layers""" + for ( + tracked_layer, + method, + ) in self.routing_scenario.tracked_layers.items(): + if getattr(da, method, None) is None: + msg = ( + f"Did not find method {method!r} in dask.array! " + f"Skipping tracked layer {tracked_layer!r}" + ) + warn(msg, revrtWarning) + continue + + if tracked_layer not in self._layers: + msg = ( + f"Did not find layer {tracked_layer!r} in cost " + f"file {str(self.routing_scenario.cost_fpath)!r}. Skipping..." + ) + warn(msg, revrtWarning) + continue + + layer = ( + self._layer_fh[tracked_layer].isel(band=0).astype(np.float32) + ) + self.tracked_layers.append( + CharacterizedLayer(tracked_layer, layer, agg_method=method) + ) + + +class CharacterizedLayer: + def __init__( + self, name, layer, is_length_invariant=False, agg_method=None + ): + self.name = name + self.layer = layer + self.is_length_invariant = is_length_invariant + self.agg_method = agg_method + + def __repr__(self): + return ( + f"CharacterizedLayer(name={self.name!r}, " + f"is_length_invariant={self.is_length_invariant}, " + f"agg_method={self.agg_method!r})" + ) + + def compute(self, route, cell_size): + """Compute aggregate value over route""" + rows, cols = np.array(route).T + layer_values = self.layer.isel( + y=xr.DataArray(rows, dims="points"), + x=xr.DataArray(cols, dims="points"), + ) + + if self.agg_method is None: + return self._compute_total_amd_length( + layer_values, route, cell_size + ) + + return self._compute_agg(layer_values) + + def _compute_total_amd_length(self, layer_values, route, cell_size): + # Use Pythagorean theorem to calculate length between cells (km) + # Use c**2 = a**2 + b**2 to determine length of individual paths + lens = np.sqrt(np.sum(np.diff(route, axis=0) ** 2, axis=1)) + + # Need to determine distance coming into and out of any cell. + # Assume paths start and end at the center of a cell. Therefore, + # distance traveled in the cell is half the distance entering it + # and half the distance exiting it. Duplicate all lengths, + # pad 0s on ends for start and end cells, and divide all + # distance by half. + lens = np.repeat(lens, 2) + lens = np.insert(np.append(lens, 0), 0, 0) + lens /= 2 + + # Group entrance and exits distance together, and add them + lens = lens.reshape((int(lens.shape[0] / 2), 2)) + lens = da.sum(lens, axis=1) + + if self.is_length_invariant: + layer_cost = da.sum(layer_values) + else: + layer_cost = da.sum(layer_values * lens) + + # Get path length in km only where layer costs are > 0 + layer_length = da.sum(lens[layer_values > 0]) * cell_size / 1000 + + return { + f"{self.name}_cost": layer_cost.astype(np.float32).compute(), + f"{self.name}_dist_km": ( + layer_length.astype(np.float32).compute() + ), + } + + def _compute_agg(self, layer_values): + aggregate = getattr(da, self.agg_method)(layer_values).astype( + np.float32 + ) + return {f"{self.name}_{self.agg_method}": aggregate.compute()} + + +class RouteResult: + """Class to compute route characteristics given layer cost maps""" + + def __init__( + self, routing_layers, route, optimized_objective, add_geom=False + ): + self._routing_layers = routing_layers + self._route = route + self._optimized_objective = optimized_objective + self._lens = self._total_path_length = None + self._by_layer_results = {} + self._add_geom = add_geom + + @property + def cell_size(self): + return abs(self._routing_layers.transform.a) + # return self._routing_layers.cell_size + + @property + def lens(self): + if self._lens is None: + self._compute_path_length() + return self._lens + + @property + def total_path_length(self): + if self._total_path_length is None: + self._compute_path_length() + return self._total_path_length + + @property + def cost(self): + rows, cols = np.array(self._route).T + cell_costs = self._routing_layers.cost.isel( + y=xr.DataArray(rows, dims="points"), + x=xr.DataArray(cols, dims="points"), + ) + + # Multiple distance travel through cell by cost and sum it! + return da.sum(cell_costs * self.lens).astype(np.float32).compute() + + @property + def end_lat(self): + row, col = self._route[-1] + return ( + self._routing_layers.latitudes.isel(y=row, x=col) + .astype(np.float32) + .compute() + .item() + ) + + @property + def end_lon(self): + row, col = self._route[-1] + return ( + self._routing_layers.longitudes.isel(y=row, x=col) + .astype(np.float32) + .compute() + .item() + ) + + @property + def geom(self): + rows, cols = np.array(self._route).T + x, y = rasterio.transform.xy( + self._routing_layers.transform, rows, cols + ) + geom = Point if len(self._route) == 1 else LineString + return geom(list(zip(x, y, strict=True))) + + def build(self): + results = { + "length_km": self.total_path_length, + "cost": self.cost, + "poi_lat": self.end_lat, + "poi_lon": self.end_lon, + "end_row": self._route[-1][0], + "end_col": self._route[-1][1], + "optimized_objective": self._optimized_objective, + } + for layer in self._routing_layers.tracked_layers: + layer_result = layer.compute(self._route, self.cell_size) + results.update(layer_result) + + if self._add_geom: + results["geometry"] = self.geom + + return results + + def _compute_path_length(self): + """Compute the total length and cell by cell length of LCP""" + # Use Pythagorean theorem to calculate length between cells (km) + # Use c**2 = a**2 + b**2 to determine length of individual paths + lens = np.sqrt(np.sum(np.diff(self._route, axis=0) ** 2, axis=1)) + self._total_path_length = np.sum(lens) * self.cell_size / 1000 + + # Need to determine distance coming into and out of any cell. + # Assume paths start and end at the center of a cell. Therefore, + # distance traveled in the cell is half the distance entering it + # and half the distance exiting it. Duplicate all lengths, + # pad 0s on ends for start and end cells, and divide all + # distance by half. + lens = np.repeat(lens, 2) + lens = np.insert(np.append(lens, 0), 0, 0) + lens /= 2 + + # Group entrance and exits distance together, and add them + lens = lens.reshape((int(lens.shape[0] / 2), 2)) + self._lens = np.sum(lens, axis=1) + + +def find_all_routes(routing_scenario, route_definitions, save_paths=False): + ts = time.monotonic() + + routing_layers = RoutingLayers(routing_scenario).build() + routes = _compute_routes( + routing_scenario, + route_definitions, + routing_layers, + save_paths=save_paths, + ) + + time_elapsed = f"{(time.monotonic() - ts) / 60:.4f} min" + logger.debug("Least Cost tie-line computed in %s", time_elapsed) + + return routes + + +def _compute_routes( + routing_scenario, route_definitions, routing_layers, save_paths +): + routes = [] + cost_layers = json.dumps({"cost_layers": routing_scenario.cost_layers}) + for start_point, end_points in route_definitions: + try: + indices, optimized_objective = find_paths( + zarr_fp=routing_scenario.cost_fpath, + cost_layers=cost_layers, + start=[start_point], + end=end_points, + )[0] + except Exception: + logger.error( + "Error computing least cost path from %s to %s", + start_point, + end_points, + ) + continue + + route = RouteResult( + routing_layers, + indices, + optimized_objective, + add_geom=save_paths, + ) + routes.append(route.build()) + + if save_paths: + if routes: + return gpd.GeoDataFrame( + routes, geometry="geometry", crs=routing_layers.cost_crs + ) + return gpd.GeoDataFrame([]) + + return pd.DataFrame(routes) From a3f4c3d8129e450b56bf25bf05b52f822119f6fd Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 21:42:51 -0700 Subject: [PATCH 11/61] More tests --- .../routing/test_routing_point_to_many.py | 67 ++++++++++++++++++- 1 file changed, 66 insertions(+), 1 deletion(-) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index 444945d0..0fc97ce4 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -160,7 +160,10 @@ def test_basic_single_route_layered_file(sample_layered_data): output = find_all_routes( scenario, - route_definitions=[((1, 1), [(2, 6)]), ((1, 2), [(2, 6)])], + route_definitions=[ + ((1, 1), [(2, 6)]), + ((1, 2), [(2, 6)]), + ], save_paths=False, ) @@ -178,5 +181,67 @@ def test_basic_single_route_layered_file(sample_layered_data): ) +def test_multi_layer_route_layered_file(sample_layered_data): + """Test routing across multiple cost layers""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[ + {"layer_name": "layer_1"}, + {"layer_name": "layer_2"}, + ], + ) + + output = find_all_routes( + scenario, + route_definitions=[((1, 1), [(2, 6)]), ((1, 2), [(2, 6)])], + save_paths=False, + ) + + assert len(output) == 2 + + first_route = output.iloc[0] + assert first_route["cost"] == pytest.approx( + 27.606602, + rel=1e-4, + ) + assert first_route["length_km"] == pytest.approx( + 0.005414, + rel=1e-4, + ) + assert first_route["layer_1_cost"] == pytest.approx( + 17.571068, + rel=1e-4, + ) + assert first_route["layer_2_cost"] == pytest.approx( + 10.035534, + rel=1e-4, + ) + assert np.isclose( + first_route["cost"], first_route["optimized_objective"], rtol=1e-6 + ) + + second_route = output.iloc[1] + assert second_route["cost"] == pytest.approx( + 25.106602, + rel=1e-4, + ) + assert second_route["length_km"] == pytest.approx( + 0.004414, + rel=1e-4, + ) + assert second_route["layer_1_cost"] == pytest.approx( + 16.071068, + rel=1e-4, + ) + assert second_route["layer_2_cost"] == pytest.approx( + 9.035534, + rel=1e-4, + ) + assert np.isclose( + second_route["cost"], second_route["optimized_objective"], rtol=1e-6 + ) + + if __name__ == "__main__": pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) From 953561bc2f5f472668d0624a73054257e7a0a440 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 22:08:27 -0700 Subject: [PATCH 12/61] Add test --- .../routing/test_routing_point_to_many.py | 72 +++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index 0fc97ce4..454eece1 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -243,5 +243,77 @@ def test_multi_layer_route_layered_file(sample_layered_data): ) +def test_multi_layer_route_with_multiplier(sample_layered_data): + """Test routing with multiple layers and a scalar multiplier""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[ + {"layer_name": "layer_1"}, + { + "layer_name": "layer_2", + "multiplier_scalar": 0.5, + }, + ], + ) + + output = find_all_routes( + scenario, + route_definitions=[ + ((1, 1), [(2, 6)]), + ((1, 2), [(2, 6)]), + ], + save_paths=False, + ) + + assert len(output) == 2 + + first_route = output.iloc[0] + assert first_route["cost"] == pytest.approx( + 22.588835, + rel=1e-4, + ) + assert first_route["length_km"] == pytest.approx( + 0.005414, + rel=1e-4, + ) + assert first_route["layer_1_cost"] == pytest.approx( + 17.571068, + rel=1e-4, + ) + assert first_route["layer_2_cost"] == pytest.approx( + 5.017767, + rel=1e-4, + ) + assert np.isclose( + first_route["cost"], + first_route["optimized_objective"], + rtol=1e-4, + ) + + second_route = output.iloc[1] + assert second_route["cost"] == pytest.approx( + 20.588835, + rel=1e-4, + ) + assert second_route["length_km"] == pytest.approx( + 0.004414, + rel=1e-4, + ) + assert second_route["layer_1_cost"] == pytest.approx( + 16.071068, + rel=1e-4, + ) + assert second_route["layer_2_cost"] == pytest.approx( + 4.517767, + rel=1e-4, + ) + assert np.isclose( + second_route["cost"], + second_route["optimized_objective"], + rtol=1e-4, + ) + + if __name__ == "__main__": pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) From 9213d5b257553af681afd048ff12b6d1278af40e Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 22:19:14 -0700 Subject: [PATCH 13/61] Add test --- .../routing/test_routing_point_to_many.py | 44 +++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index 454eece1..9d1b9f63 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -315,5 +315,49 @@ def test_multi_layer_route_with_multiplier(sample_layered_data): ) +def test_multi_layer_route_with_scalar_and_layer_multipliers( + sample_layered_data, +): + """Test routing when combining scalar and layer multipliers""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[ + {"layer_name": "layer_1"}, + {"layer_name": "layer_2", "multiplier_scalar": 0.5}, + { + "layer_name": "layer_3", + "multiplier_layer": "layer_4", + }, + { + "layer_name": "layer_5", + "multiplier_scalar": 2, + "multiplier_layer": "layer_4", + }, + ], + ) + + output = find_all_routes( + scenario, + route_definitions=[((1, 1), [(1, 2)])], + save_paths=False, + ) + + assert len(output) == 1 + + route = output.iloc[0] + assert route["cost"] == pytest.approx(2.0, rel=1e-4) + assert route["length_km"] == pytest.approx(0.001, rel=1e-4) + assert route["layer_1_cost"] == pytest.approx(1.5, rel=1e-4) + assert route["layer_2_cost"] == pytest.approx(0.5, rel=1e-4) + assert route["layer_3_cost"] == pytest.approx(0.0, abs=1e-8) + assert route["layer_5_cost"] == pytest.approx(0.0, abs=1e-8) + assert route["layer_1_dist_km"] == pytest.approx(0.001, rel=1e-4) + assert route["layer_2_dist_km"] == pytest.approx(0.001, rel=1e-4) + assert route["layer_3_dist_km"] == pytest.approx(0.001, rel=1e-4) + assert route["layer_5_dist_km"] == pytest.approx(0.001, rel=1e-4) + assert np.isclose(route["cost"], route["optimized_objective"], rtol=1e-6) + + if __name__ == "__main__": pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) From c1bc6450dd777d6178f2907640dabf850be043f8 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 22:54:16 -0700 Subject: [PATCH 14/61] Fix typo --- revrt/routing/point_to_many.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/revrt/routing/point_to_many.py b/revrt/routing/point_to_many.py index 9b41f851..904487fd 100644 --- a/revrt/routing/point_to_many.py +++ b/revrt/routing/point_to_many.py @@ -325,13 +325,13 @@ def compute(self, route, cell_size): ) if self.agg_method is None: - return self._compute_total_amd_length( + return self._compute_total_and_length( layer_values, route, cell_size ) return self._compute_agg(layer_values) - def _compute_total_amd_length(self, layer_values, route, cell_size): + def _compute_total_and_length(self, layer_values, route, cell_size): # Use Pythagorean theorem to calculate length between cells (km) # Use c**2 = a**2 + b**2 to determine length of individual paths lens = np.sqrt(np.sum(np.diff(route, axis=0) ** 2, axis=1)) From f483595d8fb89473e083fe24574a6f678054812b Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 22:59:54 -0700 Subject: [PATCH 15/61] Fix bug --- revrt/routing/point_to_many.py | 73 +++++++++---------- .../routing/test_routing_point_to_many.py | 4 +- 2 files changed, 36 insertions(+), 41 deletions(-) diff --git a/revrt/routing/point_to_many.py b/revrt/routing/point_to_many.py index 904487fd..5b7fccad 100644 --- a/revrt/routing/point_to_many.py +++ b/revrt/routing/point_to_many.py @@ -332,31 +332,18 @@ def compute(self, route, cell_size): return self._compute_agg(layer_values) def _compute_total_and_length(self, layer_values, route, cell_size): - # Use Pythagorean theorem to calculate length between cells (km) - # Use c**2 = a**2 + b**2 to determine length of individual paths - lens = np.sqrt(np.sum(np.diff(route, axis=0) ** 2, axis=1)) - - # Need to determine distance coming into and out of any cell. - # Assume paths start and end at the center of a cell. Therefore, - # distance traveled in the cell is half the distance entering it - # and half the distance exiting it. Duplicate all lengths, - # pad 0s on ends for start and end cells, and divide all - # distance by half. - lens = np.repeat(lens, 2) - lens = np.insert(np.append(lens, 0), 0, 0) - lens /= 2 - - # Group entrance and exits distance together, and add them - lens = lens.reshape((int(lens.shape[0] / 2), 2)) - lens = da.sum(lens, axis=1) + lens, __ = _compute_lens(route, cell_size) + + layer_data = getattr(layer_values, "data", layer_values) + if not isinstance(layer_data, da.Array): + layer_data = da.asarray(layer_data) if self.is_length_invariant: - layer_cost = da.sum(layer_values) + layer_cost = da.sum(layer_data) else: - layer_cost = da.sum(layer_values * lens) + layer_cost = da.sum(layer_data * lens) - # Get path length in km only where layer costs are > 0 - layer_length = da.sum(lens[layer_values > 0]) * cell_size / 1000 + layer_length = da.sum(lens[layer_data > 0]) * cell_size / 1000 return { f"{self.name}_cost": layer_cost.astype(np.float32).compute(), @@ -463,24 +450,9 @@ def build(self): def _compute_path_length(self): """Compute the total length and cell by cell length of LCP""" - # Use Pythagorean theorem to calculate length between cells (km) - # Use c**2 = a**2 + b**2 to determine length of individual paths - lens = np.sqrt(np.sum(np.diff(self._route, axis=0) ** 2, axis=1)) - self._total_path_length = np.sum(lens) * self.cell_size / 1000 - - # Need to determine distance coming into and out of any cell. - # Assume paths start and end at the center of a cell. Therefore, - # distance traveled in the cell is half the distance entering it - # and half the distance exiting it. Duplicate all lengths, - # pad 0s on ends for start and end cells, and divide all - # distance by half. - lens = np.repeat(lens, 2) - lens = np.insert(np.append(lens, 0), 0, 0) - lens /= 2 - - # Group entrance and exits distance together, and add them - lens = lens.reshape((int(lens.shape[0] / 2), 2)) - self._lens = np.sum(lens, axis=1) + self._lens, self._total_path_length = _compute_lens( + self._route, self.cell_size + ) def find_all_routes(routing_scenario, route_definitions, save_paths=False): @@ -537,3 +509,26 @@ def _compute_routes( return gpd.GeoDataFrame([]) return pd.DataFrame(routes) + + +def _compute_lens(route, cell_size): + """Compute the total length and cell by cell length of LCP""" + # Use Pythagorean theorem to calculate length between cells (km) + # Use c**2 = a**2 + b**2 to determine length of individual paths + lens = np.sqrt(np.sum(np.diff(route, axis=0) ** 2, axis=1)) + total_path_length = np.sum(lens) * cell_size / 1000 + + # Need to determine distance coming into and out of any cell. + # Assume paths start and end at the center of a cell. Therefore, + # distance traveled in the cell is half the distance entering it + # and half the distance exiting it. Duplicate all lengths, + # pad 0s on ends for start and end cells, and divide all + # distance by half. + lens = np.repeat(lens, 2) + lens = np.insert(np.append(lens, 0), 0, 0) + lens /= 2 + + # Group entrance and exits distance together, and add them + lens = lens.reshape((int(lens.shape[0] / 2), 2)) + lens = np.sum(lens, axis=1) + return lens, total_path_length diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index 9d1b9f63..155170ce 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -354,8 +354,8 @@ def test_multi_layer_route_with_scalar_and_layer_multipliers( assert route["layer_5_cost"] == pytest.approx(0.0, abs=1e-8) assert route["layer_1_dist_km"] == pytest.approx(0.001, rel=1e-4) assert route["layer_2_dist_km"] == pytest.approx(0.001, rel=1e-4) - assert route["layer_3_dist_km"] == pytest.approx(0.001, rel=1e-4) - assert route["layer_5_dist_km"] == pytest.approx(0.001, rel=1e-4) + assert route["layer_3_dist_km"] == pytest.approx(0.0, abs=1e-8) + assert route["layer_5_dist_km"] == pytest.approx(0.0, abs=1e-8) assert np.isclose(route["cost"], route["optimized_objective"], rtol=1e-6) From 08b14b31ae00846d426872529f17396c19d4357b Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 23:02:20 -0700 Subject: [PATCH 16/61] Add tracked layers test --- .../routing/test_routing_point_to_many.py | 33 +++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index 155170ce..e260480e 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -359,5 +359,38 @@ def test_multi_layer_route_with_scalar_and_layer_multipliers( assert np.isclose(route["cost"], route["optimized_objective"], rtol=1e-6) +def test_routing_with_tracked_layers(sample_layered_data): + """Tracked layers report aggregated stats alongside routing results""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + tracked_layers={ + "layer_1": "mean", + "layer_2": "max", + "layer_3": "min", + }, + ) + + output = find_all_routes( + scenario, + route_definitions=[((1, 1), [(1, 2)])], + save_paths=False, + ) + + assert len(output) == 1 + route = output.iloc[0] + + assert { + "layer_1_mean", + "layer_2_max", + "layer_3_min", + }.issubset(output.columns) + + assert route["layer_1_mean"] == pytest.approx(1.5) + assert route["layer_2_max"] == pytest.approx(1.0) + assert route["layer_3_min"] == pytest.approx(2.0) + + if __name__ == "__main__": pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) From 5bffe8ac1d67f4b0340f1f01067554bb846e031b Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 23:06:27 -0700 Subject: [PATCH 17/61] Add geometry test --- .../routing/test_routing_point_to_many.py | 54 +++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index e260480e..da23be33 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -2,6 +2,7 @@ from pathlib import Path +import geopandas as gpd import pytest import numpy as np import xarray as xr @@ -243,6 +244,59 @@ def test_multi_layer_route_layered_file(sample_layered_data): ) +def test_save_paths_returns_expected_geometry(sample_layered_data): + """Saving paths returns expected geometries for each route""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + output = find_all_routes( + scenario, + route_definitions=[ + ((1, 1), [(2, 6)]), + ((1, 2), [(2, 6)]), + ], + save_paths=True, + ) + + assert isinstance(output, gpd.GeoDataFrame) + assert len(output) == 2 + assert output.crs.to_string() == "EPSG:4326" + + expected_geometries = [ + [ + (1.5, 5.5), + (1.5, 4.5), + (2.5, 3.5), + (3.5, 2.5), + (4.5, 1.5), + (5.5, 2.5), + (5.5, 3.5), + (6.5, 4.5), + ], + [ + (2.5, 5.5), + (2.5, 4.5), + (3.5, 3.5), + (3.5, 2.5), + (4.5, 1.5), + (5.5, 2.5), + (5.5, 3.5), + (6.5, 4.5), + ], + ] + + for geom, expected_coords in zip( + output.geometry, expected_geometries, strict=True + ): + assert geom.geom_type == "LineString" + assert np.allclose( + np.asarray(geom.coords), np.asarray(expected_coords) + ) + + def test_multi_layer_route_with_multiplier(sample_layered_data): """Test routing with multiple layers and a scalar multiplier""" From 7e5df95af85aee595cbc8233d40ff05d7546118f Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Sat, 6 Dec 2025 23:12:30 -0700 Subject: [PATCH 18/61] A few more tests --- .../routing/test_routing_point_to_many.py | 41 +++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index da23be33..877fb4a2 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -3,6 +3,7 @@ from pathlib import Path import geopandas as gpd +import pandas as pd import pytest import numpy as np import xarray as xr @@ -297,6 +298,46 @@ def test_save_paths_returns_expected_geometry(sample_layered_data): ) +def test_empty_route_definitions_returns_empty_dataframe(sample_layered_data): + """Empty route definitions return an empty dataframe""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + output = find_all_routes( + scenario, + route_definitions=[], + save_paths=False, + ) + + assert isinstance(output, pd.DataFrame) + assert output.empty + assert list(output.columns) == [] + + +def test_empty_route_definitions_returns_empty_geo_dataframe( + sample_layered_data, +): + """Empty route definitions return an empty dataframe""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + output = find_all_routes( + scenario, + route_definitions=[], + save_paths=True, + ) + + assert isinstance(output, gpd.GeoDataFrame) + assert output.empty + assert list(output.columns) == [] + + def test_multi_layer_route_with_multiplier(sample_layered_data): """Test routing with multiple layers and a scalar multiplier""" From fd86ec315f4e5bd61d2ecc6674c561949ee447b0 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Mon, 8 Dec 2025 18:37:23 -0700 Subject: [PATCH 19/61] Add `__init__.py` to module --- revrt/routing/__init__.py | 1 + 1 file changed, 1 insertion(+) create mode 100644 revrt/routing/__init__.py diff --git a/revrt/routing/__init__.py b/revrt/routing/__init__.py new file mode 100644 index 00000000..def6ca7a --- /dev/null +++ b/revrt/routing/__init__.py @@ -0,0 +1 @@ +"""reVRt routing logic""" From 97f49fb28987c208412a403b689bddf6c40996ff Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Mon, 8 Dec 2025 20:16:08 -0700 Subject: [PATCH 20/61] Updates --- revrt/routing/point_to_many.py | 111 +++++++++++++++++++++++---------- 1 file changed, 79 insertions(+), 32 deletions(-) diff --git a/revrt/routing/point_to_many.py b/revrt/routing/point_to_many.py index 5b7fccad..1d1aceca 100644 --- a/revrt/routing/point_to_many.py +++ b/revrt/routing/point_to_many.py @@ -4,6 +4,7 @@ import time import logging from warnings import warn +from functools import cached_property import rasterio import numpy as np @@ -16,11 +17,11 @@ from revrt import find_paths -# from revrt.exceptions import ( -# # revrtValueError, -# # revrtInvalidStartCostError, -# revrtLeastCostPathNotFoundError, -# ) +from revrt.exceptions import ( + # revrtValueError, + # revrtInvalidStartCostError, + revrtLeastCostPathNotFoundError, +) from revrt.warn import revrtWarning logger = logging.getLogger(__name__) @@ -47,6 +48,10 @@ def __init__( self.cost_multiplier_scalar = cost_multiplier_scalar self.use_hard_barrier = use_hard_barrier + @cached_property + def cl_as_json(self): + return json.dumps({"cost_layers": self.cost_layers}) + class RoutingLayers: """Class to build a routing layer from user input""" @@ -99,6 +104,7 @@ def __init__(self, routing_scenario, chunks="auto"): self.cost = None self.li_cost = None + self.final_routing_layer = None @property def latitudes(self): @@ -176,7 +182,7 @@ def build(self): ) self._build_cost_layer() - self._build_mcp_cost() + self._build_final_routing_layer() self._build_tracked_layers() return self @@ -208,15 +214,16 @@ def _build_cost_layer(self): self.cost *= self._layer_fh[mll].isel(band=0).astype(np.float32) self.cost *= self.routing_scenario.cost_multiplier_scalar + self.li_cost += self.cost * 0 - def _build_mcp_cost(self): + def _build_final_routing_layer(self): """Build out the routing array""" - self.mcp_cost = self.cost.copy() - self.mcp_cost += self.li_cost + self.final_routing_layer = self.cost.copy() + self.final_routing_layer += self.li_cost # for li_cost_layer in self.li_cost_layer_map.values(): # # normalize by cell size for routing only - # self.mcp_cost += li_cost_layer / self.cell_size + # self.final_routing_layer += li_cost_layer / self.cell_size friction_costs = da.zeros(self._full_shape, dtype=np.float32) for layer_info in self.routing_scenario.friction_layers: @@ -234,19 +241,21 @@ def _build_mcp_cost(self): # Must happen at end of loop so that "lcp_agg_cost" # remains constant - self.mcp_cost += friction_costs + self.final_routing_layer += friction_costs - max_val = np.max(self.mcp_cost) * self.SOFT_BARRIER_MULTIPLIER - self.mcp_cost = da.where( - self.mcp_cost <= 0, - -1 if self.routing_scenario.use_hard_barrier else max_val, - self.mcp_cost, + max_val = ( + np.max(self.final_routing_layer) * self.SOFT_BARRIER_MULTIPLIER ) + self.final_routing_layer = da.where( + self.final_routing_layer <= 0, + -1 if self.routing_scenario.use_hard_barrier else max_val, + self.final_routing_layer, + ) + (self.cost * 0) # logger.debug( # "MCP cost min: %.2f, max: %.2f, median: %.2f", - # np.min(self.mcp_cost), - # np.max(self.mcp_cost), - # np.median(self.mcp_cost), + # np.min(self.final_routing_layer), + # np.max(self.final_routing_layer), + # np.median(self.final_routing_layer), # ) def _extract_and_scale_layer(self, layer_info, allow_cl=False): @@ -265,7 +274,7 @@ def _extract_and_scale_layer(self, layer_info, allow_cl=False): def _extract_layer(self, layer_name, allow_cl=False): """Extract layer based on name""" if allow_cl and layer_name == LCP_AGG_COST_LAYER_NAME: - return self.mcp_cost.copy() + return self.final_routing_layer.copy() self._verify_layer_exists(layer_name) return self._layer_fh[layer_name].isel(band=0).astype(np.float32) @@ -476,21 +485,12 @@ def _compute_routes( routing_scenario, route_definitions, routing_layers, save_paths ): routes = [] - cost_layers = json.dumps({"cost_layers": routing_scenario.cost_layers}) for start_point, end_points in route_definitions: try: - indices, optimized_objective = find_paths( - zarr_fp=routing_scenario.cost_fpath, - cost_layers=cost_layers, - start=[start_point], - end=end_points, - )[0] - except Exception: - logger.error( - "Error computing least cost path from %s to %s", - start_point, - end_points, + indices, optimized_objective = _compute_valid_path( + routing_scenario, routing_layers, start_point, end_points ) + except revrtLeastCostPathNotFoundError: continue route = RouteResult( @@ -511,6 +511,53 @@ def _compute_routes( return pd.DataFrame(routes) +def _compute_valid_path( + routing_scenario, routing_layers, start_point, end_points +): + start_row, start_col = start_point + # print(f"{routing_layers.final_routing_layer=}") + start_cost = ( + routing_layers.final_routing_layer.isel(y=start_row, x=start_col) + .compute() + .item() + ) + + if start_cost <= 0: + msg = ( + f"Start idx {start_point} does not have a valid cost: " + f"{start_cost:.2f} (must be > 0)!" + ) + raise revrtLeastCostPathNotFoundError(msg) + + end_rows, end_cols = np.array(end_points).T + end_costs = routing_layers.final_routing_layer.isel( + y=xr.DataArray(end_rows, dims="points"), + x=xr.DataArray(end_cols, dims="points"), + ) + if not np.any(end_costs.compute() > 0): + msg = ( + f"None of the end idx {end_points} have a valid cost " + f"(must be > 0)!" + ) + raise revrtLeastCostPathNotFoundError(msg) + + try: + indices, optimized_objective = find_paths( + zarr_fp=routing_scenario.cost_fpath, + cost_layers=routing_scenario.cl_as_json, + start=[start_point], + end=end_points, + )[0] + except Exception as ex: + msg = ( + f"Unable to find path from {start_point} any of {end_points}: {ex}" + ) + + raise revrtLeastCostPathNotFoundError(msg) from ex + + return indices, optimized_objective + + def _compute_lens(route, cell_size): """Compute the total length and cell by cell length of LCP""" # Use Pythagorean theorem to calculate length between cells (km) From 4e828653bd6c1ddcd408445a92afe33c84aac4f2 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Mon, 8 Dec 2025 20:58:40 -0700 Subject: [PATCH 21/61] Don't allow negative costs --- crates/revrt/src/dataset.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/crates/revrt/src/dataset.rs b/crates/revrt/src/dataset.rs index 4281f255..e55e6c2c 100644 --- a/crates/revrt/src/dataset.rs +++ b/crates/revrt/src/dataset.rs @@ -292,6 +292,7 @@ impl Dataset { let neighbors = neighbors .iter() .filter(|((ir, jr), _)| !(*ir == i && *jr == j)) // no center point + .filter(|((_, _), v)| *v > 0.) // only positive costs .map(|((ir, jr), v)| ((ir, jr), 0.5 * (v + center.1))) .map(|((ir, jr), v)| { if *ir != i && *jr != j { From 9b31040df0f7549a75499720ba2f490edd9bc415 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Mon, 8 Dec 2025 20:59:49 -0700 Subject: [PATCH 22/61] Add tests for invalid start/end costs --- .../routing/test_routing_point_to_many.py | 61 ++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index 877fb4a2..e5ab533f 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -35,7 +35,7 @@ def sample_layered_data(tmp_path_factory): layer_1 = np.array( [ [ - [7, 7, 8, 0, 9, 9, 9, 0], + [7, 7, 8, -1, 9, 9, 9, 0], [8, 1, 2, 2, 9, 9, 9, 0], [9, 1, 3, 3, 9, 1, 2, 3], [9, 1, 2, 1, 9, 1, 9, 0], @@ -487,5 +487,64 @@ def test_routing_with_tracked_layers(sample_layered_data): assert route["layer_3_min"] == pytest.approx(2.0) +def test_start_point_on_barrier_returns_no_route(sample_layered_data): + """If the start point is on a barrier (cost <= 0) no route is returned""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + # (0, 3) in layer_1 is 0 -> treated as barrier + output = find_all_routes( + scenario, route_definitions=[((0, 3), [(2, 6)])], save_paths=False + ) + + assert isinstance(output, pd.DataFrame) + assert output.empty + + +def test_some_endpoints_include_barriers_but_one_valid(sample_layered_data): + """If some end points <=0 but at least one is valid, route is found""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + # include one barrier end (0,3) and one valid end (2,6) + output = find_all_routes( + scenario, + route_definitions=[((1, 1), [(0, 3), (2, 6)])], + save_paths=False, + ) + + assert len(output) == 1 + # At least one valid endpoint must be reached and cost must be positive. + assert output.iloc[0]["cost"] > 0 + + end_row = int(output.iloc[0]["end_row"]) + end_col = int(output.iloc[0]["end_col"]) + assert (end_row, end_col) == (2, 6) + + +def test_all_endpoints_are_barriers_returns_no_route(sample_layered_data): + """If all end points are barriers, no route is returned""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + output = find_all_routes( + scenario, + route_definitions=[((1, 1), [(0, 3), (0, 7)])], + save_paths=False, + ) + + assert isinstance(output, pd.DataFrame) + assert output.empty + + if __name__ == "__main__": pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) From dce81b3c48b0d9d3315fbb3df06e9ae91232be67 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Mon, 8 Dec 2025 21:03:52 -0700 Subject: [PATCH 23/61] Check for error message logs --- .../unit/routing/test_routing_point_to_many.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index e5ab533f..adc88107 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -35,7 +35,7 @@ def sample_layered_data(tmp_path_factory): layer_1 = np.array( [ [ - [7, 7, 8, -1, 9, 9, 9, 0], + [7, 7, 8, 0, 9, 9, 9, 0], [8, 1, 2, 2, 9, 9, 9, 0], [9, 1, 3, 3, 9, 1, 2, 3], [9, 1, 2, 1, 9, 1, 9, 0], @@ -487,7 +487,9 @@ def test_routing_with_tracked_layers(sample_layered_data): assert route["layer_3_min"] == pytest.approx(2.0) -def test_start_point_on_barrier_returns_no_route(sample_layered_data): +def test_start_point_on_barrier_returns_no_route( + sample_layered_data, assert_message_was_logged +): """If the start point is on a barrier (cost <= 0) no route is returned""" scenario = RoutingScenario( @@ -499,6 +501,10 @@ def test_start_point_on_barrier_returns_no_route(sample_layered_data): output = find_all_routes( scenario, route_definitions=[((0, 3), [(2, 6)])], save_paths=False ) + assert_message_was_logged( + "Start idx (0, 3) does not have a valid cost: -1.00 (must be > 0)!", + "ERROR", + ) assert isinstance(output, pd.DataFrame) assert output.empty @@ -528,7 +534,9 @@ def test_some_endpoints_include_barriers_but_one_valid(sample_layered_data): assert (end_row, end_col) == (2, 6) -def test_all_endpoints_are_barriers_returns_no_route(sample_layered_data): +def test_all_endpoints_are_barriers_returns_no_route( + sample_layered_data, assert_message_was_logged +): """If all end points are barriers, no route is returned""" scenario = RoutingScenario( @@ -541,6 +549,10 @@ def test_all_endpoints_are_barriers_returns_no_route(sample_layered_data): route_definitions=[((1, 1), [(0, 3), (0, 7)])], save_paths=False, ) + assert_message_was_logged( + "None of the end idx [(0, 3), (0, 7)] have a valid cost (must be > 0)", + "ERROR", + ) assert isinstance(output, pd.DataFrame) assert output.empty From fe8ee72c380a90a1fd8439488c488538b5c6bff8 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Mon, 8 Dec 2025 22:09:42 -0700 Subject: [PATCH 24/61] Add routing utility tests --- .../test_routing_utilities_integration.py | 101 ++++++++ .../unit/routing/test_routing_utilities.py | 215 ++++++++++++++++++ 2 files changed, 316 insertions(+) create mode 100644 tests/python/integration/test_routing_utilities_integration.py create mode 100644 tests/python/unit/routing/test_routing_utilities.py diff --git a/tests/python/integration/test_routing_utilities_integration.py b/tests/python/integration/test_routing_utilities_integration.py new file mode 100644 index 00000000..c25cad46 --- /dev/null +++ b/tests/python/integration/test_routing_utilities_integration.py @@ -0,0 +1,101 @@ +"""reVrt tests for routing utilities""" + +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from rasterio.transform import from_origin + +from revrt.routing.utilities import map_to_costs +from revrt.utilities import LayeredFile + + +@pytest.fixture(scope="module") +def sample_layered_data(tmp_path_factory): + """Sample layered data files to use across tests""" + data_dir = tmp_path_factory.mktemp("routing_data") + + layered_fp = data_dir / "test_layered.zarr" + layer_file = LayeredFile(layered_fp) + + height, width = (7, 8) + cell_size = 1.0 + x0, y0 = 0.0, float(height) + transform = from_origin(x0, y0, cell_size, cell_size) + x_coords = ( + x0 + np.arange(width, dtype=np.float32) * cell_size + cell_size / 2 + ) + y_coords = ( + y0 - np.arange(height, dtype=np.float32) * cell_size - cell_size / 2 + ) + + layer_1 = np.array( + [ + [ + [7, 7, 8, 0, 9, 9, 9, 0], + [8, 1, 2, 2, 9, 9, 9, 0], + [9, 1, 3, 3, 9, 1, 2, 3], + [9, 1, 2, 1, 9, 1, 9, 0], + [9, 9, 9, 1, 9, 1, 9, 0], + [9, 9, 9, 1, 1, 1, 9, 0], + [9, 9, 9, 9, 9, 9, 9, 0], + ] + ], + dtype=np.float32, + ) + + da = xr.DataArray( + layer_1, + dims=("band", "y", "x"), + coords={"y": y_coords, "x": x_coords}, + ) + da = da.rio.write_crs("EPSG:4326") + da = da.rio.write_transform(transform) + + geotiff_fp = data_dir / "layer_1.tif" + da.rio.to_raster(geotiff_fp, driver="GTiff") + + layer_file.write_geotiff_to_file(geotiff_fp, "layer_1", overwrite=True) + return layered_fp + + +def test_basic_map_to_costs_integration(sample_layered_data): + """Basic integration test for mapping to costs from layered data""" + + route_table = pd.DataFrame( + { + "start_lat": 5.5, + "start_lon": [1.5, 2.5], + "end_lat": 4.5, + "end_lon": 6.5, + } + ) + + with xr.open_dataset( + sample_layered_data, consolidated=False, engine="zarr" + ) as ds: + route_table = map_to_costs( + route_table, + crs=ds.rio.crs, + transform=ds.rio.transform(), + shape=ds.rio.shape, + ) + + np.testing.assert_array_equal( + route_table["start_row"].to_numpy(), np.array([1, 1]) + ) + np.testing.assert_array_equal( + route_table["start_col"].to_numpy(), np.array([1, 2]) + ) + np.testing.assert_array_equal( + route_table["end_row"].to_numpy(), np.array([2, 2]) + ) + np.testing.assert_array_equal( + route_table["end_col"].to_numpy(), np.array([6, 6]) + ) + + +if __name__ == "__main__": + pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) diff --git a/tests/python/unit/routing/test_routing_utilities.py b/tests/python/unit/routing/test_routing_utilities.py new file mode 100644 index 00000000..20123157 --- /dev/null +++ b/tests/python/unit/routing/test_routing_utilities.py @@ -0,0 +1,215 @@ +"""reVrt tests for routing utilities""" + +from pathlib import Path +import warnings + +import numpy as np +import pandas as pd +import pytest +from rasterio.transform import from_origin, xy + +from revrt.routing import utilities as routing_utils +from revrt.warn import revrtWarning + + +@pytest.fixture +def cost_grid(): + """Simple grid metadata for routing tests""" + + height, width = (4, 5) + cell_size = 1.0 + transform = from_origin(0.0, float(height), cell_size, cell_size) + return "EPSG:4326", transform, (height, width) + + +def test_transform_lat_lon_to_row_col_expected_indices(cost_grid): + """Map lat/lon pairs to expected raster indices""" + + crs, transform, _ = cost_grid + lon_a, lat_a = xy(transform, 0, 0, offset="center") + lon_b, lat_b = xy(transform, 3, 4, offset="center") + + row, col = routing_utils._transform_lat_lon_to_row_col( + transform, crs, np.array([lat_a, lat_b]), np.array([lon_a, lon_b]) + ) + + assert isinstance(row, np.ndarray) + assert isinstance(col, np.ndarray) + np.testing.assert_array_equal(row, np.array([0, 3])) + np.testing.assert_array_equal(col, np.array([0, 4])) + + +def test_get_start_end_point_cost_indices_adds_expected_columns(cost_grid): + """Populate start/end index columns from coordinates""" + + crs, transform, _ = cost_grid + lon_start_a, lat_start_a = xy(transform, 0, 0, offset="center") + lon_start_b, lat_start_b = xy(transform, 1, 2, offset="center") + lon_end_a, lat_end_a = xy(transform, 2, 3, offset="center") + lon_end_b, lat_end_b = xy(transform, 3, 4, offset="center") + + route_points = pd.DataFrame( + { + "start_lat": [str(lat_start_a), str(lat_start_b)], + "start_lon": [str(lon_start_a), str(lon_start_b)], + "end_lat": [str(lat_end_a), str(lat_end_b)], + "end_lon": [str(lon_end_a), str(lon_end_b)], + } + ) + + updated = routing_utils._get_start_end_point_cost_indices( + route_points, crs, transform + ) + + assert updated is route_points + np.testing.assert_array_equal( + route_points["start_row"].to_numpy(), np.array([0, 1]) + ) + np.testing.assert_array_equal( + route_points["start_col"].to_numpy(), np.array([0, 2]) + ) + np.testing.assert_array_equal( + route_points["end_row"].to_numpy(), np.array([2, 3]) + ) + np.testing.assert_array_equal( + route_points["end_col"].to_numpy(), np.array([3, 4]) + ) + + +def test_filter_points_outside_cost_domain_no_warning(cost_grid): + """Do not warn when all routes remain in bounds""" + + _, _, shape = cost_grid + route_points = pd.DataFrame( + { + "start_row": [0, 1], + "start_col": [0, 2], + "end_row": [2, 3], + "end_col": [3, 4], + } + ) + expected = route_points.copy(deep=True) + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + filtered = routing_utils._filter_points_outside_cost_domain( + route_points, shape + ) + + assert not caught + pd.testing.assert_frame_equal(filtered, expected) + + +def test_filter_points_outside_cost_domain_warns_and_drops(cost_grid): + """Warn and drop routes that fall outside bounds""" + + _, _, shape = cost_grid + route_points = pd.DataFrame( + { + "start_row": [0, -1], + "start_col": [0, 1], + "end_row": [1, 5], + "end_col": [1, 6], + "label": ["valid", "invalid"], + } + ) + + with pytest.warns(revrtWarning) as record: + filtered = routing_utils._filter_points_outside_cost_domain( + route_points, shape + ) + + assert len(record) == 1 + assert "outside of the cost exclusion domain" in str(record[0].message) + assert filtered.index.tolist() == [0] + assert filtered["label"].tolist() == ["valid"] + + +def test_map_to_costs_maps_and_preserves_valid_routes(cost_grid): + """End-to-end mapping keeps valid route intact""" + + crs, transform, shape = cost_grid + lon_start, lat_start = xy(transform, 0, 0, offset="center") + lon_end, lat_end = xy(transform, 2, 3, offset="center") + + route_points = pd.DataFrame( + { + "start_lat": [lat_start], + "start_lon": [lon_start], + "end_lat": [lat_end], + "end_lon": [lon_end], + } + ) + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + mapped = routing_utils.map_to_costs( + route_points.copy(deep=True), crs, transform, shape + ) + + revrt_warnings = [w for w in caught if isinstance(w.message, revrtWarning)] + assert not revrt_warnings + + np.testing.assert_array_equal( + mapped["start_row"].to_numpy(), np.array([0]) + ) + np.testing.assert_array_equal( + mapped["start_col"].to_numpy(), np.array([0]) + ) + np.testing.assert_array_equal( + mapped["end_row"].to_numpy(), + np.array([2]), + ) + np.testing.assert_array_equal( + mapped["end_col"].to_numpy(), + np.array([3]), + ) + + +def test_map_to_costs_filters_routes_outside_cost_domain(cost_grid): + """Remove routes that leave the cost domain""" + + crs, transform, shape = cost_grid + lon_valid_start, lat_valid_start = xy(transform, 0, 0, offset="center") + lon_valid_end, lat_valid_end = xy(transform, 1, 1, offset="center") + + # Column beyond grid width keeps latitude inside range but forces filtering + lon_outside = lon_valid_start + shape[1] + lat_outside = lat_valid_start + + route_points = pd.DataFrame( + { + "start_lat": [lat_valid_start, lat_outside], + "start_lon": [lon_valid_start, lon_outside], + "end_lat": [lat_valid_end, lat_outside], + "end_lon": [lon_valid_end, lon_outside], + } + ) + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + mapped = routing_utils.map_to_costs( + route_points, crs, transform, shape + ) + + revrt_warnings = [w for w in caught if isinstance(w.message, revrtWarning)] + assert len(revrt_warnings) == 1 + np.testing.assert_array_equal( + mapped["start_row"].to_numpy(), np.array([0]) + ) + np.testing.assert_array_equal( + mapped["start_col"].to_numpy(), np.array([0]) + ) + np.testing.assert_array_equal( + mapped["end_row"].to_numpy(), + np.array([1]), + ) + np.testing.assert_array_equal( + mapped["end_col"].to_numpy(), + np.array([1]), + ) + assert mapped.index.tolist() == [0] + + +if __name__ == "__main__": + pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) From a80c677126deae67f2cf08f0e3fc5eb670ec5071 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Mon, 8 Dec 2025 22:16:53 -0700 Subject: [PATCH 25/61] Add utilities module --- revrt/routing/utilities.py | 107 +++++++++++++++++++++++++++++++++++++ 1 file changed, 107 insertions(+) create mode 100644 revrt/routing/utilities.py diff --git a/revrt/routing/utilities.py b/revrt/routing/utilities.py new file mode 100644 index 00000000..53738551 --- /dev/null +++ b/revrt/routing/utilities.py @@ -0,0 +1,107 @@ +"""reVRt routing module utilities""" + +import logging +from warnings import warn + +import rasterio +import numpy as np +import geopandas as gpd +from shapely.geometry import Point + +from revrt.warn import revrtWarning + + +logger = logging.getLogger(__name__) + + +def map_to_costs(route_points, crs, transform, shape): + """_summary_ + + Parameters + ---------- + route_points : _type_ + _description_ + crs : _type_ + _description_ + transform : _type_ + _description_ + shape : _type_ + _description_ + + Returns + ------- + _type_ + _description_ + """ + route_points = _get_start_end_point_cost_indices( + route_points, crs, transform + ) + return _filter_points_outside_cost_domain(route_points, shape) + + +def _get_start_end_point_cost_indices(route_points, cost_crs, transform): + """Map features to cost row, col indices""" + + logger.debug("Map %d routes to cost raster", len(route_points)) + logger.debug("First few routes:\n%s", route_points.head()) + logger.debug("Transform:\n%s", transform) + + start_lat = route_points["start_lat"].astype("float32") + start_lon = route_points["start_lon"].astype("float32") + start_row, start_col = _transform_lat_lon_to_row_col( + transform, cost_crs, start_lat, start_lon + ) + end_lat = route_points["end_lat"].astype("float32") + end_lon = route_points["end_lon"].astype("float32") + end_row, end_col = _transform_lat_lon_to_row_col( + transform, cost_crs, end_lat, end_lon + ) + + logger.debug("Mapping done!") + + route_points["start_row"] = start_row + route_points["start_col"] = start_col + route_points["end_row"] = end_row + route_points["end_col"] = end_col + + return route_points + + +def _filter_points_outside_cost_domain(route_points, shape): + """Remove points outside of cost domain""" + + logger.debug("Filtering out points outside cost domain...") + mask = route_points["start_row"] >= 0 + mask &= route_points["start_row"] < shape[0] + mask &= route_points["start_col"] >= 0 + mask &= route_points["start_col"] < shape[1] + mask &= route_points["end_row"] >= 0 + mask &= route_points["end_row"] < shape[0] + mask &= route_points["end_col"] >= 0 + mask &= route_points["end_col"] < shape[1] + + logger.debug("Mask computed!") + + if any(~mask): + msg = ( + "The following features are outside of the cost exclusion " + f"domain and will be dropped:\n{route_points.loc[~mask]}" + ) + warn(msg, revrtWarning) + route_points = route_points.loc[mask].reset_index(drop=True) + + return route_points + + +def _transform_lat_lon_to_row_col(transform, cost_crs, lat, lon): + """Transform lat/lon info to a row/col index into cost raster""" + feats = gpd.GeoDataFrame( + geometry=[Point(*p) for p in zip(lon, lat, strict=True)] + ) + coords = feats.set_crs("EPSG:4326").to_crs(cost_crs)["geometry"].centroid + row, col = rasterio.transform.rowcol( + transform, coords.x.values, coords.y.values + ) + row = np.array(row) + col = np.array(col) + return row, col From 24381416b2c4751fecb986d59f986060724047b7 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Mon, 8 Dec 2025 22:23:54 -0700 Subject: [PATCH 26/61] Update docstrings --- revrt/routing/utilities.py | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/revrt/routing/utilities.py b/revrt/routing/utilities.py index 53738551..2511ae7f 100644 --- a/revrt/routing/utilities.py +++ b/revrt/routing/utilities.py @@ -15,23 +15,24 @@ def map_to_costs(route_points, crs, transform, shape): - """_summary_ + """Map route table to cost indices and drop out-of-bounds rows Parameters ---------- - route_points : _type_ - _description_ - crs : _type_ - _description_ - transform : _type_ - _description_ - shape : _type_ - _description_ + route_points : pandas.DataFrame + Route definitions table with at least `start_lat`, `start_lon`, + `end_lat`, and `end_lon` coordinate columns. + crs : str or pyproj.CRS + Coordinate reference system for the cost raster. + transform : affine.Affine + Rasterio affine transform giving pixel origin and resolution. + shape : tuple[int, int] + Raster height and width for bounds checking. Returns ------- - _type_ - _description_ + pandas.DataFrame + Updated route table filtered to routes within the cost domain. """ route_points = _get_start_end_point_cost_indices( route_points, crs, transform @@ -40,7 +41,7 @@ def map_to_costs(route_points, crs, transform, shape): def _get_start_end_point_cost_indices(route_points, cost_crs, transform): - """Map features to cost row, col indices""" + """Populate start/end row and column indices for each route""" logger.debug("Map %d routes to cost raster", len(route_points)) logger.debug("First few routes:\n%s", route_points.head()) @@ -68,7 +69,7 @@ def _get_start_end_point_cost_indices(route_points, cost_crs, transform): def _filter_points_outside_cost_domain(route_points, shape): - """Remove points outside of cost domain""" + """Drop routes whose indices fall outside the cost domain""" logger.debug("Filtering out points outside cost domain...") mask = route_points["start_row"] >= 0 @@ -94,7 +95,7 @@ def _filter_points_outside_cost_domain(route_points, shape): def _transform_lat_lon_to_row_col(transform, cost_crs, lat, lon): - """Transform lat/lon info to a row/col index into cost raster""" + """Convert WGS84 coordinates to cost grid row and column arrays""" feats = gpd.GeoDataFrame( geometry=[Point(*p) for p in zip(lon, lat, strict=True)] ) From 54ca855da970f11c2cc108ac5fe60de06947f480 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 21:57:44 -0700 Subject: [PATCH 27/61] Add first pass of cli file --- revrt/routing/cli.py | 464 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 464 insertions(+) create mode 100644 revrt/routing/cli.py diff --git a/revrt/routing/cli.py b/revrt/routing/cli.py new file mode 100644 index 00000000..2677227d --- /dev/null +++ b/revrt/routing/cli.py @@ -0,0 +1,464 @@ +"""reVRt routing CLI functions""" + +import time +import logging +import contextlib +from math import ceil +from pathlib import Path +from copy import deepcopy + +import pandas as pd +import geopandas as gpd +import xarray as xr +from gaps.cli import CLICommandFromFunction + +from revrt.costs.config import parse_config +from revrt.routing.point_to_many import find_all_routes, RoutingScenario +from revrt.routing.utilities import map_to_costs + + +logger = logging.getLogger(__name__) +_MILLION_USD_PER_MILE_TO_USD_PER_PIXEL = 55923.40730136006 +"""Conversion from million dollars/mile to $/pixel + +1,000,000 [$/million dollars] +* 90 [meters/pixel] +/ 1609.344 [meters/mile] += 55923.40730136006 [$/pixel] +""" + + +def compute_lcp_routes( # noqa: PLR0913, PLR0917 + cost_fpath, + route_table, + cost_layers, + out_dir, + job_name, + friction_layers=None, + tracked_layers=None, + cost_multiplier_layer=None, + cost_multiplier_scalar=1, + transmission_config=None, + save_paths=False, + use_hard_barrier=False, + _split_params=None, +): + """_summary_ + + Parameters + ---------- + cost_fpath : path-like + Path to layered Zarr file containing cost and other required + routing layers. + route_table : path-like + Path to CSV file defining the start and + end points of all routes. Must have the following columns: + + "start_lat": Stating point latitude + "start_lon": Stating point longitude + "end_lat": Ending point latitude + "end_lon": Ending point longitude + + cost_layers : list + List of H5 layers that are summed to determine total costs + raster used for routing. Costs and distances for each individual + layer are also reported as requested (e.g. wet and dry costs). + Each dict in the list should have the following keys: + + - "layer_name": (REQUIRED) Name of layer in layered file + containing cost data. + - "multiplier_layer": (OPTIONAL) Name of layer in layered + file containing spatially explicit multiplier values to + apply to this cost layer before summing it with the + others. Default is ``None``. + - "multiplier_scalar": (OPTIONAL) Scalar value to multiply + this layer by before summing it with the others. Default + is ``1``. + - "is_invariant": (OPTIONAL) Boolean flag indicating whether + this layer is length invariant (i.e. should NOT be + multiplied by path length; values should be $). Default is + ``False``. + - "include_in_report": (OPTIONAL) Boolean flag indicating + whether the costs and distances for this layer should be + output in the final LCP table. Default is ``True``. + - "apply_row_mult": (OPTIONAL) Boolean flag indicating + whether the right-of-way width multiplier should be + applied for this layer. If ``True``, then the transmission + config should have a "row_width" dictionary that maps + voltages to right-of-way width multipliers. Also, the + routing table input should have a "voltage" entry for + every route. Every "voltage" value in the routing table + must be given in the "row_width" dictionary in the + transmission config, otherwise an error will be thrown. + Default is ``False``. + - "apply_polarity_mult": (OPTIONAL) Boolean flag indicating + whether the polarity multiplier should be applied for this + layer. If ``True``, then the transmission config should + have a "voltage_polarity_mult" dictionary that maps + voltages to a new dictionary, the latter mapping + polarities to multipliers. For example, a valid + "voltage_polarity_mult" dictionary might be + ``{"138": {"ac": 1.15, "dc": 2}}``. + In addition, the routing table input should have a + "voltage" **and** a "polarity" entry for every route. + Every "voltage" + "polarity" combination in the routing + table must be given in the "voltage_polarity_mult" + dictionary in the transmission config, otherwise an error + will be thrown. + + ..IMPORTANT:: + The multiplier in this config is assumed to be in units + of million $ per mile and will be converted to + $ per pixel before being applied to the layer!. + + Default is ``False``. + + out_dir : path-like + Directory to store routing outputs. + job_name : str + Name to use for this routing job. This name will be used to + construct output file names. + friction_layers : list, optional + Layers to be added to costs to influence routing but NOT + reported in final cost (i.e. friction, barriers, etc.). Each + item in this list should be a dictionary containing the + following keys: + + - "layer_name": (REQUIRED) Name of layer in layered file + containing routing data. This can also be "lcp_agg_costs", + which represents the layer built out using the + `cost_layers` input. + - "multiplier_layer": (OPTIONAL) Name of layer in layered + file containing spatially explicit multiplier values to + apply to this routing layer before summing it with the + others. Default is ``None``. This can also be + "lcp_agg_costs", which represents the layer built out + using the `cost_layers` input. + - "multiplier_scalar": (OPTIONAL) Scalar value to multiply + this layer by before summing it with the others. Default + is ``1``. + - "include_in_report": (OPTIONAL) Boolean flag indicating + whether the routing and distances for this layer should be + output in the final LCP table. Default is ``False``. + - "apply_row_mult": (OPTIONAL) Boolean flag indicating + whether the right-of-way width multiplier should be + applied for this layer. If ``True``, then the transmission + config should have a "row_width" dictionary that maps + voltages to right-of-way width multipliers. Also, the + routing table input should have a "voltage" entry for + every route. Every "voltage" value in the routing table + must be given in the "row_width" dictionary in the + transmission config, otherwise an error will be thrown. + Default is ``False``. + - "apply_polarity_mult": (OPTIONAL) Boolean flag indicating + whether the polarity multiplier should be applied for this + layer. If ``True``, then the transmission config should + have a "voltage_polarity_mult" dictionary that maps + voltages to a new dictionary, the latter mapping + polarities to multipliers. For example, a valid + "voltage_polarity_mult" dictionary might be + ``{"138": {"ac": 1.15, "dc": 2}}``. + In addition, the routing table input should have a + "voltage" **and** a "polarity" entry for every route. + Every "voltage" + "polarity" combination in the routing + table must be given in the "voltage_polarity_mult" + dictionary in the transmission config, otherwise an error + will be thrown. + + ..IMPORTANT:: + The multiplier in this config is assumed to be in units + of million $ per mile and will be converted to + $ per pixel before being applied to the layer!. + + Default is ``False``. + + By default, ``None``. + tracked_layers : dict, optional + Dictionary mapping layer names to strings, where the strings are + dask (which supports almost all numpy-equivalent) methods that + should be applied to the layer along the LCP to be included as a + characterization column in the output. By default, ``None``. + cost_multiplier_layer : str, optional + Name of layer in layered Zarr file containing final cost layer + spatial multipliers. By default, ``None``. + cost_multiplier_scalar : int, default=1 + Final cost layer multiplier. By default, ``1``. + transmission_config : path-like or dict, optional + Dictionary of transmission cost configuration values, or + path to JSON/JSON5 file containing this dictionary. The + dictionary should have a subset of the following keys: + + - base_line_costs + - iso_lookup + - iso_multipliers + - land_use_classes + - new_substation_costs + - power_classes + - power_to_voltage + - transformer_costs + - upgrade_substation_costs + - voltage_polarity_mult + - row_width + + Each of these keys should point to another dictionary or + path to JSON/JSON5 file containing a dictionary of + configurations for each section. For the expected contents + of each dictionary, see the default config. If ``None``, + values from the default config are used. + By default, ``None``. + save_paths : bool, default=False + Boolean flag to save the least cost path geometry in the output + (i.e. the output file will be a GeoPackage instead of a CSV). + By default, ``False``. + use_hard_barrier : bool, optional + Optional flag to treat any cost values of <= 0 as a hard barrier + (i.e. no paths can ever cross this). If ``False``, cost values + of <= 0 are set to a large value to simulate a strong but + permeable barrier. By default, ``False``. + + Returns + ------- + str or None + Path to file containing route data, if any routes were + completed. + """ + + start_time = time.time() + + out_dir = Path(out_dir) + out_dir.mkdir(parents=True, exist_ok=True) + + logger.debug("Tracked layers input: %r", tracked_layers) + logger.debug("Transmission config input: %r", transmission_config) + + transmission_config = parse_config(config=transmission_config) + logger.debug("Transmission Config: %s", transmission_config) + + inds = _get_indices( + route_table, + sort_cols=["start_lat", "start_lon"], + split_params=_split_params, + ) + if len(inds) == 0: + logger.info("No indices to process: %s", inds) + return None + + with contextlib.suppress(TypeError, UnicodeDecodeError): + route_points = pd.read_csv(route_table) + + if inds is not None: + route_points = route_points.loc[inds] + + out_fp = ( + out_dir / f"{job_name}.gpkg" + if save_paths + else out_dir / f"{job_name}.csv" + ) + + _run_lcp( + cost_fpath, + route_points, + cost_layers, + out_fp=out_fp, + transmission_config=transmission_config, + cost_multiplier_layer=cost_multiplier_layer, + cost_multiplier_scalar=cost_multiplier_scalar, + friction_layers=friction_layers, + tracked_layers=tracked_layers, + use_hard_barrier=use_hard_barrier, + ) + + elapsed_time = (time.time() - start_time) / 60 + logger.info("Processing took %.2f minutes", elapsed_time) + + return str(out_fp) + + +def _run_lcp( + cost_fpath, + route_points, + cost_layers, + out_fp, + transmission_config=None, + cost_multiplier_layer=None, + cost_multiplier_scalar=1, + friction_layers=None, + tracked_layers=None, + use_hard_barrier=True, +): + ts = time.monotonic() + out_fp = Path(out_fp) + save_paths = out_fp.suffix.lower() == ".gpkg" + + if route_points.empty: + logger.info("Found no paths to compute!") + + return + + with xr.open_dataset(cost_fpath, consolidated=False, engine="zarr") as ds: + route_points = map_to_costs( + route_points, + crs=ds.rio.crs, + transform=ds.rio.transform(), + shape=ds.rio.shape, + ) + + logger.info("Computing Least Cost Paths") + num_computed = 0 + for polarity, voltage, routes in _paths_to_compute(route_points, out_fp): + route_cl = _update_multipliers( + cost_layers, polarity, voltage, transmission_config + ) + route_fl = _update_multipliers( + friction_layers or [], polarity, voltage, transmission_config + ) + route_definitions = [ + ( + (info["start_row"], info["start_col"]), + [(info["end_row"], info["end_col"])], + info.to_dict(), + ) + for __, info in routes.iterrows() + ] + + scenario = RoutingScenario( + cost_fpath=cost_fpath, + cost_layers=route_cl, + friction_layers=route_fl, + tracked_layers=tracked_layers, + cost_multiplier_layer=cost_multiplier_layer, + cost_multiplier_scalar=cost_multiplier_scalar, + use_hard_barrier=use_hard_barrier, + ) + + out = find_all_routes( + scenario, + route_definitions=route_definitions, + save_paths=save_paths, + ) + + num_computed += len(out) + + if save_paths: + gpd.GeoDataFrame(out).to_file(out_fp, driver="GPKG", mode="a") + else: + pd.DataFrame(out).to_csv(out_fp, mode="a") + + time_elapsed = f"{(time.monotonic() - ts) / 3600:.4f} hour(s)" + logger.info("%d paths were computed in %s", num_computed, time_elapsed) + + +def _paths_to_compute(route_points, out_fp): + """Iterate over the paths that should be computed + + Filters out any routes that are already present in `out_fp`. If all + routes in a group already exist, the group is skipped. + """ + existing_routes = _collect_existing_routes(out_fp) + + group_cols = ["start_row", "start_col", "polarity", "voltage"] + for check_col in ["polarity", "voltage"]: + if check_col not in route_points.columns: + route_points[check_col] = "unknown" + + for group_info, routes in route_points.groupby(group_cols): + if existing_routes: + mask = routes.apply( + lambda row: ( + int(row["start_row"]), + int(row["start_col"]), + int(row["end_row"]), + int(row["end_col"]), + str(row.get("polarity", "unknown")), + str(row.get("voltage", "unknown")), + ) + not in existing_routes, + axis=1, + ) + routes = routes[mask] # noqa: PLW2901 + + if routes.empty: + continue + + *__, polarity, voltage = group_info + yield polarity, voltage, routes + + +def _collect_existing_routes(out_fp): + if out_fp is None or not out_fp.exists(): + return set() + + if out_fp.suffix.lower() == ".gpkg": + existing_df = gpd.read_file(out_fp) + else: + existing_df = pd.read_csv(out_fp) + + return { + ( + int(row["start_row"]), + int(row["start_col"]), + int(row["end_row"]), + int(row["end_col"]), + str(row.get("polarity", "unknown")), + str(row.get("voltage", "unknown")), + ) + for __, row in existing_df.iterrows() + } + + +def _update_multipliers(layers, polarity, voltage, transmission_config): + """update layer multipliers based on user input""" + output_layers = deepcopy(layers) + polarity = str(polarity) + voltage = str(int(voltage)) if voltage != "unknown" else "unknown" + + for layer in output_layers: + if layer.pop("apply_row_mult", False): + row_multiplier = transmission_config["row_width"][voltage] + layer["multiplier_scalar"] = ( + layer.get("multiplier_scalar", 1) * row_multiplier + ) + + if layer.pop("apply_polarity_mult", False): + polarity_config = transmission_config["voltage_polarity_mult"] + polarity_multiplier = polarity_config[voltage][polarity] + layer["multiplier_scalar"] = ( + layer.get("multiplier_scalar", 1) + * polarity_multiplier + * _MILLION_USD_PER_MILE_TO_USD_PER_PIXEL + ) + + return output_layers + + +def _get_indices(features_fp, sort_cols, split_params): + """Get indices of points that are sorted by location""" + features = gpd.read_file(features_fp) + features = features.sort_values(sort_cols) + + start_ind, n_chunks = split_params or (0, 1) + chunk_size = ceil(len(features) / n_chunks) + return features.index[ + start_ind * chunk_size : (start_ind + 1) * chunk_size + ] + + +def _split_routes(config): + """Compute route split params inside of config""" + exec_control = config.get("execution_control", {}) + if exec_control.get("option") == "local": + num_nodes = 1 + else: + num_nodes = exec_control.pop("nodes", 1) + + config["_split_params"] = [(i, num_nodes) for i in range(num_nodes)] + return config + + +route_points_command = CLICommandFromFunction( + compute_lcp_routes, + name="route-points", + add_collect=False, + split_keys={"_split_params"}, + config_preprocessor=_split_routes, +) From 57f558e22e35b859d4a3594b4caa5986a9e5dbae Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:04:04 -0700 Subject: [PATCH 28/61] Minor updates --- revrt/routing/cli.py | 36 +++++++++++++++++------------------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/revrt/routing/cli.py b/revrt/routing/cli.py index 2677227d..bb7af50c 100644 --- a/revrt/routing/cli.py +++ b/revrt/routing/cli.py @@ -43,7 +43,7 @@ def compute_lcp_routes( # noqa: PLR0913, PLR0917 use_hard_barrier=False, _split_params=None, ): - """_summary_ + """Run least-cost path routing for pairs of points Parameters ---------- @@ -114,10 +114,9 @@ def compute_lcp_routes( # noqa: PLR0913, PLR0917 Default is ``False``. out_dir : path-like - Directory to store routing outputs. + Directory where routing outputs should be written. job_name : str - Name to use for this routing job. This name will be used to - construct output file names. + Label used to name the generated output file. friction_layers : list, optional Layers to be added to costs to influence routing but NOT reported in final cost (i.e. friction, barriers, etc.). Each @@ -175,14 +174,15 @@ def compute_lcp_routes( # noqa: PLR0913, PLR0917 By default, ``None``. tracked_layers : dict, optional Dictionary mapping layer names to strings, where the strings are - dask (which supports almost all numpy-equivalent) methods that + dask aggregation methods (similar to what numpy has) that should be applied to the layer along the LCP to be included as a characterization column in the output. By default, ``None``. cost_multiplier_layer : str, optional - Name of layer in layered Zarr file containing final cost layer - spatial multipliers. By default, ``None``. + Name of the spatial multiplier layer applied to final costs. + By default, ``None``. cost_multiplier_scalar : int, default=1 - Final cost layer multiplier. By default, ``1``. + Scalar multiplier applied to the final cost surface. + By default, ``1``. transmission_config : path-like or dict, optional Dictionary of transmission cost configuration values, or path to JSON/JSON5 file containing this dictionary. The @@ -207,9 +207,8 @@ def compute_lcp_routes( # noqa: PLR0913, PLR0917 values from the default config are used. By default, ``None``. save_paths : bool, default=False - Boolean flag to save the least cost path geometry in the output - (i.e. the output file will be a GeoPackage instead of a CSV). - By default, ``False``. + Save outputs as a GeoPackage with path geometries when ``True``. + Defaults to ``False``. use_hard_barrier : bool, optional Optional flag to treat any cost values of <= 0 as a hard barrier (i.e. no paths can ever cross this). If ``False``, cost values @@ -219,8 +218,7 @@ def compute_lcp_routes( # noqa: PLR0913, PLR0917 Returns ------- str or None - Path to file containing route data, if any routes were - completed. + Path to the output table if any routes were computed. """ start_time = time.time() @@ -286,6 +284,8 @@ def _run_lcp( tracked_layers=None, use_hard_barrier=True, ): + """Execute least-cost path routing for the prepared route subset""" + ts = time.monotonic() out_fp = Path(out_fp) save_paths = out_fp.suffix.lower() == ".gpkg" @@ -349,11 +349,7 @@ def _run_lcp( def _paths_to_compute(route_points, out_fp): - """Iterate over the paths that should be computed - - Filters out any routes that are already present in `out_fp`. If all - routes in a group already exist, the group is skipped. - """ + """Yield route groups that still require computation""" existing_routes = _collect_existing_routes(out_fp) group_cols = ["start_row", "start_col", "polarity", "voltage"] @@ -385,6 +381,8 @@ def _paths_to_compute(route_points, out_fp): def _collect_existing_routes(out_fp): + """Collect already computed routes from an existing output file""" + if out_fp is None or not out_fp.exists(): return set() @@ -407,7 +405,7 @@ def _collect_existing_routes(out_fp): def _update_multipliers(layers, polarity, voltage, transmission_config): - """update layer multipliers based on user input""" + """Update layer multipliers based on user input""" output_layers = deepcopy(layers) polarity = str(polarity) voltage = str(int(voltage)) if voltage != "unknown" else "unknown" From d82b8ae62691bb18cead2507f7ca9f552f01a4e3 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:04:12 -0700 Subject: [PATCH 29/61] Connect first routing CLI --- revrt/_cli.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/revrt/_cli.py b/revrt/_cli.py index cb017d5a..2849120e 100644 --- a/revrt/_cli.py +++ b/revrt/_cli.py @@ -7,6 +7,7 @@ from revrt import __version__ from revrt.spatial_characterization.cli import route_characterizations_command from revrt.costs.cli import build_routing_layers_command +from revrt.routing.cli import route_points_command from revrt.utilities.cli import ( layers_to_file_command, layers_from_file_command, @@ -21,6 +22,7 @@ layers_from_file_command, build_routing_layers_command, route_characterizations_command, + route_points_command, ] main = make_cli(commands, info={"name": "reVRt", "version": __version__}) From fe8f4aaa203a65641e4dd04e8b1d5e7d26e11e6b Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:09:05 -0700 Subject: [PATCH 30/61] Add changes to points logic --- revrt/routing/point_to_many.py | 101 ++++++++++++++++++++------------- 1 file changed, 63 insertions(+), 38 deletions(-) diff --git a/revrt/routing/point_to_many.py b/revrt/routing/point_to_many.py index 1d1aceca..ea45c123 100644 --- a/revrt/routing/point_to_many.py +++ b/revrt/routing/point_to_many.py @@ -18,8 +18,7 @@ from revrt import find_paths from revrt.exceptions import ( - # revrtValueError, - # revrtInvalidStartCostError, + revrtKeyError, revrtLeastCostPathNotFoundError, ) from revrt.warn import revrtWarning @@ -121,8 +120,7 @@ def _verify_layer_exists(self, layer_name): f"Did not find layer {layer_name!r} in cost " f"file {str(self.routing_scenario.cost_fpath)!r}" ) - logger.error(msg) - raise KeyError(msg) + raise revrtKeyError(msg) def close(self): """Close the underlying xarray file handle""" @@ -372,7 +370,12 @@ class RouteResult: """Class to compute route characteristics given layer cost maps""" def __init__( - self, routing_layers, route, optimized_objective, add_geom=False + self, + routing_layers, + route, + optimized_objective, + add_geom=False, + attrs=None, ): self._routing_layers = routing_layers self._route = route @@ -380,11 +383,11 @@ def __init__( self._lens = self._total_path_length = None self._by_layer_results = {} self._add_geom = add_geom + self._attrs = attrs or {} @property def cell_size(self): return abs(self._routing_layers.transform.a) - # return self._routing_layers.cell_size @property def lens(self): @@ -444,10 +447,25 @@ def build(self): "cost": self.cost, "poi_lat": self.end_lat, "poi_lon": self.end_lon, + "start_row": self._route[0][0], + "start_col": self._route[0][1], "end_row": self._route[-1][0], "end_col": self._route[-1][1], "optimized_objective": self._optimized_objective, } + for check_key in ["start_row", "start_col", "end_row", "end_col"]: + if ( + check_key in self._attrs + and results[check_key] != self._attrs[check_key] + ): + msg = ( + f"Computed {check_key}={results[check_key]} does " + f"not match expected {check_key}=" + f"{self._attrs[check_key]}!" + ) + warn(msg, revrtWarning) + + results.update(self._attrs) for layer in self._routing_layers.tracked_layers: layer_result = layer.compute(self._route, self.cell_size) results.update(layer_result) @@ -465,15 +483,21 @@ def _compute_path_length(self): def find_all_routes(routing_scenario, route_definitions, save_paths=False): + if not route_definitions: + return [] + ts = time.monotonic() routing_layers = RoutingLayers(routing_scenario).build() - routes = _compute_routes( - routing_scenario, - route_definitions, - routing_layers, - save_paths=save_paths, - ) + try: + routes = _compute_routes( + routing_scenario, + route_definitions, + routing_layers, + save_paths=save_paths, + ) + finally: + routing_layers.close() time_elapsed = f"{(time.monotonic() - ts) / 60:.4f} min" logger.debug("Least Cost tie-line computed in %s", time_elapsed) @@ -485,7 +509,7 @@ def _compute_routes( routing_scenario, route_definitions, routing_layers, save_paths ): routes = [] - for start_point, end_points in route_definitions: + for start_point, end_points, attrs in route_definitions: try: indices, optimized_objective = _compute_valid_path( routing_scenario, routing_layers, start_point, end_points @@ -498,24 +522,39 @@ def _compute_routes( indices, optimized_objective, add_geom=save_paths, + attrs=attrs, ) routes.append(route.build()) - if save_paths: - if routes: - return gpd.GeoDataFrame( - routes, geometry="geometry", crs=routing_layers.cost_crs - ) - return gpd.GeoDataFrame([]) - - return pd.DataFrame(routes) + return routes def _compute_valid_path( routing_scenario, routing_layers, start_point, end_points ): + _validate_starting_point(routing_layers, start_point) + _validate_end_points(routing_layers, end_points) + + try: + route_result = find_paths( + zarr_fp=routing_scenario.cost_fpath, + cost_layers=routing_scenario.cl_as_json, + start=[start_point], + end=end_points, + )[0] + except Exception as ex: + msg = ( + f"Unable to find path from {start_point} any of {end_points}: {ex}" + ) + logger.exception(msg) + + raise revrtLeastCostPathNotFoundError(msg) from ex + + return route_result + + +def _validate_starting_point(routing_layers, start_point): start_row, start_col = start_point - # print(f"{routing_layers.final_routing_layer=}") start_cost = ( routing_layers.final_routing_layer.isel(y=start_row, x=start_col) .compute() @@ -529,6 +568,8 @@ def _compute_valid_path( ) raise revrtLeastCostPathNotFoundError(msg) + +def _validate_end_points(routing_layers, end_points): end_rows, end_cols = np.array(end_points).T end_costs = routing_layers.final_routing_layer.isel( y=xr.DataArray(end_rows, dims="points"), @@ -541,22 +582,6 @@ def _compute_valid_path( ) raise revrtLeastCostPathNotFoundError(msg) - try: - indices, optimized_objective = find_paths( - zarr_fp=routing_scenario.cost_fpath, - cost_layers=routing_scenario.cl_as_json, - start=[start_point], - end=end_points, - )[0] - except Exception as ex: - msg = ( - f"Unable to find path from {start_point} any of {end_points}: {ex}" - ) - - raise revrtLeastCostPathNotFoundError(msg) from ex - - return indices, optimized_objective - def _compute_lens(route, cell_size): """Compute the total length and cell by cell length of LCP""" From fa64388344922819261735a4984317c2ee2ca186 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:09:15 -0700 Subject: [PATCH 31/61] Minor docstring updates --- revrt/costs/config/config.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/revrt/costs/config/config.py b/revrt/costs/config/config.py index 3c1317cc..768eef43 100644 --- a/revrt/costs/config/config.py +++ b/revrt/costs/config/config.py @@ -54,7 +54,7 @@ def __init__(self, config=None): config : str or dict, optional Dictionary of transmission cost configuration values, or path to JSON/JSON5 file containing this dictionary. The - dictionary should have the following keys: + dictionary should have a subset of the following keys: - base_line_costs - iso_lookup @@ -65,6 +65,8 @@ def __init__(self, config=None): - power_to_voltage - transformer_costs - upgrade_substation_costs + - voltage_polarity_mult + - row_width Each of these keys should point to another dictionary or path to JSON/JSON5 file containing a dictionary of From 2faaa31e83452a4e89caa42c2440e1b0fcd1e45b Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:15:59 -0700 Subject: [PATCH 32/61] Update --- revrt/routing/cli.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/revrt/routing/cli.py b/revrt/routing/cli.py index bb7af50c..e9d8476d 100644 --- a/revrt/routing/cli.py +++ b/revrt/routing/cli.py @@ -340,7 +340,9 @@ def _run_lcp( num_computed += len(out) if save_paths: - gpd.GeoDataFrame(out).to_file(out_fp, driver="GPKG", mode="a") + gpd.GeoDataFrame(out, geometry="geometry").to_file( + out_fp, driver="GPKG", mode="a" + ) else: pd.DataFrame(out).to_csv(out_fp, mode="a") From b16bef050a8d62b03af7d84d575c36728f1087b0 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:16:25 -0700 Subject: [PATCH 33/61] Partially fix tests --- .../routing/test_routing_point_to_many.py | 66 ++++++++++++------- 1 file changed, 41 insertions(+), 25 deletions(-) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index adc88107..59f1ab00 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -143,7 +143,11 @@ def test_basic_single_route_layered_file_short_path(sample_layered_data): ) output = find_all_routes( - scenario, route_definitions=[((1, 1), [(1, 2)])], save_paths=False + scenario, + route_definitions=[ + ((1, 1), [(1, 2)], {}), + ], + save_paths=False, ) assert len(output) == 1 @@ -163,8 +167,8 @@ def test_basic_single_route_layered_file(sample_layered_data): output = find_all_routes( scenario, route_definitions=[ - ((1, 1), [(2, 6)]), - ((1, 2), [(2, 6)]), + ((1, 1), [(2, 6)], {}), + ((1, 2), [(2, 6)], {}), ], save_paths=False, ) @@ -196,7 +200,10 @@ def test_multi_layer_route_layered_file(sample_layered_data): output = find_all_routes( scenario, - route_definitions=[((1, 1), [(2, 6)]), ((1, 2), [(2, 6)])], + route_definitions=[ + ((1, 1), [(2, 6)], {}), + ((1, 2), [(2, 6)], {}), + ], save_paths=False, ) @@ -256,15 +263,14 @@ def test_save_paths_returns_expected_geometry(sample_layered_data): output = find_all_routes( scenario, route_definitions=[ - ((1, 1), [(2, 6)]), - ((1, 2), [(2, 6)]), + ((1, 1), [(2, 6)], {}), + ((1, 2), [(2, 6)], {}), ], save_paths=True, ) - assert isinstance(output, gpd.GeoDataFrame) + assert isinstance(output, list) assert len(output) == 2 - assert output.crs.to_string() == "EPSG:4326" expected_geometries = [ [ @@ -312,9 +318,8 @@ def test_empty_route_definitions_returns_empty_dataframe(sample_layered_data): save_paths=False, ) - assert isinstance(output, pd.DataFrame) - assert output.empty - assert list(output.columns) == [] + assert isinstance(output, list) + assert not output def test_empty_route_definitions_returns_empty_geo_dataframe( @@ -333,9 +338,8 @@ def test_empty_route_definitions_returns_empty_geo_dataframe( save_paths=True, ) - assert isinstance(output, gpd.GeoDataFrame) - assert output.empty - assert list(output.columns) == [] + assert isinstance(output, list) + assert not output def test_multi_layer_route_with_multiplier(sample_layered_data): @@ -355,8 +359,8 @@ def test_multi_layer_route_with_multiplier(sample_layered_data): output = find_all_routes( scenario, route_definitions=[ - ((1, 1), [(2, 6)]), - ((1, 2), [(2, 6)]), + ((1, 1), [(2, 6)], {}), + ((1, 2), [(2, 6)], {}), ], save_paths=False, ) @@ -434,7 +438,9 @@ def test_multi_layer_route_with_scalar_and_layer_multipliers( output = find_all_routes( scenario, - route_definitions=[((1, 1), [(1, 2)])], + route_definitions=[ + ((1, 1), [(1, 2)], {}), + ], save_paths=False, ) @@ -469,7 +475,9 @@ def test_routing_with_tracked_layers(sample_layered_data): output = find_all_routes( scenario, - route_definitions=[((1, 1), [(1, 2)])], + route_definitions=[ + ((1, 1), [(1, 2)], {}), + ], save_paths=False, ) @@ -499,15 +507,19 @@ def test_start_point_on_barrier_returns_no_route( # (0, 3) in layer_1 is 0 -> treated as barrier output = find_all_routes( - scenario, route_definitions=[((0, 3), [(2, 6)])], save_paths=False + scenario, + route_definitions=[ + ((0, 3), [(2, 6)], {}), + ], + save_paths=False, ) assert_message_was_logged( "Start idx (0, 3) does not have a valid cost: -1.00 (must be > 0)!", "ERROR", ) - assert isinstance(output, pd.DataFrame) - assert output.empty + assert isinstance(output, list) + assert not output def test_some_endpoints_include_barriers_but_one_valid(sample_layered_data): @@ -521,7 +533,9 @@ def test_some_endpoints_include_barriers_but_one_valid(sample_layered_data): # include one barrier end (0,3) and one valid end (2,6) output = find_all_routes( scenario, - route_definitions=[((1, 1), [(0, 3), (2, 6)])], + route_definitions=[ + ((1, 1), [(0, 3), (2, 6)], {}), + ], save_paths=False, ) @@ -546,7 +560,9 @@ def test_all_endpoints_are_barriers_returns_no_route( output = find_all_routes( scenario, - route_definitions=[((1, 1), [(0, 3), (0, 7)])], + route_definitions=[ + ((1, 1), [(0, 3), (0, 7)], {}), + ], save_paths=False, ) assert_message_was_logged( @@ -554,8 +570,8 @@ def test_all_endpoints_are_barriers_returns_no_route( "ERROR", ) - assert isinstance(output, pd.DataFrame) - assert output.empty + assert isinstance(output, list) + assert not output if __name__ == "__main__": From e9fa52120fe4e0de38d78477ce926a1dec0c3ab3 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:19:38 -0700 Subject: [PATCH 34/61] Fix tests --- .../routing/test_routing_point_to_many.py | 52 ++++++++++--------- 1 file changed, 27 insertions(+), 25 deletions(-) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index 59f1ab00..a0c54329 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -2,8 +2,6 @@ from pathlib import Path -import geopandas as gpd -import pandas as pd import pytest import numpy as np import xarray as xr @@ -151,9 +149,10 @@ def test_basic_single_route_layered_file_short_path(sample_layered_data): ) assert len(output) == 1 - assert output.iloc[0]["cost"] == pytest.approx((1 + 2) / 2) - assert output.iloc[0]["length_km"] == 1 / 1000 - assert output.iloc[0]["cost"] == output.iloc[0]["optimized_objective"] + route = output[0] + assert route["cost"] == pytest.approx((1 + 2) / 2) + assert route["length_km"] == 1 / 1000 + assert route["cost"] == route["optimized_objective"] def test_basic_single_route_layered_file(sample_layered_data): @@ -174,16 +173,16 @@ def test_basic_single_route_layered_file(sample_layered_data): ) assert len(output) == 2 - assert output.iloc[0]["cost"] == pytest.approx(11.192389) - assert output.iloc[0]["length_km"] == pytest.approx(0.0090710678) + first_route = output[0] + assert first_route["cost"] == pytest.approx(11.192389) + assert first_route["length_km"] == pytest.approx(0.0090710678) + assert np.isclose(first_route["cost"], first_route["optimized_objective"]) + + second_route = output[1] + assert second_route["cost"] == pytest.approx(12.278174) + assert second_route["length_km"] == pytest.approx(0.008656854) assert np.isclose( - output.iloc[0]["cost"], output.iloc[0]["optimized_objective"] - ) - - assert output.iloc[1]["cost"] == pytest.approx(12.278174) - assert output.iloc[1]["length_km"] == pytest.approx(0.008656854) - assert np.isclose( - output.iloc[1]["cost"], output.iloc[1]["optimized_objective"] + second_route["cost"], second_route["optimized_objective"] ) @@ -209,7 +208,7 @@ def test_multi_layer_route_layered_file(sample_layered_data): assert len(output) == 2 - first_route = output.iloc[0] + first_route = output[0] assert first_route["cost"] == pytest.approx( 27.606602, rel=1e-4, @@ -230,7 +229,7 @@ def test_multi_layer_route_layered_file(sample_layered_data): first_route["cost"], first_route["optimized_objective"], rtol=1e-6 ) - second_route = output.iloc[1] + second_route = output[1] assert second_route["cost"] == pytest.approx( 25.106602, rel=1e-4, @@ -272,6 +271,8 @@ def test_save_paths_returns_expected_geometry(sample_layered_data): assert isinstance(output, list) assert len(output) == 2 + route_geoms = [route["geometry"] for route in output] + expected_geometries = [ [ (1.5, 5.5), @@ -296,7 +297,7 @@ def test_save_paths_returns_expected_geometry(sample_layered_data): ] for geom, expected_coords in zip( - output.geometry, expected_geometries, strict=True + route_geoms, expected_geometries, strict=True ): assert geom.geom_type == "LineString" assert np.allclose( @@ -367,7 +368,7 @@ def test_multi_layer_route_with_multiplier(sample_layered_data): assert len(output) == 2 - first_route = output.iloc[0] + first_route = output[0] assert first_route["cost"] == pytest.approx( 22.588835, rel=1e-4, @@ -390,7 +391,7 @@ def test_multi_layer_route_with_multiplier(sample_layered_data): rtol=1e-4, ) - second_route = output.iloc[1] + second_route = output[1] assert second_route["cost"] == pytest.approx( 20.588835, rel=1e-4, @@ -446,7 +447,7 @@ def test_multi_layer_route_with_scalar_and_layer_multipliers( assert len(output) == 1 - route = output.iloc[0] + route = output[0] assert route["cost"] == pytest.approx(2.0, rel=1e-4) assert route["length_km"] == pytest.approx(0.001, rel=1e-4) assert route["layer_1_cost"] == pytest.approx(1.5, rel=1e-4) @@ -482,13 +483,13 @@ def test_routing_with_tracked_layers(sample_layered_data): ) assert len(output) == 1 - route = output.iloc[0] + route = output[0] assert { "layer_1_mean", "layer_2_max", "layer_3_min", - }.issubset(output.columns) + }.issubset(route.keys()) assert route["layer_1_mean"] == pytest.approx(1.5) assert route["layer_2_max"] == pytest.approx(1.0) @@ -541,10 +542,11 @@ def test_some_endpoints_include_barriers_but_one_valid(sample_layered_data): assert len(output) == 1 # At least one valid endpoint must be reached and cost must be positive. - assert output.iloc[0]["cost"] > 0 + route = output[0] + assert route["cost"] > 0 - end_row = int(output.iloc[0]["end_row"]) - end_col = int(output.iloc[0]["end_col"]) + end_row = int(route["end_row"]) + end_col = int(route["end_col"]) assert (end_row, end_col) == (2, 6) From ed9162709987e21922799370106b4a6aa93410c9 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:42:04 -0700 Subject: [PATCH 35/61] Add docstrings --- revrt/routing/point_to_many.py | 216 ++++++++++++++++++++------------- 1 file changed, 130 insertions(+), 86 deletions(-) diff --git a/revrt/routing/point_to_many.py b/revrt/routing/point_to_many.py index ea45c123..282c4628 100644 --- a/revrt/routing/point_to_many.py +++ b/revrt/routing/point_to_many.py @@ -8,10 +8,8 @@ import rasterio import numpy as np -import pandas as pd import xarray as xr import dask.array as da -import geopandas as gpd from shapely.geometry import Point from shapely.geometry.linestring import LineString @@ -29,6 +27,8 @@ class RoutingScenario: + """Container for routing scenario configuration""" + def __init__( self, cost_fpath, @@ -39,6 +39,28 @@ def __init__( cost_multiplier_scalar=1, use_hard_barrier=True, ): + """ + + Parameters + ---------- + cost_fpath : path-like + Path to the cost layer Zarr store used for routing. + cost_layers : list + List of dictionaries containing layer definitions + contributing to the summed routing cost. + friction_layers : list, optional + List of dictionaries defining layers that influence routing + but are excluded from reports. + tracked_layers : dict, optional + Layers to summarize along the path, mapped to aggregation + names. + cost_multiplier_layer : str, optional + Layer name providing spatial multipliers for total cost. + cost_multiplier_scalar : int or float, optional + Scalar multiplier applied to the final cost surface. + use_hard_barrier : bool, optional + Flag indicating whether non-positive costs block traversal. + """ self.cost_fpath = cost_fpath self.cost_layers = cost_layers self.friction_layers = friction_layers or [] @@ -47,8 +69,18 @@ def __init__( self.cost_multiplier_scalar = cost_multiplier_scalar self.use_hard_barrier = use_hard_barrier + def __repr__(self): + return ( + "RoutingScenario:" + f"\n\t- cost_layers: {self.cost_layers}" + f"\n\t- friction_layers: {self.friction_layers}" + f"\n\t- cost_multiplier_layer: {self.cost_multiplier_layer}" + f"\n\t- cost_multiplier_scalar: {self.cost_multiplier_scalar}" + ) + @cached_property def cl_as_json(self): + """str: JSON string describing configured cost layers""" return json.dumps({"cost_layers": self.cost_layers}) @@ -66,22 +98,11 @@ def __init__(self, routing_scenario, chunks="auto"): Parameters ---------- - cost_fpath : str - Full path to HDF5 file containing cost arrays. The cost data - layers should be named ``"tie_line_costs_{capacity}MW"``, - where ``capacity`` is an integer value representing the - capacity of the line (the integer values must matches at - least one of the integer capacity values represented by the - keys in the ``power_classes`` portion of the transmission - config). - cell_size : int, optional - Side length of each cell, in meters. Cells are assumed to be - square. By default, :obj:`CELL_SIZE`. - use_hard_barrier : bool, optional - Optional flag to treat any cost values of <= 0 as a hard - barrier (i.e. no paths can ever cross this). If ``False``, - cost values of <= 0 are set to a large value to simulate a - strong but permeable barrier. By default, ``True``. + routing_scenario : RoutingScenario + Scenario containing cost, friction, and tracking metadata. + chunks : str or mapping, default="auto" + Chunk specification forwarded to ``xarray.open_dataset``. + By default, ``"auto"``. """ self.routing_scenario = routing_scenario self._layer_fh = xr.open_dataset( @@ -90,27 +111,29 @@ def __init__(self, routing_scenario, chunks="auto"): consolidated=False, engine="zarr", ) - # self.cost_layer_map = {} self.li_cost_layer_map = {} - # self.friction_layer_map = {} self.tracked_layers = [] self.transform = self._layer_fh.rio.transform() self._full_shape = self._layer_fh.rio.shape self.cost_crs = self._layer_fh.rio.crs - # self.cell_size = abs(self.transform.a) self._layers = set(self._layer_fh.variables) self.cost = None self.li_cost = None self.final_routing_layer = None + def __repr__(self): + return f"RoutingLayers for {self.routing_scenario!r}" + @property def latitudes(self): + """array-like: Latitude coordinates for the cost grid""" return self._layer_fh["latitude"] @property def longitudes(self): + """array-like: Longitude coordinates for the cost grid""" return self._layer_fh["longitude"] def _verify_layer_exists(self, layer_name): @@ -127,58 +150,9 @@ def close(self): self._layer_fh.close() def build(self): - """Extract clipped cost arrays from exclusion files - - Parameters - ---------- - cost_layers : List[str] - List of layers in H5 that are summed to determine total - costs raster used for routing. Costs and distances for each - individual layer are also reported (e.g. wet and dry costs). - determining path using main cost layer. - friction_layers : List[dict] | None, optional - List of layers in H5 to be added to the cost raster to - influence routing but NOT reported in final cost. Should - have the same format as the `cost_layers` input. - By default, ``None``. - tracked_layers : dict, optional - Dictionary mapping layer names to strings, where the strings - are numpy methods that should be applied to the layer along - the LCP. For example, - ``tracked_layers={'layer_1': 'mean', 'layer_2': 'max}`` - would report the average of ``layer_1`` values along the - least cost path and the max of ``layer_2`` values along the - least cost path. Examples of numpy methods (non-exhaustive): - - - mean - - max - - min - - mode - - median - - std - - By default, ``None``, which does not track any extra layers. - cost_multiplier_layer : str, optional - Name of layer containing final cost layer spatial - multipliers. By default, ``None``. - cost_multiplier_scalar : int | float, optional - Final cost layer multiplier. By default, ``1``. - """ - - logger.debug( - # TODO: Turn this input a repr of routing scenario and use - # that here - "Building cost layer with the following inputs:" - "\n\t- cost_layers: %r" - "\n\t- friction_layers: %r" - "\n\t- cost_multiplier_layer: %r" - "\n\t- cost_multiplier_scalar: %r", - self.routing_scenario.cost_layers, - self.routing_scenario.friction_layers, - self.routing_scenario.cost_multiplier_layer, - self.routing_scenario.cost_multiplier_scalar, - ) + """Build lazy routing arrays from layered file""" + logger.debug("Building %r", self) self._build_cost_layer() self._build_final_routing_layer() self._build_tracked_layers() @@ -219,17 +193,12 @@ def _build_final_routing_layer(self): self.final_routing_layer = self.cost.copy() self.final_routing_layer += self.li_cost - # for li_cost_layer in self.li_cost_layer_map.values(): - # # normalize by cell size for routing only - # self.final_routing_layer += li_cost_layer / self.cell_size - friction_costs = da.zeros(self._full_shape, dtype=np.float32) for layer_info in self.routing_scenario.friction_layers: layer_name = layer_info["layer_name"] friction_layer = self._extract_and_scale_layer( layer_info, allow_cl=True ) - # self.friction_layer_map[layer_name] = friction_layer if layer_info.get("include_in_report", False): self.tracked_layers.append( CharacterizedLayer(layer_name, friction_layer) @@ -249,12 +218,6 @@ def _build_final_routing_layer(self): -1 if self.routing_scenario.use_hard_barrier else max_val, self.final_routing_layer, ) + (self.cost * 0) - # logger.debug( - # "MCP cost min: %.2f, max: %.2f, median: %.2f", - # np.min(self.final_routing_layer), - # np.max(self.final_routing_layer), - # np.median(self.final_routing_layer), - # ) def _extract_and_scale_layer(self, layer_info, allow_cl=False): """Extract layer based on name and scale according to input""" @@ -293,8 +256,9 @@ def _build_tracked_layers(self): if tracked_layer not in self._layers: msg = ( - f"Did not find layer {tracked_layer!r} in cost " - f"file {str(self.routing_scenario.cost_fpath)!r}. Skipping..." + f"Did not find layer {tracked_layer!r} in cost file " + f"{str(self.routing_scenario.cost_fpath)!r}. " + "Skipping..." ) warn(msg, revrtWarning) continue @@ -308,9 +272,26 @@ def _build_tracked_layers(self): class CharacterizedLayer: + """Encapsulate tracked routing layer metadata""" + def __init__( self, name, layer, is_length_invariant=False, agg_method=None ): + """ + + Parameters + ---------- + name : str + Identifier used when reporting layer-derived metrics. + layer : xarray.DataArray or dask.array.Array + Data values sampled from the routing stack. + is_length_invariant : bool, default=False + Flag signaling cost values should ignore segment length. + By default, ``False``. + agg_method : str, optional + Name of dask aggregation used when tracking statistics. + By default, ``None``. + """ self.name = name self.layer = layer self.is_length_invariant = is_length_invariant @@ -324,7 +305,20 @@ def __repr__(self): ) def compute(self, route, cell_size): - """Compute aggregate value over route""" + """Compute layer metrics along a route + + Parameters + ---------- + route : sequence + Ordered ``(row, col)`` indices describing the path. + cell_size : float + Raster cell size in meters for distance calculations. + + Returns + ------- + dict + Mapping of metric names to aggregated layer values. + """ rows, cols = np.array(route).T layer_values = self.layer.isel( y=xr.DataArray(rows, dims="points"), @@ -339,6 +333,7 @@ def compute(self, route, cell_size): return self._compute_agg(layer_values) def _compute_total_and_length(self, layer_values, route, cell_size): + """Compute total cost and length metrics for the layer""" lens, __ = _compute_lens(route, cell_size) layer_data = getattr(layer_values, "data", layer_values) @@ -360,6 +355,7 @@ def _compute_total_and_length(self, layer_values, route, cell_size): } def _compute_agg(self, layer_values): + """Compute aggregated statistic for tracked layer""" aggregate = getattr(da, self.agg_method)(layer_values).astype( np.float32 ) @@ -377,6 +373,23 @@ def __init__( add_geom=False, attrs=None, ): + """ + + Parameters + ---------- + routing_layers : RoutingLayers + Routing layer manager containing cost and tracker arrays. + route : list + Ordered row and column indices defining the path. + optimized_objective : float + Objective value returned by the routing solver. + add_geom : bool, default=False + Include shapely geometry in the output when ``True``. + By default, ``False``. + attrs : dict, optional + Additional metadata merged into the result dictionary. + By default, ``None``. + """ self._routing_layers = routing_layers self._route = route self._optimized_objective = optimized_objective @@ -387,22 +400,26 @@ def __init__( @property def cell_size(self): + """float: Raster cell size in meters""" return abs(self._routing_layers.transform.a) @property def lens(self): + """array-like: Cached per-cell travel distances""" if self._lens is None: self._compute_path_length() return self._lens @property def total_path_length(self): + """float: Total path length in kilometers""" if self._total_path_length is None: self._compute_path_length() return self._total_path_length @property def cost(self): + """float: Optimized objective evaluated over the route""" rows, cols = np.array(self._route).T cell_costs = self._routing_layers.cost.isel( y=xr.DataArray(rows, dims="points"), @@ -414,6 +431,7 @@ def cost(self): @property def end_lat(self): + """float: Latitude of the terminal path cell""" row, col = self._route[-1] return ( self._routing_layers.latitudes.isel(y=row, x=col) @@ -424,6 +442,7 @@ def end_lat(self): @property def end_lon(self): + """float: Longitude of the terminal path cell""" row, col = self._route[-1] return ( self._routing_layers.longitudes.isel(y=row, x=col) @@ -434,6 +453,7 @@ def end_lon(self): @property def geom(self): + """shapely.GeometryType: Geometry for the computed path""" rows, cols = np.array(self._route).T x, y = rasterio.transform.xy( self._routing_layers.transform, rows, cols @@ -442,6 +462,7 @@ def geom(self): return geom(list(zip(x, y, strict=True))) def build(self): + """Assemble route metrics and optional geometry payload""" results = { "length_km": self.total_path_length, "cost": self.cost, @@ -483,6 +504,25 @@ def _compute_path_length(self): def find_all_routes(routing_scenario, route_definitions, save_paths=False): + """Compute least-cost routes for each start and destination set + + Parameters + ---------- + routing_scenario : RoutingScenario + Scenario describing the cost layers and routing options. + route_definitions : Iterable + Sequence of ``(start_point, end_points, attrs)`` tuples + defining which points to route between. The `attrs` dictionary + will be appended to the output for each route. + save_paths : bool, default=False + Include shapely geometries in the output when ``True``. + By default, ``False``. + + Returns + ------- + list + Route summaries for each successfully computed path. + """ if not route_definitions: return [] @@ -508,6 +548,7 @@ def find_all_routes(routing_scenario, route_definitions, save_paths=False): def _compute_routes( routing_scenario, route_definitions, routing_layers, save_paths ): + """Evaluate route definitions and build result records""" routes = [] for start_point, end_points, attrs in route_definitions: try: @@ -532,6 +573,7 @@ def _compute_routes( def _compute_valid_path( routing_scenario, routing_layers, start_point, end_points ): + """Validate provided indices then solve for the least-cost path""" _validate_starting_point(routing_layers, start_point) _validate_end_points(routing_layers, end_points) @@ -554,6 +596,7 @@ def _compute_valid_path( def _validate_starting_point(routing_layers, start_point): + """Raise when the starting cell lacks a positive traversal cost""" start_row, start_col = start_point start_cost = ( routing_layers.final_routing_layer.isel(y=start_row, x=start_col) @@ -570,6 +613,7 @@ def _validate_starting_point(routing_layers, start_point): def _validate_end_points(routing_layers, end_points): + """Raise when no end cell provides a positive traversal cost""" end_rows, end_cols = np.array(end_points).T end_costs = routing_layers.final_routing_layer.isel( y=xr.DataArray(end_rows, dims="points"), From 2fe36fbc42ed86d20f917464cf4888338e37ddb5 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:54:11 -0700 Subject: [PATCH 36/61] Minor updates --- revrt/routing/point_to_many.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/revrt/routing/point_to_many.py b/revrt/routing/point_to_many.py index 282c4628..6bb14922 100644 --- a/revrt/routing/point_to_many.py +++ b/revrt/routing/point_to_many.py @@ -128,12 +128,12 @@ def __repr__(self): @property def latitudes(self): - """array-like: Latitude coordinates for the cost grid""" + """xarray.DataArray: Latitude coordinates for the cost grid""" return self._layer_fh["latitude"] @property def longitudes(self): - """array-like: Longitude coordinates for the cost grid""" + """xarray.DataArray: Longitude coordinates for the cost grid""" return self._layer_fh["longitude"] def _verify_layer_exists(self, layer_name): @@ -404,7 +404,7 @@ def cell_size(self): return abs(self._routing_layers.transform.a) @property - def lens(self): + def _lens(self): """array-like: Cached per-cell travel distances""" if self._lens is None: self._compute_path_length() @@ -427,7 +427,7 @@ def cost(self): ) # Multiple distance travel through cell by cost and sum it! - return da.sum(cell_costs * self.lens).astype(np.float32).compute() + return da.sum(cell_costs * self._lens).astype(np.float32).compute() @property def end_lat(self): From 31171277fa4f378cd181e5457a594dba5c15017c Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:54:15 -0700 Subject: [PATCH 37/61] Fix docs --- revrt/routing/utilities.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/revrt/routing/utilities.py b/revrt/routing/utilities.py index 2511ae7f..6101fbea 100644 --- a/revrt/routing/utilities.py +++ b/revrt/routing/utilities.py @@ -22,11 +22,11 @@ def map_to_costs(route_points, crs, transform, shape): route_points : pandas.DataFrame Route definitions table with at least `start_lat`, `start_lon`, `end_lat`, and `end_lon` coordinate columns. - crs : str or pyproj.CRS + crs : str or pyproj.crs.CRS Coordinate reference system for the cost raster. transform : affine.Affine Rasterio affine transform giving pixel origin and resolution. - shape : tuple[int, int] + shape : tuple Raster height and width for bounds checking. Returns From b12b5e6c555eaa71779a4cea804286f8d18bfabe Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:54:33 -0700 Subject: [PATCH 38/61] No extra log --- revrt/routing/cli.py | 1 - 1 file changed, 1 deletion(-) diff --git a/revrt/routing/cli.py b/revrt/routing/cli.py index e9d8476d..3af8d990 100644 --- a/revrt/routing/cli.py +++ b/revrt/routing/cli.py @@ -230,7 +230,6 @@ def compute_lcp_routes( # noqa: PLR0913, PLR0917 logger.debug("Transmission config input: %r", transmission_config) transmission_config = parse_config(config=transmission_config) - logger.debug("Transmission Config: %s", transmission_config) inds = _get_indices( route_table, From 73dd97401478bf48d99631979fff626585de99d3 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 22:58:03 -0700 Subject: [PATCH 39/61] Fix recursion error --- revrt/routing/point_to_many.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/revrt/routing/point_to_many.py b/revrt/routing/point_to_many.py index 6bb14922..7f4ec666 100644 --- a/revrt/routing/point_to_many.py +++ b/revrt/routing/point_to_many.py @@ -393,7 +393,7 @@ def __init__( self._routing_layers = routing_layers self._route = route self._optimized_objective = optimized_objective - self._lens = self._total_path_length = None + self.__lens = self._total_path_length = None self._by_layer_results = {} self._add_geom = add_geom self._attrs = attrs or {} @@ -406,9 +406,9 @@ def cell_size(self): @property def _lens(self): """array-like: Cached per-cell travel distances""" - if self._lens is None: + if self.__lens is None: self._compute_path_length() - return self._lens + return self.__lens @property def total_path_length(self): @@ -498,7 +498,7 @@ def build(self): def _compute_path_length(self): """Compute the total length and cell by cell length of LCP""" - self._lens, self._total_path_length = _compute_lens( + self.__lens, self._total_path_length = _compute_lens( self._route, self.cell_size ) From fbec8933bde172592e0466893251733684341d21 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 23:19:43 -0700 Subject: [PATCH 40/61] Add tests --- .../routing/test_routing_point_to_many.py | 403 +++++++++++++++++- 1 file changed, 402 insertions(+), 1 deletion(-) diff --git a/tests/python/unit/routing/test_routing_point_to_many.py b/tests/python/unit/routing/test_routing_point_to_many.py index a0c54329..de138d27 100644 --- a/tests/python/unit/routing/test_routing_point_to_many.py +++ b/tests/python/unit/routing/test_routing_point_to_many.py @@ -7,8 +7,17 @@ import xarray as xr from rasterio.transform import from_origin +from revrt.routing import point_to_many from revrt.utilities import LayeredFile -from revrt.routing.point_to_many import find_all_routes, RoutingScenario +from revrt.routing.point_to_many import ( + find_all_routes, + LCP_AGG_COST_LAYER_NAME, + RouteResult, + RoutingLayers, + RoutingScenario, +) +from revrt.exceptions import revrtKeyError, revrtLeastCostPathNotFoundError +from revrt.warn import revrtWarning @pytest.fixture(scope="module") @@ -576,5 +585,397 @@ def test_all_endpoints_are_barriers_returns_no_route( assert not output +def test_routing_scenario_repr_contains_fields(sample_layered_data): + """RoutingScenario repr surfaces configured layer metadata""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + friction_layers=[{"layer_name": "layer_2"}], + cost_multiplier_layer="layer_3", + cost_multiplier_scalar=1.5, + ) + + representation = repr(scenario) + + assert "layer_1" in representation + assert "layer_2" in representation + assert "cost_multiplier_scalar: 1.5" in representation + + +def test_missing_cost_layer_raises_key_error(sample_layered_data): + """Missing layers surface a revrtKeyError during build""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "not_there"}], + ) + + with pytest.raises( + revrtKeyError, match="Did not find layer 'not_there' in cost file" + ): + find_all_routes( + scenario, + route_definitions=[ + ((1, 1), [(1, 2)], {}), + ], + save_paths=False, + ) + + +def test_cost_multiplier_layer_and_scalar_applied(sample_layered_data): + """Cost multipliers scale base costs before routing aggregation""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + cost_multiplier_layer="layer_3", + cost_multiplier_scalar=2.0, + ) + + routing_layers = RoutingLayers(scenario).build() + try: + cost_val = routing_layers.cost.isel(y=1, x=1).compute().item() + layer_one = ( + routing_layers._layer_fh["layer_1"] + .isel(band=0, y=1, x=1) + .compute() + .item() + ) + layer_three = ( + routing_layers._layer_fh["layer_3"] + .isel(band=0, y=1, x=1) + .compute() + .item() + ) + expected = layer_one * layer_three * scenario.cost_multiplier_scalar + + assert cost_val == pytest.approx(expected) + finally: + routing_layers.close() + + +def test_length_invariant_layer_costs_ignore_path_length( + sample_layered_data, +): + """Length invariant cost layers ignore per-cell distances""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[ + {"layer_name": "layer_1"}, + {"layer_name": "layer_2", "is_invariant": True}, + ], + ) + + routing_layers = RoutingLayers(scenario).build() + try: + route = [(1, 1), (1, 2)] + result = RouteResult( + routing_layers, + route, + optimized_objective=0.0, + ).build() + + layer_two = ( + routing_layers._layer_fh["layer_2"].isel(band=0).compute().values + ) + expected = sum(layer_two[row, col] for row, col in route) + + assert result["layer_2_cost"] == pytest.approx(expected) + finally: + routing_layers.close() + + +def test_soft_barrier_setting_controls_barrier_value(sample_layered_data): + """Soft barriers convert impassable cells to large positive costs""" + + hard_scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + use_hard_barrier=True, + ) + hard_layers = RoutingLayers(hard_scenario).build() + try: + hard_value = ( + hard_layers.final_routing_layer.isel(y=0, x=3).compute().item() + ) + finally: + hard_layers.close() + + soft_scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + use_hard_barrier=False, + ) + soft_layers = RoutingLayers(soft_scenario).build() + try: + soft_value = ( + soft_layers.final_routing_layer.isel(y=0, x=3).compute().item() + ) + assert hard_value == -1 + assert soft_value > 0 + assert soft_value > abs(hard_value) + finally: + soft_layers.close() + + +def test_tracked_layers_invalid_configs_warn( + sample_layered_data, assert_message_was_logged +): + """Tracked layer config issues emit revrtWarning messages""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + tracked_layers={ + "layer_1": "does_not_exist", + "missing_layer": "mean", + }, + ) + + with pytest.warns(revrtWarning) as warning_record: + routing_layers = RoutingLayers(scenario).build() + + assert_message_was_logged("Did not find layer", "WARNING") + assert_message_was_logged("Did not find method", "WARNING") + + try: + assert len(warning_record) == 2 + finally: + routing_layers.close() + + +def test_friction_layers_and_lcp_agg_costs(sample_layered_data): + """Friction layers may include cost stack and tracked layer toggles""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[ + {"layer_name": "layer_1", "include_in_report": False}, + ], + friction_layers=[ + { + "layer_name": "layer_2", + "multiplier_scalar": 0.5, + "include_in_report": True, + }, + { + "layer_name": LCP_AGG_COST_LAYER_NAME, + "multiplier_scalar": 0.1, + }, + ], + ) + + routing_layers = RoutingLayers(scenario).build() + try: + tracked_names = {layer.name for layer in routing_layers.tracked_layers} + assert "layer_1" not in tracked_names + assert "layer_2" in tracked_names + + base_value = routing_layers.cost.isel(y=1, x=1).compute().item() + final_value = ( + routing_layers.final_routing_layer.isel(y=1, x=1).compute().item() + ) + + assert final_value > base_value + finally: + routing_layers.close() + + +def test_route_result_build_warns_on_attr_mismatch( + sample_layered_data, assert_message_was_logged +): + """RouteResult build warns when provided attrs contradict results""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + routing_layers = RoutingLayers(scenario).build() + try: + route = [(1, 1), (1, 2)] + with pytest.warns(revrtWarning): + result = RouteResult( + routing_layers, + route, + optimized_objective=0.0, + add_geom=True, + attrs={"start_row": 0}, + ).build() + + assert_message_was_logged("does not match", "WARNING") + + assert result["geometry"].geom_type == "LineString" + assert result["start_row"] == 0 + finally: + routing_layers.close() + + +def test_route_result_geom_returns_point_for_single_cell(sample_layered_data): + """RouteResult.geom returns a Point geometry for single-cell routes""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + routing_layers = RoutingLayers(scenario).build() + try: + route = [(1, 1)] + result = RouteResult( + routing_layers, + route, + optimized_objective=0.0, + ) + + assert result.geom.geom_type == "Point" + finally: + routing_layers.close() + + +def test_characterized_layer_length_metric_uses_positive_mask( + sample_layered_data, +): + """CharacterizedLayer uses positive-value mask when summing lengths""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + routing_layers = RoutingLayers(scenario).build() + try: + layer = next( + tracked + for tracked in routing_layers.tracked_layers + if tracked.name == "layer_1" + ) + route = [(1, 1), (1, 2), (2, 3)] + metrics = layer.compute(route, abs(routing_layers.transform.a)) + + assert metrics["layer_1_dist_km"] >= 0 + finally: + routing_layers.close() + + +def test_route_result_cached_properties_reuse_computed_values( + sample_layered_data, +): + """RouteResult caches per-route lengths after first computation""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + routing_layers = RoutingLayers(scenario).build() + try: + route = [(1, 1), (1, 2), (2, 3)] + result = RouteResult( + routing_layers, + route, + optimized_objective=0.0, + ) + + first_length = result.total_path_length + assert isinstance(first_length, float) + second_length = result.total_path_length + assert second_length == first_length + + first_lens = result._lens + assert np.allclose(result._lens, first_lens) + finally: + routing_layers.close() + + +def test_route_result_cost_property_returns_value(sample_layered_data): + """RouteResult.cost multiplies cell costs by cached travel lengths""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + routing_layers = RoutingLayers(scenario).build() + try: + route = [(1, 1), (1, 2), (2, 3)] + result = RouteResult( + routing_layers, + route, + optimized_objective=0.0, + ) + + assert result.cost > 0 + finally: + routing_layers.close() + + +def test_characterized_layer_total_length_computation(sample_layered_data): + """CharacterizedLayer computes length-weighted costs for eager data""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + routing_layers = RoutingLayers(scenario, chunks=None).build() + try: + layer = next( + tracked + for tracked in routing_layers.tracked_layers + if tracked.name == "layer_1" + ) + route = [(1, 1), (1, 2), (2, 3)] + metrics = layer.compute(route, abs(routing_layers.transform.a)) + + assert metrics["layer_1_cost"] > 0 + assert metrics["layer_1_dist_km"] >= 0 + finally: + routing_layers.close() + + +def test_find_paths_exception_yields_no_routes( + sample_layered_data, monkeypatch +): + """find_paths failures propagate as revrtLeastCostPathNotFoundError""" + + scenario = RoutingScenario( + cost_fpath=sample_layered_data, + cost_layers=[{"layer_name": "layer_1"}], + ) + + def boom(**_): + msg = "boom" + raise RuntimeError(msg) + + monkeypatch.setattr("revrt.routing.point_to_many.find_paths", boom) + + output = find_all_routes( + scenario, + route_definitions=[ + ((1, 1), [(1, 2)], {}), + ], + save_paths=False, + ) + + assert output == [] + + routing_layers = RoutingLayers(scenario).build() + try: + with pytest.raises( + revrtLeastCostPathNotFoundError, match="Unable to find path" + ): + point_to_many._compute_valid_path( + scenario, + routing_layers, + (1, 1), + [(1, 2)], + ) + finally: + routing_layers.close() + + if __name__ == "__main__": pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) From 7cd8bd40984e3cb0d95698bf583327d2d8ade258 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 23:44:19 -0700 Subject: [PATCH 41/61] MVP routing CLI tests --- tests/python/unit/routing/test_routing_cli.py | 519 ++++++++++++++++++ 1 file changed, 519 insertions(+) create mode 100644 tests/python/unit/routing/test_routing_cli.py diff --git a/tests/python/unit/routing/test_routing_cli.py b/tests/python/unit/routing/test_routing_cli.py new file mode 100644 index 00000000..a311e5a9 --- /dev/null +++ b/tests/python/unit/routing/test_routing_cli.py @@ -0,0 +1,519 @@ +"""reVRt routing CLI unit tests""" + +from pathlib import Path +import math + +import numpy as np +import pandas as pd +import pytest +import xarray as xr +import geopandas as gpd +from shapely.geometry import Point +from rasterio.transform import from_origin + +from revrt.utilities import LayeredFile +from revrt.routing.utilities import map_to_costs +from revrt.routing.cli import ( + compute_lcp_routes, + _run_lcp, + _collect_existing_routes, + _update_multipliers, + _route_points_subset, + _paths_to_compute, + _split_routes, + _MILLION_USD_PER_MILE_TO_USD_PER_PIXEL, +) + + +@pytest.fixture(scope="module") +def sample_layered_data(tmp_path_factory): + """Create layered routing data mimicking point_to_many tests""" + + data_dir = tmp_path_factory.mktemp("routing_cli_data") + + layered_fp = data_dir / "test_layered.zarr" + layer_file = LayeredFile(layered_fp) + + height, width = (7, 8) + cell_size = 1.0 + x0, y0 = 0.0, float(height) + transform = from_origin(x0, y0, cell_size, cell_size) + x_coords = ( + x0 + np.arange(width, dtype=np.float32) * cell_size + cell_size / 2 + ) + y_coords = ( + y0 - np.arange(height, dtype=np.float32) * cell_size - cell_size / 2 + ) + + layer_values = [ + np.array( + [ + [ + [7, 7, 8, 0, 9, 9, 9, 0], + [8, 1, 2, 2, 9, 9, 9, 0], + [9, 1, 3, 3, 9, 1, 2, 3], + [9, 1, 2, 1, 9, 1, 9, 0], + [9, 9, 9, 1, 9, 1, 9, 0], + [9, 9, 9, 1, 1, 1, 9, 0], + [9, 9, 9, 9, 9, 9, 9, 0], + ] + ], + dtype=np.float32, + ), + np.array( + [ + [ + [8, 7, 6, 5, 5, 6, 7, 9], + [7, 1, 1, 2, 3, 3, 2, 8], + [6, 2, 9, 6, 5, 2, 1, 7], + [7, 3, 8, 1, 2, 3, 2, 6], + [8, 4, 7, 2, 8, 4, 3, 5], + [9, 5, 6, 3, 4, 4, 3, 4], + [9, 6, 7, 4, 5, 5, 4, 3], + ] + ], + dtype=np.float32, + ), + np.array( + [ + [ + [6, 6, 6, 6, 6, 7, 8, 9], + [5, 2, 2, 3, 4, 5, 6, 8], + [4, 3, 7, 7, 6, 4, 5, 7], + [5, 4, 6, 2, 3, 4, 4, 6], + [6, 5, 5, 3, 7, 5, 5, 5], + [7, 6, 6, 4, 5, 5, 4, 4], + [8, 7, 7, 5, 6, 5, 4, 3], + ] + ], + dtype=np.float32, + ), + ] + + for ind, routing_layer in enumerate(layer_values, start=1): + da = xr.DataArray( + routing_layer, + dims=("band", "y", "x"), + coords={"y": y_coords, "x": x_coords}, + ) + da = da.rio.write_crs("EPSG:4326") + da = da.rio.write_transform(transform) + + geotiff_fp = data_dir / f"layer_{ind}.tif" + da.rio.to_raster(geotiff_fp, driver="GTiff") + + layer_file.write_geotiff_to_file( + geotiff_fp, f"layer_{ind}", overwrite=True + ) + + return layered_fp + + +def _build_route_table(layered_fp, rows_cols): + """Helper to construct route tables with CRS-aligned coordinates""" + + with xr.open_dataset(layered_fp, consolidated=False, engine="zarr") as ds: + latitudes = ds["latitude"].to_numpy() + longitudes = ds["longitude"].to_numpy() + + records = [] + for idx, (start, end) in enumerate(rows_cols): + s_row, s_col = start + e_row, e_col = end + records.append( + { + "route_id": f"route_{idx}", + "start_lat": float(latitudes[s_row, s_col]), + "start_lon": float(longitudes[s_row, s_col]), + "end_lat": float(latitudes[e_row, e_col]), + "end_lon": float(longitudes[e_row, e_col]), + "voltage": 138, + "polarity": "ac", + } + ) + + return pd.DataFrame.from_records(records) + + +def test_compute_lcp_routes_generates_csv( + sample_layered_data, tmp_path, monkeypatch +): + """compute_lcp_routes should map points and persist CSV outputs""" + + routes = _build_route_table( + sample_layered_data, + rows_cols=[((1, 1), (2, 3)), ((0, 0), (3, 4))], + ) + route_table_fp = tmp_path / "route_table.csv" + routes.to_csv(route_table_fp, index=False) + + monkeypatch.setattr( + "revrt.routing.cli._route_points_subset", + lambda *_args, **_kwargs: routes, + ) + monkeypatch.setattr( + "revrt.routing.cli.parse_config", lambda config=None: config + ) + + captured_scenarios = [] + + def fake_find_all_routes(scenario, route_definitions, save_paths): + captured_scenarios.append(scenario) + outputs = [] + for start, end_points, info in route_definitions: + end_row, end_col = end_points[0] + out = info.copy() + out.update( + { + "start_row": start[0], + "start_col": start[1], + "end_row": end_row, + "end_col": end_col, + "cost": sum( + layer.get("multiplier_scalar", 1) + for layer in scenario.cost_layers + ), + "length_km": 0.001, + "optimized_objective": 0.001, + } + ) + outputs.append(out) + return outputs + + monkeypatch.setattr( + "revrt.routing.cli.find_all_routes", fake_find_all_routes + ) + + out_dir = tmp_path / "routing_outputs" + + cost_layers = [ + { + "layer_name": "layer_1", + "apply_row_mult": True, + "multiplier_scalar": 1, + }, + { + "layer_name": "layer_2", + "apply_polarity_mult": True, + "multiplier_scalar": 1, + }, + ] + + transmission_config = { + "row_width": {"138": 1.5}, + "voltage_polarity_mult": {"138": {"ac": 2e-5}}, + } + + result_fp = compute_lcp_routes( + cost_fpath=sample_layered_data, + route_table=route_table_fp, + cost_layers=cost_layers, + out_dir=out_dir, + job_name="run", + transmission_config=transmission_config, + save_paths=False, + ) + + output_path = Path(result_fp) + assert output_path.exists() + + df = pd.read_csv(output_path) + df = df[df["route_id"] != "route_id"].reset_index(drop=True) + df["cost"] = pd.to_numeric(df["cost"]) + assert len(df) == len(routes) + expected_scalar = 1.5 + 2e-5 * _MILLION_USD_PER_MILE_TO_USD_PER_PIXEL + assert pytest.approx(expected_scalar) == df["cost"].iloc[0] + + updated_layers = captured_scenarios[0].cost_layers + assert updated_layers[0]["multiplier_scalar"] == pytest.approx(1.5) + assert updated_layers[1]["multiplier_scalar"] == pytest.approx( + 2e-5 * _MILLION_USD_PER_MILE_TO_USD_PER_PIXEL + ) + + +def test_compute_lcp_routes_returns_none_on_empty_indices( + sample_layered_data, tmp_path, monkeypatch +): + """Ensure compute_lcp_routes short-circuits when route points are empty""" + + route_table_fp = tmp_path / "routes.csv" + _build_route_table(sample_layered_data, [((1, 1), (2, 2))]).to_csv( + route_table_fp, index=False + ) + + monkeypatch.setattr( + "revrt.routing.cli._route_points_subset", + lambda *_, **__: pd.DataFrame([]), + ) + monkeypatch.setattr( + "revrt.routing.cli.parse_config", lambda config=None: config + ) + + result = compute_lcp_routes( + cost_fpath=sample_layered_data, + route_table=route_table_fp, + cost_layers=[{"layer_name": "layer_1"}], + out_dir=tmp_path, + job_name="no_routes", + ) + + assert result is None + assert not (tmp_path / "no_routes.csv").exists() + + +def test_run_lcp_with_save_paths_filters_existing_routes( + sample_layered_data, tmp_path, monkeypatch +): + """_run_lcp should skip already processed routes and append geometries""" + + routes = _build_route_table( + sample_layered_data, + rows_cols=[((1, 1), (2, 2)), ((2, 2), (4, 4))], + ) + + with xr.open_dataset( + sample_layered_data, consolidated=False, engine="zarr" + ) as ds: + mapped_routes = map_to_costs( + routes.copy(), ds.rio.crs, ds.rio.transform(), ds.rio.shape + ) + + existing_tuple = ( + int(mapped_routes.iloc[0]["start_row"]), + int(mapped_routes.iloc[0]["start_col"]), + int(mapped_routes.iloc[0]["end_row"]), + int(mapped_routes.iloc[0]["end_col"]), + routes.iloc[0]["polarity"], + str(routes.iloc[0]["voltage"]), + ) + + monkeypatch.setattr( + "revrt.routing.cli._collect_existing_routes", + lambda _: {existing_tuple}, + ) + + captured_scenarios = [] + + def fake_find_all_routes(scenario, route_definitions, save_paths): + captured_scenarios.append(scenario) + assert len(route_definitions) == 1 + (start_row, start_col), end_points, info = route_definitions[0] + end_row, end_col = end_points[0] + return [ + { + "route_id": info["route_id"], + "start_row": start_row, + "start_col": start_col, + "end_row": end_row, + "end_col": end_col, + "geometry": Point(0, 0), + "cost": math.pi, + "length_km": 0.002, + "optimized_objective": 0.002, + } + ] + + monkeypatch.setattr( + "revrt.routing.cli.find_all_routes", fake_find_all_routes + ) + + saved_calls = [] + + def fake_to_file(self, path, driver=None, mode=None, **_kwargs): + saved_calls.append((path, driver, mode, self.copy(deep=True))) + + monkeypatch.setattr( + "revrt.routing.cli.gpd.GeoDataFrame.to_file", fake_to_file + ) + + out_fp = tmp_path / "routes.gpkg" + + _run_lcp( + cost_fpath=sample_layered_data, + route_points=routes, + cost_layers=[{"layer_name": "layer_1"}], + out_fp=out_fp, + transmission_config={ + "row_width": {"138": 1.0}, + "voltage_polarity_mult": {"138": {"ac": 1.0}}, + }, + cost_multiplier_scalar=1, + friction_layers=[{"layer_name": "layer_2", "apply_row_mult": True}], + tracked_layers={"layer_3": "max"}, + use_hard_barrier=True, + ) + + assert len(saved_calls) == 1 + saved_path, driver, mode, saved_gdf = saved_calls[0] + assert saved_path == out_fp + assert driver == "GPKG" + assert mode == "a" + assert saved_gdf["route_id"].iloc[0] == routes.iloc[1]["route_id"] + + updated_layers = captured_scenarios[0].cost_layers + assert updated_layers[0]["layer_name"] == "layer_1" + + +def test_run_lcp_returns_immediately_when_no_routes(tmp_path): + """_run_lcp should exit early when route_points is empty""" + + _run_lcp( + cost_fpath="unused", # cost file is ignored in this branch + route_points=pd.DataFrame(), + cost_layers=[], + out_fp=tmp_path / "unused.csv", + ) + + +def test_collect_existing_routes_csv(tmp_path): + """_collect_existing_routes should read CSV outputs""" + + data = pd.DataFrame( + [ + { + "start_row": 0, + "start_col": 1, + "end_row": 2, + "end_col": 3, + "polarity": "ac", + "voltage": "230", + } + ] + ) + csv_fp = tmp_path / "routes.csv" + data.to_csv(csv_fp, index=False) + + result = _collect_existing_routes(csv_fp) + assert result == {(0, 1, 2, 3, "ac", "230")} + + +def test_collect_existing_routes_gpkg(monkeypatch, tmp_path): + """_collect_existing_routes should support GeoPackage outputs""" + + gpkg_fp = tmp_path / "routes.gpkg" + gpkg_fp.touch() + + gdf = gpd.GeoDataFrame( + { + "start_row": [1], + "start_col": [2], + "end_row": [3], + "end_col": [4], + "polarity": ["unknown"], + "voltage": ["unknown"], + "geometry": [Point(0, 0)], + } + ) + + monkeypatch.setattr("revrt.routing.cli.gpd.read_file", lambda _: gdf) + + result = _collect_existing_routes(gpkg_fp) + assert result == {(1, 2, 3, 4, "unknown", "unknown")} + + +def test_collect_existing_routes_when_missing(tmp_path): + """Missing outputs should result in an empty existing route set""" + + assert _collect_existing_routes(None) == set() + assert _collect_existing_routes(tmp_path / "missing.csv") == set() + + +def test_route_points_subset_with_chunking(monkeypatch): + """_route_points_subset should slice sorted features by chunk""" + + features = pd.DataFrame( + { + "start_lat": [5.0, 1.0, 3.0, 7.0], + "start_lon": [0.0, 1.0, 2.0, 3.0], + } + ) + + monkeypatch.setattr("revrt.routing.cli.pd.read_csv", lambda _: features) + + first_chunk = _route_points_subset( + "dummy", ["start_lat", "start_lon"], (0, 2) + ) + assert first_chunk["start_lat"].tolist() == [1.0, 3.0] + assert first_chunk["start_lon"].tolist() == [1.0, 2.0] + + second_chunk = _route_points_subset( + "dummy", ["start_lat", "start_lon"], (1, 2) + ) + assert second_chunk["start_lat"].tolist() == [5.0, 7.0] + assert second_chunk["start_lon"].tolist() == [0.0, 3.0] + + +def test_paths_to_compute_inserts_missing_columns(tmp_path): + """_paths_to_compute should back-fill missing polarity/voltage columns""" + + route_points = pd.DataFrame( + { + "start_row": [0], + "start_col": [1], + "end_row": [2], + "end_col": [3], + } + ) + + groups = list(_paths_to_compute(route_points, tmp_path / "not_there.csv")) + assert groups + polarity, voltage, grouped_routes = groups[0] + assert polarity == "unknown" + assert voltage == "unknown" + assert grouped_routes.iloc[0]["start_row"] == 0 + + +def test_split_routes_handles_local_and_cluster(): + """_split_routes should configure chunking for local and cluster modes""" + + local_config = {"execution_control": {"option": "local", "nodes": 4}} + result_local = _split_routes(local_config) + assert result_local["_split_params"] == [(0, 1)] + assert result_local["execution_control"]["nodes"] == 4 + + cluster_config = {"execution_control": {"nodes": 3}} + result_cluster = _split_routes(cluster_config) + assert result_cluster["_split_params"] == [(0, 3), (1, 3), (2, 3)] + assert "nodes" not in result_cluster["execution_control"] + + +def test_update_multipliers_applies_row_and_polarity(): + """_update_multipliers should apply configured scalar adjustments""" + + layers = [ + { + "layer_name": "layer_1", + "multiplier_scalar": 2, + "apply_row_mult": True, + }, + {"layer_name": "layer_2", "apply_polarity_mult": True}, + ] + + transmission_config = { + "row_width": {"138": 1.5}, + "voltage_polarity_mult": {"138": {"ac": 0.5}}, + } + + updated = _update_multipliers( + layers, + polarity="ac", + voltage=138, + transmission_config=transmission_config, + ) + + # original input remains unchanged + assert layers[0]["apply_row_mult"] is True + assert updated[0]["multiplier_scalar"] == pytest.approx(3) + assert updated[1]["multiplier_scalar"] == pytest.approx( + 0.5 * _MILLION_USD_PER_MILE_TO_USD_PER_PIXEL + ) + + # Voltage marked as unknown should skip multiplier lookups + unchanged = _update_multipliers( + [{"layer_name": "layer_3"}], "dc", "unknown", transmission_config + ) + assert unchanged[0]["layer_name"] == "layer_3" + + +if __name__ == "__main__": + pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) From c78f4f632bb3127632e9d75778eff82a17d0a84f Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 23:45:11 -0700 Subject: [PATCH 42/61] Minor refactor --- revrt/routing/cli.py | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/revrt/routing/cli.py b/revrt/routing/cli.py index 3af8d990..c161f49c 100644 --- a/revrt/routing/cli.py +++ b/revrt/routing/cli.py @@ -231,21 +231,15 @@ def compute_lcp_routes( # noqa: PLR0913, PLR0917 transmission_config = parse_config(config=transmission_config) - inds = _get_indices( + route_points = _route_points_subset( route_table, sort_cols=["start_lat", "start_lon"], split_params=_split_params, ) - if len(inds) == 0: - logger.info("No indices to process: %s", inds) + if len(route_points) == 0: + logger.info("No routes to process!") return None - with contextlib.suppress(TypeError, UnicodeDecodeError): - route_points = pd.read_csv(route_table) - - if inds is not None: - route_points = route_points.loc[inds] - out_fp = ( out_dir / f"{job_name}.gpkg" if save_paths @@ -430,14 +424,17 @@ def _update_multipliers(layers, polarity, voltage, transmission_config): return output_layers -def _get_indices(features_fp, sort_cols, split_params): +def _route_points_subset(route_table, sort_cols, split_params): """Get indices of points that are sorted by location""" - features = gpd.read_file(features_fp) - features = features.sort_values(sort_cols) + + with contextlib.suppress(TypeError, UnicodeDecodeError): + route_points = pd.read_csv(route_table) + + route_points = route_points.sort_values(sort_cols).reset_index(drop=True) start_ind, n_chunks = split_params or (0, 1) - chunk_size = ceil(len(features) / n_chunks) - return features.index[ + chunk_size = ceil(len(route_points) / n_chunks) + return route_points.iloc[ start_ind * chunk_size : (start_ind + 1) * chunk_size ] From e614fa68bdf7c4999cb8a2902e9119a4a5bf498c Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 23:46:53 -0700 Subject: [PATCH 43/61] Minor formatting --- revrt/routing/cli.py | 1 - 1 file changed, 1 deletion(-) diff --git a/revrt/routing/cli.py b/revrt/routing/cli.py index c161f49c..4bf9e840 100644 --- a/revrt/routing/cli.py +++ b/revrt/routing/cli.py @@ -285,7 +285,6 @@ def _run_lcp( if route_points.empty: logger.info("Found no paths to compute!") - return with xr.open_dataset(cost_fpath, consolidated=False, engine="zarr") as ds: From 15cd004a369a05ec532163f27de96a25da003501 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Tue, 9 Dec 2025 23:55:51 -0700 Subject: [PATCH 44/61] update tests --- tests/python/unit/routing/test_routing_cli.py | 133 +++++++++--------- 1 file changed, 68 insertions(+), 65 deletions(-) diff --git a/tests/python/unit/routing/test_routing_cli.py b/tests/python/unit/routing/test_routing_cli.py index a311e5a9..690c6565 100644 --- a/tests/python/unit/routing/test_routing_cli.py +++ b/tests/python/unit/routing/test_routing_cli.py @@ -1,14 +1,13 @@ """reVRt routing CLI unit tests""" from pathlib import Path -import math import numpy as np import pandas as pd import pytest import xarray as xr import geopandas as gpd -from shapely.geometry import Point +from shapely.geometry import Point, LineString from rasterio.transform import from_origin from revrt.utilities import LayeredFile @@ -155,35 +154,6 @@ def test_compute_lcp_routes_generates_csv( "revrt.routing.cli.parse_config", lambda config=None: config ) - captured_scenarios = [] - - def fake_find_all_routes(scenario, route_definitions, save_paths): - captured_scenarios.append(scenario) - outputs = [] - for start, end_points, info in route_definitions: - end_row, end_col = end_points[0] - out = info.copy() - out.update( - { - "start_row": start[0], - "start_col": start[1], - "end_row": end_row, - "end_col": end_col, - "cost": sum( - layer.get("multiplier_scalar", 1) - for layer in scenario.cost_layers - ), - "length_km": 0.001, - "optimized_objective": 0.001, - } - ) - outputs.append(out) - return outputs - - monkeypatch.setattr( - "revrt.routing.cli.find_all_routes", fake_find_all_routes - ) - out_dir = tmp_path / "routing_outputs" cost_layers = [ @@ -219,15 +189,57 @@ def fake_find_all_routes(scenario, route_definitions, save_paths): df = pd.read_csv(output_path) df = df[df["route_id"] != "route_id"].reset_index(drop=True) - df["cost"] = pd.to_numeric(df["cost"]) + df = df.drop(columns=[c for c in df.columns if c.startswith("Unnamed")]) + + numeric_cols = [ + "cost", + "length_km", + "optimized_objective", + "layer_1_cost", + "layer_1_dist_km", + "layer_2_cost", + "layer_2_dist_km", + ] + for col in numeric_cols: + if col in df.columns: + df[col] = pd.to_numeric(df[col]) + assert len(df) == len(routes) - expected_scalar = 1.5 + 2e-5 * _MILLION_USD_PER_MILE_TO_USD_PER_PIXEL - assert pytest.approx(expected_scalar) == df["cost"].iloc[0] + assert set(df["route_id"]) == set(routes["route_id"]) - updated_layers = captured_scenarios[0].cost_layers - assert updated_layers[0]["multiplier_scalar"] == pytest.approx(1.5) - assert updated_layers[1]["multiplier_scalar"] == pytest.approx( - 2e-5 * _MILLION_USD_PER_MILE_TO_USD_PER_PIXEL + with xr.open_dataset( + sample_layered_data, consolidated=False, engine="zarr" + ) as ds: + mapped_routes = map_to_costs( + routes.copy(), ds.rio.crs, ds.rio.transform(), ds.rio.shape + ) + + merged = df.merge( + mapped_routes[ + [ + "route_id", + "start_row", + "start_col", + "end_row", + "end_col", + ] + ], + on="route_id", + how="left", + suffixes=("", "_expected"), + ) + + for col in ["start_row", "start_col", "end_row", "end_col"]: + assert ( + merged[col].astype(int) == merged[f"{col}_expected"].astype(int) + ).all() + + assert np.allclose( + merged["cost"], merged["layer_1_cost"] + merged["layer_2_cost"] + ) + assert np.all(merged["length_km"] > 0) + assert np.allclose( + merged["cost"], merged["optimized_objective"], rtol=1e-5 ) @@ -292,31 +304,6 @@ def test_run_lcp_with_save_paths_filters_existing_routes( lambda _: {existing_tuple}, ) - captured_scenarios = [] - - def fake_find_all_routes(scenario, route_definitions, save_paths): - captured_scenarios.append(scenario) - assert len(route_definitions) == 1 - (start_row, start_col), end_points, info = route_definitions[0] - end_row, end_col = end_points[0] - return [ - { - "route_id": info["route_id"], - "start_row": start_row, - "start_col": start_col, - "end_row": end_row, - "end_col": end_col, - "geometry": Point(0, 0), - "cost": math.pi, - "length_km": 0.002, - "optimized_objective": 0.002, - } - ] - - monkeypatch.setattr( - "revrt.routing.cli.find_all_routes", fake_find_all_routes - ) - saved_calls = [] def fake_to_file(self, path, driver=None, mode=None, **_kwargs): @@ -348,10 +335,26 @@ def fake_to_file(self, path, driver=None, mode=None, **_kwargs): assert saved_path == out_fp assert driver == "GPKG" assert mode == "a" + assert len(saved_gdf) == 1 assert saved_gdf["route_id"].iloc[0] == routes.iloc[1]["route_id"] - updated_layers = captured_scenarios[0].cost_layers - assert updated_layers[0]["layer_name"] == "layer_1" + expected = mapped_routes.iloc[1] + assert int(saved_gdf["start_row"].iloc[0]) == int(expected["start_row"]) + assert int(saved_gdf["start_col"].iloc[0]) == int(expected["start_col"]) + assert int(saved_gdf["end_row"].iloc[0]) == int(expected["end_row"]) + assert int(saved_gdf["end_col"].iloc[0]) == int(expected["end_col"]) + + cost_val = float(saved_gdf["cost"].iloc[0]) + objective_val = float(saved_gdf["optimized_objective"].iloc[0]) + length_val = float(saved_gdf["length_km"].iloc[0]) + + assert cost_val > 0 + assert length_val > 0 + assert objective_val == pytest.approx(cost_val, rel=1e-5) + + geom = saved_gdf.geometry.iloc[0] + assert isinstance(geom, LineString) + assert len(geom.coords) >= 2 def test_run_lcp_returns_immediately_when_no_routes(tmp_path): From 024e837500222c74954be11c43229b641d4f4e56 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 00:10:08 -0700 Subject: [PATCH 45/61] Update tests --- tests/python/unit/routing/test_routing_cli.py | 41 +++++-------------- 1 file changed, 11 insertions(+), 30 deletions(-) diff --git a/tests/python/unit/routing/test_routing_cli.py b/tests/python/unit/routing/test_routing_cli.py index 690c6565..279fe758 100644 --- a/tests/python/unit/routing/test_routing_cli.py +++ b/tests/python/unit/routing/test_routing_cli.py @@ -134,9 +134,7 @@ def _build_route_table(layered_fp, rows_cols): return pd.DataFrame.from_records(records) -def test_compute_lcp_routes_generates_csv( - sample_layered_data, tmp_path, monkeypatch -): +def test_compute_lcp_routes_generates_csv(sample_layered_data, tmp_path): """compute_lcp_routes should map points and persist CSV outputs""" routes = _build_route_table( @@ -146,16 +144,7 @@ def test_compute_lcp_routes_generates_csv( route_table_fp = tmp_path / "route_table.csv" routes.to_csv(route_table_fp, index=False) - monkeypatch.setattr( - "revrt.routing.cli._route_points_subset", - lambda *_args, **_kwargs: routes, - ) - monkeypatch.setattr( - "revrt.routing.cli.parse_config", lambda config=None: config - ) - out_dir = tmp_path / "routing_outputs" - cost_layers = [ { "layer_name": "layer_1", @@ -182,6 +171,7 @@ def test_compute_lcp_routes_generates_csv( job_name="run", transmission_config=transmission_config, save_paths=False, + _split_params=(0, 1), ) output_path = Path(result_fp) @@ -244,7 +234,7 @@ def test_compute_lcp_routes_generates_csv( def test_compute_lcp_routes_returns_none_on_empty_indices( - sample_layered_data, tmp_path, monkeypatch + sample_layered_data, tmp_path ): """Ensure compute_lcp_routes short-circuits when route points are empty""" @@ -253,20 +243,13 @@ def test_compute_lcp_routes_returns_none_on_empty_indices( route_table_fp, index=False ) - monkeypatch.setattr( - "revrt.routing.cli._route_points_subset", - lambda *_, **__: pd.DataFrame([]), - ) - monkeypatch.setattr( - "revrt.routing.cli.parse_config", lambda config=None: config - ) - result = compute_lcp_routes( cost_fpath=sample_layered_data, route_table=route_table_fp, cost_layers=[{"layer_name": "layer_1"}], out_dir=tmp_path, job_name="no_routes", + _split_params=(1000, 1), ) assert result is None @@ -390,12 +373,10 @@ def test_collect_existing_routes_csv(tmp_path): assert result == {(0, 1, 2, 3, "ac", "230")} -def test_collect_existing_routes_gpkg(monkeypatch, tmp_path): +def test_collect_existing_routes_gpkg(tmp_path): """_collect_existing_routes should support GeoPackage outputs""" gpkg_fp = tmp_path / "routes.gpkg" - gpkg_fp.touch() - gdf = gpd.GeoDataFrame( { "start_row": [1], @@ -407,8 +388,7 @@ def test_collect_existing_routes_gpkg(monkeypatch, tmp_path): "geometry": [Point(0, 0)], } ) - - monkeypatch.setattr("revrt.routing.cli.gpd.read_file", lambda _: gdf) + gdf.to_file(gpkg_fp, driver="GPKG") result = _collect_existing_routes(gpkg_fp) assert result == {(1, 2, 3, 4, "unknown", "unknown")} @@ -421,9 +401,10 @@ def test_collect_existing_routes_when_missing(tmp_path): assert _collect_existing_routes(tmp_path / "missing.csv") == set() -def test_route_points_subset_with_chunking(monkeypatch): +def test_route_points_subset_with_chunking(tmp_path): """_route_points_subset should slice sorted features by chunk""" + test_fp = tmp_path / "features.csv" features = pd.DataFrame( { "start_lat": [5.0, 1.0, 3.0, 7.0], @@ -431,16 +412,16 @@ def test_route_points_subset_with_chunking(monkeypatch): } ) - monkeypatch.setattr("revrt.routing.cli.pd.read_csv", lambda _: features) + features.to_csv(test_fp, index=False) first_chunk = _route_points_subset( - "dummy", ["start_lat", "start_lon"], (0, 2) + test_fp, ["start_lat", "start_lon"], (0, 2) ) assert first_chunk["start_lat"].tolist() == [1.0, 3.0] assert first_chunk["start_lon"].tolist() == [1.0, 2.0] second_chunk = _route_points_subset( - "dummy", ["start_lat", "start_lon"], (1, 2) + test_fp, ["start_lat", "start_lon"], (1, 2) ) assert second_chunk["start_lat"].tolist() == [5.0, 7.0] assert second_chunk["start_lon"].tolist() == [0.0, 3.0] From 45e5d393e317a1f6fd0c555fb45e389d2944e916 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 00:18:43 -0700 Subject: [PATCH 46/61] Update lockfile --- pixi.lock | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/pixi.lock b/pixi.lock index dc8a008e..c968dddb 100644 --- a/pixi.lock +++ b/pixi.lock @@ -2345,6 +2345,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/httpcore-1.0.9-pyh29332c3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/httpx-0.28.1-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hyperframe-6.1.0-pyhd8ed1ab_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/hypothesis-6.148.7-pyha770c72_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/icu-75.1-he02047a_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/id-1.5.0-pyh29332c3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/idna-3.11-pyhd8ed1ab_0.conda @@ -2841,6 +2842,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/httpcore-1.0.9-pyh29332c3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/httpx-0.28.1-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hyperframe-6.1.0-pyhd8ed1ab_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/hypothesis-6.148.7-pyha770c72_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/icu-75.1-h120a0e1_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/id-1.5.0-pyh29332c3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/idna-3.11-pyhd8ed1ab_0.conda @@ -3282,6 +3284,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/httpcore-1.0.9-pyh29332c3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/httpx-0.28.1-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hyperframe-6.1.0-pyhd8ed1ab_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/hypothesis-6.148.7-pyha770c72_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/icu-75.1-hfee45f7_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/id-1.5.0-pyh29332c3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/idna-3.11-pyhd8ed1ab_0.conda @@ -3726,6 +3729,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/httpcore-1.0.9-pyh29332c3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/httpx-0.28.1-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hyperframe-6.1.0-pyhd8ed1ab_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/hypothesis-6.148.7-pyha770c72_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/icu-75.1-he0c23c2_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/id-1.5.0-pyh29332c3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/idna-3.11-pyhd8ed1ab_0.conda @@ -5241,6 +5245,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/h2-4.3.0-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hpack-4.1.0-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hyperframe-6.1.0-pyhd8ed1ab_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/hypothesis-6.148.7-pyha770c72_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/icu-75.1-he02047a_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/idna-3.11-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/linux-64/imagecodecs-2025.11.11-py313hf092b87_1.conda @@ -5558,6 +5563,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/h2-4.3.0-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hpack-4.1.0-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hyperframe-6.1.0-pyhd8ed1ab_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/hypothesis-6.148.7-pyha770c72_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/icu-75.1-h120a0e1_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/idna-3.11-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-64/imagecodecs-2025.11.11-py313hca676d3_1.conda @@ -5863,6 +5869,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/h2-4.3.0-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hpack-4.1.0-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hyperframe-6.1.0-pyhd8ed1ab_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/hypothesis-6.148.7-pyha770c72_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/icu-75.1-hfee45f7_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/idna-3.11-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/osx-arm64/imagecodecs-2025.11.11-py313h58bcf60_1.conda @@ -6161,6 +6168,7 @@ environments: - conda: https://conda.anaconda.org/conda-forge/noarch/h2-4.3.0-pyhcf101f3_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hpack-4.1.0-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/hyperframe-6.1.0-pyhd8ed1ab_0.conda + - conda: https://conda.anaconda.org/conda-forge/noarch/hypothesis-6.148.7-pyha770c72_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/icu-75.1-he0c23c2_0.conda - conda: https://conda.anaconda.org/conda-forge/noarch/idna-3.11-pyhd8ed1ab_0.conda - conda: https://conda.anaconda.org/conda-forge/win-64/imagecodecs-2025.11.11-py313h7cfe00e_1.conda @@ -10382,6 +10390,22 @@ packages: - pkg:pypi/hyperframe?source=hash-mapping size: 17397 timestamp: 1737618427549 +- conda: https://conda.anaconda.org/conda-forge/noarch/hypothesis-6.148.7-pyha770c72_0.conda + sha256: 43d83347c867a016c7529a0e934158ef8cc6d7df7dc4092dfa42f280cd0c2e52 + md5: 1058d2fd766a1df97d22f998a25f6542 + depends: + - attrs >=22.2.0 + - click >=7.0 + - exceptiongroup >=1.0.0 + - python >=3.10 + - setuptools + - sortedcontainers >=2.1.0,<3.0.0 + license: MPL-2.0 + license_family: MOZILLA + purls: + - pkg:pypi/hypothesis?source=hash-mapping + size: 382777 + timestamp: 1764903748476 - conda: https://conda.anaconda.org/conda-forge/linux-64/icu-75.1-he02047a_0.conda sha256: 71e750d509f5fa3421087ba88ef9a7b9be11c53174af3aa4d06aff4c18b38e8e md5: 8b189310083baabfb622af68fd9d3ae3 @@ -16889,7 +16913,7 @@ packages: - pypi: ./ name: nrel-revrt version: 0.2.1 - sha256: 8db001fdf27ec3ffa8b7a17346a3e9965d791d1cacdff0da92b63d0544a8a25a + sha256: 102f68ffc4003a9215315d66814d795bce44460e1cbaf43a66ca5ee2707d1932 requires_dist: - affine>=2.4.0,<3 - dask>=2025.4.1,<2026 @@ -16910,6 +16934,7 @@ packages: - pipreqs>=0.4.13,<0.5 ; extra == 'dev' - ruff>=0.14.0,<0.15 ; extra == 'dev' - ruff-lsp>=0.0.62,<0.0.63 ; extra == 'dev' + - hypothesis>=6.0,<7 ; extra == 'test' - pytest>=8.3.5,<9 ; extra == 'test' - pytest-cases>=3.8.6,<4 ; extra == 'test' - pytest-cov>=6.1.1,<7 ; extra == 'test' From 06dad47849b78190b5e4ee1439004dd7e23d9a81 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 12:11:39 -0700 Subject: [PATCH 47/61] MInor formatting --- tests/python/unit/test_exceptions.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/tests/python/unit/test_exceptions.py b/tests/python/unit/test_exceptions.py index 951a7e28..c4dc443c 100644 --- a/tests/python/unit/test_exceptions.py +++ b/tests/python/unit/test_exceptions.py @@ -25,7 +25,7 @@ def test_exceptions_log_error(caplog, assert_message_was_logged): - """Test that a raised exception logs message, if any.""" + """Test that a raised exception logs message, if any""" try: raise revrtError # noqa: TRY301 @@ -102,8 +102,10 @@ def test_exceptions_log_uncaught_error(assert_message_was_logged): (revrtValueError, [revrtError, revrtValueError, ValueError]), ], ) -def test_catching_error_by_type(raise_type, catch_types, assert_message_was_logged): - """Test that gaps exceptions are caught correctly.""" +def test_catching_error_by_type( + raise_type, catch_types, assert_message_was_logged +): + """Test that gaps exceptions are caught correctly""" for catch_type in catch_types: with pytest.raises(catch_type) as exc_info: raise raise_type(BASIC_ERROR_MESSAGE) From b53749ee6d4920ae243a4344afb04971b4d837ac Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 12:18:56 -0700 Subject: [PATCH 48/61] Update instructions --- .github/copilot-instructions.md | 53 +++++++++++++++++++++++++-------- 1 file changed, 40 insertions(+), 13 deletions(-) diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md index 8b1ae1d1..9738b93b 100644 --- a/.github/copilot-instructions.md +++ b/.github/copilot-instructions.md @@ -87,25 +87,44 @@ that is only used for development (tests, linting, docs, etc.). - Follow `docs/source/dev/README.rst` for style: maintain numpy-style docstrings, avoid type hints unless pre-existing, and keep module-level imports in the documented order (`numpy as np`, `xarray as xr`, etc.). -- Respect Ruff defaults (line length 79; docstrings 72). Run - `pixi run -e dev ruff check .` and `pixi run -e dev ruff format .` before - opening a PR. -- Favor descriptive names instead of extra comments; only add comments for - non-obvious behavior (e.g., tricky array ops or concurrency concerns). +- Follow Ruff configuration (79 char lines, 72 char lines for docstrings, + double quotes, numpy docstyle). Run locally: `pixi run -e dev ruff check .` + and `pixi run -e dev ruff format .` +- Do not add comments unless they are absolutely critical to understanding; + always prefer descriptive function and variable names instead. If some functionality needs further explanation, add this to the class/function/method docstring under the "Notes" section. - Use absolute imports under `revrt.`. - Surface warnings/errors through `revrt.warn` and `revrt.exceptions` (e.g., raise `revrtValueError` and emit `revrtWarning`) to ensure logging hooks fire. +- Logging: Use `logging.getLogger(__name__)`. Log at appropriate levels + (DEBUG, INFO, WARNING, ERROR, CRITICAL). Avoid print statements. +- Version: Never edit `_version.py` manually (auto-generated by `setuptools_scm`). + Tag releases `vX.Y.Z`. - When touching PyO3 bindings, update both Python shims in `revrt/_rust.py` (if needed) and the Rust exports to keep signatures aligned. ## 6. Docstring Guidelines (Python) - Use numpy-style docstrings; first line must omit a trailing period. -- Keep docstring lines ≤72 characters and avoid short summaries on - `__init__` methods. -- Document parameters in the method/function docstring, not the class-level - docstring; protected helpers (`_name`) should have single-sentence docs. +- Avoid type hints. +- Keep docstring length to 72 characters per line. +- Never include a period (".") at the end of the first line of docstrings. +- Do not add a short summary to __init__ methods. Instead, keep the line blank + and start the "Parameters" section after a second newline. +- Do not document parameters in the class docstring - do that in the __init__ + docstring instead. +- Do not add docstring to dunder methods (e.g., __str__, __repr__, etc.) + unless absolutely necessary. +- All @property and @cached_property method documentation should be one line + long and should start with the return type followed by a colon + (e.g. `"""str: My string property"""`). +- If the default value for a parameter is **not** `None`, document it using + the format: `param_name : type, default=`. If the default + value for a parameter **is** `None`, use the format : `param_name : type, optional`. +- "Protected" functions and methods (i.e. starting with an underscore) + should always be documented using **only** one-line summary docstrings. +- To exclude functions or classes from the public API documentation, start + the docstring with the token ``[NOT PUBLIC API]``. - Maintain intersphinx references where possible (see dev guide for mappings). ## 7. Coding Guidelines (Rust) @@ -128,16 +147,24 @@ that is only used for development (tests, linting, docs, etc.). `tests/conftest.py` for shared setup. Keep all test data small (individual files <1MB). The data doesn't need to fully reproduce realistic analysis cases - it just needs to include characteristics of a realistic case. +- All python test files (e.g. ``test_scenario.py``) should end with the + following block of code: + + .. code-block:: python + + if __name__ == "__main__": + pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) + + This allows the (single) file to be executed, running only the tests contained + within, which is extremely useful when updating/modifying/adding tests in the file. - Pytest options for parallel execution (`-n auto`) are supported; prefer `pixi run -e dev pytest -n auto` for heavier suites. ## 9. Logging, Errors, and Warnings - Do not log-and-raise manually; custom exceptions/warnings already emit log records. -- Prefer `revrt.utilities.log_mem` for memory-sensitive workflows to keep log - output consistent. -- CLI commands should rely on the logging configuration provided in - `revrt._cli.configure_logging` to avoid duplicate handlers. +- CLI commands should surface actionable messages and exit codes without + relying on hidden background logging. ## 10. Common Pitfalls & Gotchas - The PyO3 module must be rebuilt if Rust code changes; run a Pixi task (tests From d4968ac8d57bd11373f7d485cfb11336f6ff92d3 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 12:32:48 -0700 Subject: [PATCH 49/61] Break out dictionary extraction --- revrt/routing/cli.py | 65 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 62 insertions(+), 3 deletions(-) diff --git a/revrt/routing/cli.py b/revrt/routing/cli.py index 4bf9e840..11f6a43d 100644 --- a/revrt/routing/cli.py +++ b/revrt/routing/cli.py @@ -15,6 +15,7 @@ from revrt.costs.config import parse_config from revrt.routing.point_to_many import find_all_routes, RoutingScenario from revrt.routing.utilities import map_to_costs +from revrt.exceptions import revrtKeyError logger = logging.getLogger(__name__) @@ -406,14 +407,15 @@ def _update_multipliers(layers, polarity, voltage, transmission_config): for layer in output_layers: if layer.pop("apply_row_mult", False): - row_multiplier = transmission_config["row_width"][voltage] + row_multiplier = _get_row_multiplier(transmission_config, voltage) layer["multiplier_scalar"] = ( layer.get("multiplier_scalar", 1) * row_multiplier ) if layer.pop("apply_polarity_mult", False): - polarity_config = transmission_config["voltage_polarity_mult"] - polarity_multiplier = polarity_config[voltage][polarity] + polarity_multiplier = _get_polarity_multiplier( + transmission_config, voltage, polarity + ) layer["multiplier_scalar"] = ( layer.get("multiplier_scalar", 1) * polarity_multiplier @@ -423,6 +425,63 @@ def _update_multipliers(layers, polarity, voltage, transmission_config): return output_layers +def _get_row_multiplier(transmission_config, voltage): + try: + row_widths = transmission_config["row_width"] + except KeyError as e: + msg = ( + "`apply_row_mult` was set to `True`, but 'row_width' " + "not found in transmission config!" + ) + raise revrtKeyError(msg) from e + + try: + row_multiplier = row_widths["voltage"] + except KeyError as e: + msg = ( + "`apply_row_mult` was set to `True`, but voltage ' " + f"{voltage}' not found in transmission config " + "'row_width' settings. Available voltages: " + f"{list(row_multiplier)}" + ) + raise revrtKeyError(msg) from e + + return row_multiplier + + +def _get_polarity_multiplier(transmission_config, voltage, polarity): + try: + polarity_config = transmission_config["voltage_polarity_mult"] + except KeyError as e: + msg = ( + "`apply_polarity_mult` was set to `True`, but " + "'voltage_polarity_mult' not found in transmission config!" + ) + raise revrtKeyError(msg) from e + + try: + polarity_voltages = polarity_config[voltage] + except KeyError as e: + msg = ( + "`apply_polarity_mult` was set to `True`, but voltage ' " + f"{voltage}' not found in polarity config. Available voltages: " + f"{list(polarity_config)}" + ) + raise revrtKeyError(msg) from e + + try: + polarity_multiplier = polarity_voltages[polarity] + except KeyError as e: + msg = ( + "`apply_polarity_mult` was set to `True`, but polarity ' " + f"{polarity}' not found in voltage config. Available polarities: " + f"{list(polarity_voltages)}" + ) + raise revrtKeyError(msg) from e + + return polarity_multiplier + + def _route_points_subset(route_table, sort_cols, split_params): """Get indices of points that are sorted by location""" From eb0fe31f011466167cbbeb9b37f5a7710f9874ac Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 12:32:53 -0700 Subject: [PATCH 50/61] Minor update --- revrt/exceptions.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/revrt/exceptions.py b/revrt/exceptions.py index 3f4e28be..97271824 100644 --- a/revrt/exceptions.py +++ b/revrt/exceptions.py @@ -35,18 +35,10 @@ class revrtFileNotFoundError(revrtError, FileNotFoundError): """revrt FileNotFoundError""" -class revrtInvalidStartCostError(revrtError, ValueError): - """revrt InvalidStartCostError""" - - class revrtKeyError(revrtError, KeyError): """revrt KeyError""" -class revrtLeastCostPathNotFoundError(revrtError, RuntimeError): - """revrt LeastCostPathNotFoundError""" - - class revrtNotImplementedError(revrtError, NotImplementedError): """revrt NotImplementedError""" From bc2a0fc3c7f68ee707b3222f7f7e4ac3f94682eb Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 12:34:39 -0700 Subject: [PATCH 51/61] Docstrings --- revrt/routing/cli.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/revrt/routing/cli.py b/revrt/routing/cli.py index 11f6a43d..bf76e240 100644 --- a/revrt/routing/cli.py +++ b/revrt/routing/cli.py @@ -426,6 +426,7 @@ def _update_multipliers(layers, polarity, voltage, transmission_config): def _get_row_multiplier(transmission_config, voltage): + """Get right-of-way width multiplier for a given voltage""" try: row_widths = transmission_config["row_width"] except KeyError as e: @@ -450,6 +451,7 @@ def _get_row_multiplier(transmission_config, voltage): def _get_polarity_multiplier(transmission_config, voltage, polarity): + """Get multiplier for a given voltage and polarity""" try: polarity_config = transmission_config["voltage_polarity_mult"] except KeyError as e: From bb49f8d6bb6c8bd10d8266ca2df8a821699556a5 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 12:38:59 -0700 Subject: [PATCH 52/61] Add a few tests --- tests/python/unit/routing/test_routing_cli.py | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) diff --git a/tests/python/unit/routing/test_routing_cli.py b/tests/python/unit/routing/test_routing_cli.py index 279fe758..514fbf89 100644 --- a/tests/python/unit/routing/test_routing_cli.py +++ b/tests/python/unit/routing/test_routing_cli.py @@ -12,6 +12,7 @@ from revrt.utilities import LayeredFile from revrt.routing.utilities import map_to_costs +from revrt.exceptions import revrtKeyError from revrt.routing.cli import ( compute_lcp_routes, _run_lcp, @@ -20,6 +21,8 @@ _route_points_subset, _paths_to_compute, _split_routes, + _get_row_multiplier, + _get_polarity_multiplier, _MILLION_USD_PER_MILE_TO_USD_PER_PIXEL, ) @@ -499,5 +502,49 @@ def test_update_multipliers_applies_row_and_polarity(): assert unchanged[0]["layer_name"] == "layer_3" +def test_get_row_multiplier_missing_config(): + """_get_row_multiplier should raise when configuration keys are absent""" + + with pytest.raises(revrtKeyError) as excinfo: + _get_row_multiplier({}, "138") + assert "'row_width'" in str(excinfo.value) + + +def test_get_row_multiplier_unknown_voltage(): + """_get_row_multiplier should surface available voltages on failure""" + + config = {"row_width": {"230": 1.2}} + with pytest.raises(revrtKeyError) as excinfo: + _get_row_multiplier(config, "138") + assert "Available voltages" in str(excinfo.value) + assert "230" in str(excinfo.value) + + +def test_get_polarity_multiplier_missing_config(): + """_get_polarity_multiplier should raise when multiplier section missing""" + + with pytest.raises(revrtKeyError) as excinfo: + _get_polarity_multiplier({}, "138", "ac") + assert "voltage_polarity_mult" in str(excinfo.value) + + +def test_get_polarity_multiplier_unknown_voltage(): + """_get_polarity_multiplier should guard against unknown voltages""" + + config = {"voltage_polarity_mult": {"230": {"ac": 1.0}}} + with pytest.raises(revrtKeyError) as excinfo: + _get_polarity_multiplier(config, "138", "ac") + assert "Available voltages" in str(excinfo.value) + + +def test_get_polarity_multiplier_unknown_polarity(): + """_get_polarity_multiplier should guard against unknown polarities""" + + config = {"voltage_polarity_mult": {"138": {"dc": 1.0}}} + with pytest.raises(revrtKeyError) as excinfo: + _get_polarity_multiplier(config, "138", "ac") + assert "Available polarities" in str(excinfo.value) + + if __name__ == "__main__": pytest.main(["-q", "--show-capture=all", Path(__file__), "-rapP"]) From dbaad90624936d9dd635665178be403acf839846 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 12:41:10 -0700 Subject: [PATCH 53/61] Update --- tests/python/unit/routing/test_routing_cli.py | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/tests/python/unit/routing/test_routing_cli.py b/tests/python/unit/routing/test_routing_cli.py index 514fbf89..eecb8429 100644 --- a/tests/python/unit/routing/test_routing_cli.py +++ b/tests/python/unit/routing/test_routing_cli.py @@ -505,45 +505,39 @@ def test_update_multipliers_applies_row_and_polarity(): def test_get_row_multiplier_missing_config(): """_get_row_multiplier should raise when configuration keys are absent""" - with pytest.raises(revrtKeyError) as excinfo: + with pytest.raises(revrtKeyError, match="'row_width'"): _get_row_multiplier({}, "138") - assert "'row_width'" in str(excinfo.value) def test_get_row_multiplier_unknown_voltage(): """_get_row_multiplier should surface available voltages on failure""" config = {"row_width": {"230": 1.2}} - with pytest.raises(revrtKeyError) as excinfo: + with pytest.raises(revrtKeyError, match=r"Available voltages.*230"): _get_row_multiplier(config, "138") - assert "Available voltages" in str(excinfo.value) - assert "230" in str(excinfo.value) def test_get_polarity_multiplier_missing_config(): """_get_polarity_multiplier should raise when multiplier section missing""" - with pytest.raises(revrtKeyError) as excinfo: + with pytest.raises(revrtKeyError, match="voltage_polarity_mult"): _get_polarity_multiplier({}, "138", "ac") - assert "voltage_polarity_mult" in str(excinfo.value) def test_get_polarity_multiplier_unknown_voltage(): """_get_polarity_multiplier should guard against unknown voltages""" config = {"voltage_polarity_mult": {"230": {"ac": 1.0}}} - with pytest.raises(revrtKeyError) as excinfo: + with pytest.raises(revrtKeyError, match="Available voltages"): _get_polarity_multiplier(config, "138", "ac") - assert "Available voltages" in str(excinfo.value) def test_get_polarity_multiplier_unknown_polarity(): """_get_polarity_multiplier should guard against unknown polarities""" config = {"voltage_polarity_mult": {"138": {"dc": 1.0}}} - with pytest.raises(revrtKeyError) as excinfo: + with pytest.raises(revrtKeyError, match="Available polarities"): _get_polarity_multiplier(config, "138", "ac") - assert "Available polarities" in str(excinfo.value) if __name__ == "__main__": From cd2d55980335565c5822f77b6b58016e6f75f5c1 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 12:42:46 -0700 Subject: [PATCH 54/61] Minor fix --- revrt/routing/cli.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/revrt/routing/cli.py b/revrt/routing/cli.py index bf76e240..da5fc16e 100644 --- a/revrt/routing/cli.py +++ b/revrt/routing/cli.py @@ -437,13 +437,13 @@ def _get_row_multiplier(transmission_config, voltage): raise revrtKeyError(msg) from e try: - row_multiplier = row_widths["voltage"] + row_multiplier = row_widths[voltage] except KeyError as e: msg = ( "`apply_row_mult` was set to `True`, but voltage ' " f"{voltage}' not found in transmission config " "'row_width' settings. Available voltages: " - f"{list(row_multiplier)}" + f"{list(row_widths)}" ) raise revrtKeyError(msg) from e From 21fe9b6e47cef72acac42476bb0ee8da5eedc9d3 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 12:42:54 -0700 Subject: [PATCH 55/61] More messages --- tests/python/unit/routing/test_routing_cli.py | 41 ++++++++++++++++--- 1 file changed, 36 insertions(+), 5 deletions(-) diff --git a/tests/python/unit/routing/test_routing_cli.py b/tests/python/unit/routing/test_routing_cli.py index eecb8429..a44ef3df 100644 --- a/tests/python/unit/routing/test_routing_cli.py +++ b/tests/python/unit/routing/test_routing_cli.py @@ -505,7 +505,13 @@ def test_update_multipliers_applies_row_and_polarity(): def test_get_row_multiplier_missing_config(): """_get_row_multiplier should raise when configuration keys are absent""" - with pytest.raises(revrtKeyError, match="'row_width'"): + with pytest.raises( + revrtKeyError, + match=( + r"`apply_row_mult` was set to `True`, but 'row_width' not found " + r"in transmission config" + ), + ): _get_row_multiplier({}, "138") @@ -513,14 +519,27 @@ def test_get_row_multiplier_unknown_voltage(): """_get_row_multiplier should surface available voltages on failure""" config = {"row_width": {"230": 1.2}} - with pytest.raises(revrtKeyError, match=r"Available voltages.*230"): + with pytest.raises( + revrtKeyError, + match=( + r"`apply_row_mult` was set to `True`, but voltage '\s*138' not " + r"found in transmission config 'row_width' settings. " + r"Available voltages: \['230'\]" + ), + ): _get_row_multiplier(config, "138") def test_get_polarity_multiplier_missing_config(): """_get_polarity_multiplier should raise when multiplier section missing""" - with pytest.raises(revrtKeyError, match="voltage_polarity_mult"): + with pytest.raises( + revrtKeyError, + match=( + r"`apply_polarity_mult` was set to `True`, but " + r"'voltage_polarity_mult' not found in transmission config" + ), + ): _get_polarity_multiplier({}, "138", "ac") @@ -528,7 +547,13 @@ def test_get_polarity_multiplier_unknown_voltage(): """_get_polarity_multiplier should guard against unknown voltages""" config = {"voltage_polarity_mult": {"230": {"ac": 1.0}}} - with pytest.raises(revrtKeyError, match="Available voltages"): + with pytest.raises( + revrtKeyError, + match=( + r"`apply_polarity_mult` was set to `True`, but voltage '\s*138' " + r"not found in polarity config. Available voltages: \['230'\]" + ), + ): _get_polarity_multiplier(config, "138", "ac") @@ -536,7 +561,13 @@ def test_get_polarity_multiplier_unknown_polarity(): """_get_polarity_multiplier should guard against unknown polarities""" config = {"voltage_polarity_mult": {"138": {"dc": 1.0}}} - with pytest.raises(revrtKeyError, match="Available polarities"): + with pytest.raises( + revrtKeyError, + match=( + r"`apply_polarity_mult` was set to `True`, but polarity '\s*ac' " + r"not found in voltage config. Available polarities: \['dc'\]" + ), + ): _get_polarity_multiplier(config, "138", "ac") From 891b822789545b9c26c57b0650b3a82858ae9201 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 12:56:18 -0700 Subject: [PATCH 56/61] Use `da.max` --- revrt/routing/point_to_many.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/revrt/routing/point_to_many.py b/revrt/routing/point_to_many.py index 7f4ec666..348de390 100644 --- a/revrt/routing/point_to_many.py +++ b/revrt/routing/point_to_many.py @@ -211,7 +211,7 @@ def _build_final_routing_layer(self): self.final_routing_layer += friction_costs max_val = ( - np.max(self.final_routing_layer) * self.SOFT_BARRIER_MULTIPLIER + da.max(self.final_routing_layer) * self.SOFT_BARRIER_MULTIPLIER ) self.final_routing_layer = da.where( self.final_routing_layer <= 0, From 90dc0c307875d60b0173cadb8485c266d7db3c89 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Wed, 10 Dec 2025 13:23:00 -0700 Subject: [PATCH 57/61] Fix rust tests --- crates/revrt/src/dataset.rs | 58 ++++++++++++++++++++++++++++++++----- crates/revrt/src/lib.rs | 1 - 2 files changed, 50 insertions(+), 9 deletions(-) diff --git a/crates/revrt/src/dataset.rs b/crates/revrt/src/dataset.rs index e55e6c2c..7c2ddb8b 100644 --- a/crates/revrt/src/dataset.rs +++ b/crates/revrt/src/dataset.rs @@ -342,6 +342,13 @@ mod tests { for point in test_points { let results = dataset.get_3x3(&point); + // index 0, 0 has a cost of 0 and should therefore be filtered out + assert!( + !results + .iter() + .any(|(ArrayIndex { i, j }, _)| *i == 0 && *j == 0) + ); + for (ArrayIndex { i, j }, val) in results { let subset = zarrs::array_subset::ArraySubset::new_with_ranges(&[i..(i + 1), j..(j + 1)]); @@ -368,6 +375,13 @@ mod tests { for point in test_points { let results = dataset.get_3x3(&point); + // index 0, 0 has a cost of 0 and should therefore be filtered out + assert!( + !results + .iter() + .any(|(ArrayIndex { i, j }, _)| *i == 0 && *j == 0) + ); + for (ArrayIndex { i, j }, val) in results { let subset = zarrs::array_subset::ArraySubset::new_with_ranges(&[i..(i + 1), j..(j + 1)]); @@ -407,13 +421,20 @@ mod tests { let results = dataset.get_3x3(&ArrayIndex { i: 0, j: 0 }); + // index 0, 0 has a cost of 0 and should therefore be filtered out + assert!( + !results + .iter() + .any(|(ArrayIndex { i, j }, _)| *i == 0 && *j == 0) + ); + assert_eq!(results, vec![]); } #[test_case((0, 0), vec![(0, 1, 0.5), (1, 0, 1.0), (1, 1, 1.5 * SQRT_2)] ; "top left corner")] - #[test_case((0, 1), vec![(0, 0, 0.5), (1, 0, 1.5 * SQRT_2), (1, 1, 2.)] ; "top right corner")] - #[test_case((1, 0), vec![(0, 0, 1.), (0, 1, 1.5 * SQRT_2), (1, 1, 2.5)] ; "bottom left corner")] - #[test_case((1, 1), vec![(0, 0, 1.5 * SQRT_2), (0, 1, 2.), (1, 0, 2.5)] ; "bottom right corner")] + #[test_case((0, 1), vec![(1, 0, 1.5 * SQRT_2), (1, 1, 2.)] ; "top right corner")] + #[test_case((1, 0), vec![(0, 1, 1.5 * SQRT_2), (1, 1, 2.5)] ; "bottom left corner")] + #[test_case((1, 1), vec![(0, 1, 2.), (1, 0, 2.5)] ; "bottom right corner")] fn test_get_3x3_two_by_two_array((si, sj): (u64, u64), expected_output: Vec<(u64, u64, f32)>) { let path = samples::cost_as_index_zarr((2, 2), (2, 2)); let cost_function = @@ -423,6 +444,13 @@ mod tests { let results = dataset.get_3x3(&ArrayIndex { i: si, j: sj }); + // index 0, 0 has a cost of 0 and should therefore be filtered out + assert!( + !results + .iter() + .any(|(ArrayIndex { i, j }, _)| *i == 0 && *j == 0) + ); + assert_eq!( results, expected_output @@ -433,10 +461,10 @@ mod tests { } #[test_case((0, 0), vec![(0, 1, 0.5), (1, 0, 1.5), (1, 1, 2.0 * SQRT_2)] ; "top left corner")] - #[test_case((0, 1), vec![(0, 0, 0.5), (0, 2, 1.5), (1, 0, 2.0 * SQRT_2), (1, 1, 2.5), (1, 2, 3. * SQRT_2)] ; "top middle")] + #[test_case((0, 1), vec![(0, 2, 1.5), (1, 0, 2.0 * SQRT_2), (1, 1, 2.5), (1, 2, 3. * SQRT_2)] ; "top middle")] #[test_case((0, 2), vec![(0, 1, 1.5), (1, 1, 3.0 * SQRT_2), (1, 2, 3.5)] ; "top right corner")] - #[test_case((1, 0), vec![(0, 0, 1.5), (0, 1, 2.0 * SQRT_2), (1, 1, 3.5), (2, 0, 4.5), (2, 1, 5.0 * SQRT_2)] ; "middle left")] - #[test_case((1, 1), vec![(0, 0, 2.0 * SQRT_2), (0, 1, 2.5), (0, 2, 3.0 * SQRT_2), (1, 0, 3.5), (1, 2, 4.5), (2, 0, 5.0 * SQRT_2), (2, 1, 5.5), (2, 2, 6.0 * SQRT_2)] ; "middle middle")] + #[test_case((1, 0), vec![(0, 1, 2.0 * SQRT_2), (1, 1, 3.5), (2, 0, 4.5), (2, 1, 5.0 * SQRT_2)] ; "middle left")] + #[test_case((1, 1), vec![(0, 1, 2.5), (0, 2, 3.0 * SQRT_2), (1, 0, 3.5), (1, 2, 4.5), (2, 0, 5.0 * SQRT_2), (2, 1, 5.5), (2, 2, 6.0 * SQRT_2)] ; "middle middle")] #[test_case((1, 2), vec![(0, 1, 3.0 * SQRT_2), (0, 2, 3.5), (1, 1, 4.5), (2, 1, 6.0 * SQRT_2), (2, 2, 6.5)] ; "middle right")] #[test_case((2, 0), vec![(1, 0, 4.5), (1, 1, 5.0 * SQRT_2), (2, 1, 6.5)] ; "bottom left corner")] #[test_case((2, 1), vec![(1, 0, 5.0 * SQRT_2), (1, 1, 5.5), (1, 2, 6.0 * SQRT_2), (2, 0, 6.5), (2, 2, 7.5)] ; "bottom middle")] @@ -453,6 +481,13 @@ mod tests { let results = dataset.get_3x3(&ArrayIndex { i: si, j: sj }); + // index 0, 0 has a cost of 0 and should therefore be filtered out + assert!( + !results + .iter() + .any(|(ArrayIndex { i, j }, _)| *i == 0 && *j == 0) + ); + assert_eq!( results, expected_output @@ -463,10 +498,10 @@ mod tests { } #[test_case((0, 0), vec![(0, 1, 0.5), (1, 0, 2.), (1, 1, 2.5 * SQRT_2)] ; "top left corner")] - #[test_case((0, 1), vec![(0, 0, 0.5), (0, 2, 1.5), (1, 0, 2.5 * SQRT_2), (1, 1, 3.), (1, 2, 3.5 * SQRT_2)] ; "top left edge")] + #[test_case((0, 1), vec![(0, 2, 1.5), (1, 0, 2.5 * SQRT_2), (1, 1, 3.), (1, 2, 3.5 * SQRT_2)] ; "top left edge")] #[test_case((0, 2), vec![(0, 1, 1.5), (0, 3, 2.5), (1, 1, 3.5 * SQRT_2), (1, 2, 4.), (1, 3, 4.5 * SQRT_2)] ; "top right edge")] #[test_case((0, 3), vec![(0, 2, 2.5), (1, 2, 4.5 * SQRT_2), (1, 3, 5.)] ; "top right corner")] - #[test_case((1, 0), vec![(0, 0, 2.), (0, 1, 2.5 * SQRT_2), (1, 1, 4.5), (2, 0, 6.), (2, 1, 6.5 * SQRT_2)] ; "left top edge")] + #[test_case((1, 0), vec![(0, 1, 2.5 * SQRT_2), (1, 1, 4.5), (2, 0, 6.), (2, 1, 6.5 * SQRT_2)] ; "left top edge")] #[test_case((1, 3), vec![(0, 2, 4.5 * SQRT_2), (0, 3, 5.), (1, 2, 6.5), (2, 2, 8.5 * SQRT_2), (2, 3, 9.)] ; "right top edge")] #[test_case((2, 0), vec![(1, 0, 6.), (1, 1, 6.5 * SQRT_2), (2, 1, 8.5), (3, 0, 10.), (3, 1, 10.5 * SQRT_2)] ; "left bottom edge")] #[test_case((2, 3), vec![(1, 2, 8.5 * SQRT_2), (1, 3, 9.), (2, 2, 10.5), (3, 2, 12.5 * SQRT_2), (3, 3, 13.)] ; "right bottom edge")] @@ -486,6 +521,13 @@ mod tests { let results = dataset.get_3x3(&ArrayIndex { i: si, j: sj }); + // index 0, 0 has a cost of 0 and should therefore be filtered out + assert!( + !results + .iter() + .any(|(ArrayIndex { i, j }, _)| *i == 0 && *j == 0) + ); + assert_eq!( results, expected_output diff --git a/crates/revrt/src/lib.rs b/crates/revrt/src/lib.rs index 8ff8b2e3..3251f7cf 100644 --- a/crates/revrt/src/lib.rs +++ b/crates/revrt/src/lib.rs @@ -185,7 +185,6 @@ mod tests { #[test_case((1, 1), vec![(1, 3), (3, 1)], 1.; "horizontal and vertical")] #[test_case((3, 3), vec![(3, 5), (1, 1), (3, 1)], 1.; "horizontal")] #[test_case((3, 3), vec![(5, 3), (5, 5), (1, 3)], 1.; "vertical")] - #[test_case((3, 3), vec![(3, 1), (3, 4)], 0.; "zero costs")] fn routing_one_point_to_many_same_cost_and_length( (si, sj): (u64, u64), endpoints: Vec<(u64, u64)>, From 91d72d74215c186b7ab0e6183448dfe76aaa8116 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Thu, 11 Dec 2025 16:23:02 -0700 Subject: [PATCH 58/61] Minor refactor --- crates/revrt/src/dataset.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/crates/revrt/src/dataset.rs b/crates/revrt/src/dataset.rs index 6c9dc483..59fa67a4 100644 --- a/crates/revrt/src/dataset.rs +++ b/crates/revrt/src/dataset.rs @@ -276,6 +276,7 @@ impl Dataset { // Calculate the average with center point (half grid + other half grid). // Also, apply the diagonal factor for the extra distance. + // Finally, add any invariant costs. let cost_to_neighbors = neighbors .iter() .zip(invariant_neighbors.iter()) @@ -291,9 +292,8 @@ impl Dataset { } else { v }; - ((ir, jr), scaled + inv_cost) + (ArrayIndex { i: *ir, j: *jr }, scaled + inv_cost) }) - .map(|((ir, jr), v)| (ArrayIndex { i: *ir, j: *jr }, v)) .collect::>(); trace!("Neighbors {:?}", cost_to_neighbors); From 0dd2217157495c25c6f2a97964d49d5a5d2fc84b Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Thu, 11 Dec 2025 16:24:35 -0700 Subject: [PATCH 59/61] Much lower memory limit --- crates/revrt/src/dataset.rs | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/crates/revrt/src/dataset.rs b/crates/revrt/src/dataset.rs index 59fa67a4..8a4ca740 100644 --- a/crates/revrt/src/dataset.rs +++ b/crates/revrt/src/dataset.rs @@ -396,8 +396,7 @@ mod tests { let path = samples::multi_variable_zarr(); let cost_function = CostFunction::from_json(r#"{"cost_layers": [{"layer_name": "A"}]}"#).unwrap(); - let dataset = - Dataset::open(path, cost_function, 250_000_000).expect("Error opening dataset"); + let dataset = Dataset::open(path, cost_function, 1_000).expect("Error opening dataset"); let test_points = [ArrayIndex { i: 3, j: 1 }, ArrayIndex { i: 2, j: 2 }]; let array = zarrs::array::Array::open(dataset.source.clone(), "/A").unwrap(); From 1188678fa4c755305cb86ccc8413ecc895be8c28 Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Thu, 11 Dec 2025 17:03:08 -0700 Subject: [PATCH 60/61] `add_layer_to_data` now returns `Result` object --- crates/revrt/src/dataset.rs | 23 +++++++++++------------ crates/revrt/src/error.rs | 6 ++++++ crates/revrt/src/ffi.rs | 2 ++ 3 files changed, 19 insertions(+), 12 deletions(-) diff --git a/crates/revrt/src/dataset.rs b/crates/revrt/src/dataset.rs index 8a4ca740..fa2cba23 100644 --- a/crates/revrt/src/dataset.rs +++ b/crates/revrt/src/dataset.rs @@ -98,8 +98,8 @@ impl Dataset { let chunk_grid = tmp.chunk_grid(); debug!("Chunk grid info: {:?}", &chunk_grid); - add_layer_to_data("cost_invariant", chunk_grid, &swap); - add_layer_to_data("cost", chunk_grid, &swap); + add_layer_to_data("cost_invariant", chunk_grid, &swap)?; + add_layer_to_data("cost", chunk_grid, &swap)?; let cost_chunk_idx = ndarray::Array2::from_elem( ( @@ -360,29 +360,28 @@ fn add_layer_to_data( layer_name: &str, chunk_shape: &ChunkGrid, swap: &ReadableWritableListableStorage, -) { +) -> Result<()> { trace!("Creating an empty {} array", layer_name); let dataset_path = format!("/{layer_name}"); - zarrs::array::ArrayBuilder::new_with_chunk_grid( - // cost_shape, + let builder = zarrs::array::ArrayBuilder::new_with_chunk_grid( chunk_shape.clone(), zarrs::array::DataType::Float32, zarrs::array::FillValue::from(zarrs::array::ZARR_NAN_F32), - ) - .build(swap.clone(), &dataset_path) - .unwrap() - .store_metadata() - .unwrap(); + ); + + let built = builder.build(swap.clone(), &dataset_path)?; + built.store_metadata()?; - let array = zarrs::array::Array::open(swap.clone(), &dataset_path).unwrap(); + let array = zarrs::array::Array::open(swap.clone(), &dataset_path)?; trace!("'{}' shape: {:?}", layer_name, array.shape().to_vec()); trace!("'{}' chunk shape: {:?}", layer_name, array.chunk_grid()); trace!( "Dataset contents after '{}' creation: {:?}", layer_name, - swap.list().unwrap() + swap.list()? ); + Ok(()) } #[cfg(test)] diff --git a/crates/revrt/src/error.rs b/crates/revrt/src/error.rs index ee696053..cc9281c7 100644 --- a/crates/revrt/src/error.rs +++ b/crates/revrt/src/error.rs @@ -5,6 +5,12 @@ pub enum Error { #[error(transparent)] IO(#[from] std::io::Error), + #[error(transparent)] + ZarrsArrayCreate(#[from] zarrs::array::ArrayCreateError), + + #[error(transparent)] + ZarrsStorage(#[from] zarrs::storage::StorageError), + #[allow(dead_code)] #[error("Undefined error")] // Used during development while it is not clear a category of error diff --git a/crates/revrt/src/ffi.rs b/crates/revrt/src/ffi.rs index d060593a..36e87894 100644 --- a/crates/revrt/src/ffi.rs +++ b/crates/revrt/src/ffi.rs @@ -12,6 +12,8 @@ impl From for PyErr { fn from(err: Error) -> PyErr { match err { Error::IO(msg) => PyIOError::new_err(msg), + Error::ZarrsArrayCreate(e) => PyIOError::new_err(e.to_string()), + Error::ZarrsStorage(e) => PyIOError::new_err(e.to_string()), Error::Undefined(msg) => revrtRustError::new_err(msg), } } From 22e2519ddd659064ec0db62ab637bba13227709d Mon Sep 17 00:00:00 2001 From: ppinchuk Date: Thu, 11 Dec 2025 17:30:36 -0700 Subject: [PATCH 61/61] More error types --- crates/revrt/src/dataset.rs | 39 ++++++++++++++++++++++--------------- crates/revrt/src/error.rs | 3 +++ crates/revrt/src/ffi.rs | 1 + 3 files changed, 27 insertions(+), 16 deletions(-) diff --git a/crates/revrt/src/dataset.rs b/crates/revrt/src/dataset.rs index fa2cba23..78c66ab2 100644 --- a/crates/revrt/src/dataset.rs +++ b/crates/revrt/src/dataset.rs @@ -14,7 +14,7 @@ use zarrs::storage::{ use crate::ArrayIndex; use crate::cost::CostFunction; -use crate::error::Result; +use crate::error::{Error, Result}; pub(crate) use lazy_subset::LazySubset; /// Manages the features datasets and calculated total cost @@ -65,15 +65,13 @@ impl Dataset { trace!("Creating a new group for the cost dataset"); zarrs::group::GroupBuilder::new() - .build(swap.clone(), "/") - .unwrap() - .store_metadata() - .unwrap(); + .build(swap.clone(), "/")? + .store_metadata()?; let entries = source .list() .expect("failed to list variables in source dataset"); - let first_entry = entries + let first_entry_opt = entries .into_iter() .map(|entry| entry.to_string()) .find(|entry| { @@ -82,19 +80,28 @@ impl Dataset { const EXCLUDES: [&str; 6] = ["latitude", "longitude", "band", "x", "y", "spatial_ref"]; !name.ends_with(".json") && !EXCLUDES.iter().any(|needle| name == *needle) - }) - .unwrap_or_else(|| { - panic!( - "no non-coordinate variables found in source dataset: {:?}", - source - .list() - .unwrap_or_else(|err| panic!("failed to list dataset entries: {err}")) - ) }); + let first_entry = match first_entry_opt { + Some(e) => e, + None => { + return Err(Error::IO(std::io::Error::other(format!( + "no non-coordinate variables found in source dataset: {:?}", + source.list().ok() + )))); + } + }; + // Skip coordinate axes when selecting a representative variable for cost storage. - let varname = first_entry.split('/').next().unwrap().to_string(); + let varname = match first_entry.split('/').next() { + Some(name) => name, + None => { + return Err(Error::IO(std::io::Error::other( + "Could not determine any variable names from source dataset", + ))); + } + }; debug!("Using '{}' to determine shape of cost data", varname); - let tmp = zarrs::array::Array::open(source.clone(), &format!("/{varname}")).unwrap(); + let tmp = zarrs::array::Array::open(source.clone(), &format!("/{varname}"))?; let chunk_grid = tmp.chunk_grid(); debug!("Chunk grid info: {:?}", &chunk_grid); diff --git a/crates/revrt/src/error.rs b/crates/revrt/src/error.rs index cc9281c7..5d455031 100644 --- a/crates/revrt/src/error.rs +++ b/crates/revrt/src/error.rs @@ -8,6 +8,9 @@ pub enum Error { #[error(transparent)] ZarrsArrayCreate(#[from] zarrs::array::ArrayCreateError), + #[error(transparent)] + ZarrsGroupCreate(#[from] zarrs::group::GroupCreateError), + #[error(transparent)] ZarrsStorage(#[from] zarrs::storage::StorageError), diff --git a/crates/revrt/src/ffi.rs b/crates/revrt/src/ffi.rs index 36e87894..f31a1c09 100644 --- a/crates/revrt/src/ffi.rs +++ b/crates/revrt/src/ffi.rs @@ -14,6 +14,7 @@ impl From for PyErr { Error::IO(msg) => PyIOError::new_err(msg), Error::ZarrsArrayCreate(e) => PyIOError::new_err(e.to_string()), Error::ZarrsStorage(e) => PyIOError::new_err(e.to_string()), + Error::ZarrsGroupCreate(e) => PyIOError::new_err(e.to_string()), Error::Undefined(msg) => revrtRustError::new_err(msg), } }