diff --git a/pyforestscan/calculate.py b/pyforestscan/calculate.py index d25347e..030c96f 100644 --- a/pyforestscan/calculate.py +++ b/pyforestscan/calculate.py @@ -192,7 +192,9 @@ def calculate_canopy_cover(pad: np.ndarray, Returns: np.ndarray: 2D array (X, Y) of canopy cover values in [0, 1], with NaN where - PAD is entirely missing for the integration range. + PAD is entirely missing for the integration range. If the requested + integration range is empty (e.g., min_height >= available max height), + returns a zeros array (no canopy above the threshold). Raises: ValueError: If parameters are invalid (e.g., non-positive voxel_height, k < 0, @@ -203,16 +205,18 @@ def calculate_canopy_cover(pad: np.ndarray, if k < 0: raise ValueError(f"k must be >= 0 (got {k})") - # Compute PAI integrated from min_height up to max_height/top + # Determine effective max height and handle empty integration range + effective_max_height = max_height if max_height is not None else pad.shape[2] * voxel_height + if min_height >= effective_max_height: + # No foliage above threshold: cover is zero everywhere + return np.zeros((pad.shape[0], pad.shape[1]), dtype=float) + + # Compute PAI integrated from min_height up to effective_max_height/top pai_above = calculate_pai(pad, voxel_height, min_height=min_height, max_height=max_height) # Identify columns that are entirely NaN within the integration range - if max_height is None: - max_height = pad.shape[2] * voxel_height - if min_height >= max_height: - raise ValueError("Minimum height index must be less than maximum height index.") start_idx = int(np.ceil(min_height / voxel_height)) - end_idx = int(np.floor(max_height / voxel_height)) + end_idx = int(np.floor(effective_max_height / voxel_height)) range_slice = pad[:, :, start_idx:end_idx] all_nan_mask = np.all(np.isnan(range_slice), axis=2) diff --git a/tests/test_calculate.py b/tests/test_calculate.py index 613c374..603c0ef 100644 --- a/tests/test_calculate.py +++ b/tests/test_calculate.py @@ -646,3 +646,19 @@ def test_calculate_canopy_cover_monotonic_with_height(): cov5 = calculate_canopy_cover(pad, voxel_height=1.0, min_height=5.0, k=0.5) assert np.all(cov0 >= cov2) assert np.all(cov2 >= cov5) + + +def test_calculate_canopy_cover_empty_range_above_top_returns_zero(): + # Z extent is 0..1 m (2 layers with dz=1); threshold at 2 m leaves no range + pad = np.random.rand(4, 4, 2) + cov = calculate_canopy_cover(pad, voxel_height=1.0, min_height=2.0, k=0.5) + assert cov.shape == (4, 4) + assert np.allclose(cov, 0.0) + + +def test_calculate_canopy_cover_min_ge_max_returns_zero(): + pad = np.random.rand(3, 3, 5) + # Explicitly set max_height below min_height so integration is empty + cov = calculate_canopy_cover(pad, voxel_height=1.0, min_height=3.0, max_height=2.0, k=0.5) + assert cov.shape == (3, 3) + assert np.allclose(cov, 0.0)