diff --git a/.dockerignore b/.dockerignore index 396c99f..9df1f12 100644 --- a/.dockerignore +++ b/.dockerignore @@ -1,4 +1,5 @@ .git +.github .Rproj.user .Rhistory .RData @@ -14,4 +15,6 @@ .vscode __pycache__ data/ -tests/testthat/_snaps/ \ No newline at end of file +tests +tests/testthat/_snaps/ +inst/extdata diff --git a/.github/workflows/container.yml b/.github/workflows/container.yml new file mode 100644 index 0000000..bb1fd6b --- /dev/null +++ b/.github/workflows/container.yml @@ -0,0 +1,37 @@ +name: Container + +on: + push: + branches: [main] + tags: ["v*"] + workflow_dispatch: + +permissions: + contents: read + packages: write + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: docker/setup-buildx-action@v3 + - uses: docker/login-action@v3 + with: + registry: ghcr.io + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + - uses: docker/metadata-action@v5 + id: metadata + with: + images: ghcr.io/3dtrees-earth/3dtrees_csp_standsegmentation + tags: | + type=raw,value=0.2.0,enable=${{ github.ref == 'refs/heads/main' }} + type=ref,event=tag + type=sha + - uses: docker/build-push-action@v6 + with: + context: . + push: true + tags: ${{ steps.metadata.outputs.tags }} + labels: ${{ steps.metadata.outputs.labels }} diff --git a/.github/workflows/r.yml b/.github/workflows/r.yml deleted file mode 100644 index d9957c6..0000000 --- a/.github/workflows/r.yml +++ /dev/null @@ -1,62 +0,0 @@ -# Workflow derived from https://github.com/r-lib/actions/tree/master/examples -# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help -on: - push: - branches: [main] - pull_request: - branches: [main] - -name: R-CMD-check - -jobs: - R-CMD-check: - runs-on: ${{ matrix.config.os }} - - name: ${{ matrix.config.os }} (${{ matrix.config.r }}) - - strategy: - fail-fast: false - matrix: - config: - - {os: macOS-latest, r: 'release'} - - {os: windows-latest, r: 'release'} - - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} - - {os: ubuntu-latest, r: 'release'} - - {os: ubuntu-latest, r: 'oldrel-1'} - - env: - GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - R_KEEP_PKG_SOURCE: yes - - steps: - - uses: actions/checkout@v4 - - - uses: r-lib/actions/setup-pandoc@v2 - - - name: Install macOS system dependencies - if: runner.os == 'macos' - run: brew install gdal proj - - - uses: r-lib/actions/setup-r@v2 - with: - r-version: ${{ matrix.config.r }} - http-user-agent: ${{ matrix.config.http-user-agent }} - use-public-rspm: true - - - uses: r-lib/actions/setup-r-dependencies@v2 - with: - extra-packages: rcmdcheck - - - uses: r-lib/actions/check-r-package@v2 - - - name: Show testthat output - if: always() - run: find check -name 'testthat.Rout*' -exec cat '{}' \; || true - shell: bash - - - name: Upload check results - if: failure() - uses: actions/upload-artifact@v4 - with: - name: ${{ runner.os }}-r${{ matrix.config.r }}-results - path: check diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..398f372 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,56 @@ +FROM rocker/geospatial:4.4.3 + +ENV DEBIAN_FRONTEND=noninteractive \ + RGL_USE_NULL=TRUE + +ARG RLAS_COMMIT=82cbba42f158d1dfc91efda3207923260a052564 + +RUN apt-get update \ + && apt-get install --yes --no-install-recommends \ + libgl1-mesa-dev \ + libglu1-mesa-dev \ + libx11-dev \ + libxt-dev \ + curl \ + patch \ + python3-pip \ + && rm -rf /var/lib/apt/lists/* + +RUN python3 -m pip install --break-system-packages --no-cache-dir --no-deps \ + laspy==2.6.1 \ + lazrs==0.8.1 + +RUN install2.r --error --skipinstalled --ncpus -1 \ + BH \ + RANN \ + RCSF \ + Rcpp \ + RcppArmadillo \ + colorspace \ + conicfit \ + data.table \ + dbscan \ + doParallel \ + foreach \ + geometry \ + igraph \ + jsonlite \ + lidR \ + magrittr \ + rgl \ + testthat + +COPY patches/rlas-read-all-extrabytes.patch /tmp/rlas-read-all-extrabytes.patch +RUN curl -L -sS "https://github.com/r-lidar/rlas/archive/${RLAS_COMMIT}.tar.gz" -o /tmp/rlas.tar.gz \ + && mkdir -p /tmp/rlas-src \ + && tar -xzf /tmp/rlas.tar.gz -C /tmp/rlas-src --strip-components=1 \ + && patch -d /tmp/rlas-src -p1 < /tmp/rlas-read-all-extrabytes.patch \ + && R CMD INSTALL /tmp/rlas-src \ + && rm -rf /tmp/rlas-src /tmp/rlas.tar.gz /tmp/rlas-read-all-extrabytes.patch + +WORKDIR /opt/CspStandSegmentation +COPY . /opt/CspStandSegmentation +RUN chmod -R a+rX /opt/CspStandSegmentation \ + && R CMD INSTALL --no-multiarch --with-keep.source /opt/CspStandSegmentation + +CMD ["Rscript", "/opt/CspStandSegmentation/exec/run.R", "--help"] diff --git a/R/forest_inventory.R b/R/forest_inventory.R index 36eb979..a585910 100644 --- a/R/forest_inventory.R +++ b/R/forest_inventory.R @@ -37,6 +37,62 @@ suppress_cat <- function(f, ...) { f(...) # Call the function and capture its return value } +# Quiet, vectorized equivalent of conicfit::CircleFitByPratt(). The upstream +# implementation emits diagnostics with cat(), which previously required +# opening and sinking /dev/null for every RANSAC iteration. Forest inventory +# calls this thousands of times, making output suppression more expensive than +# the fit itself. +.circle_fit_by_pratt_quiet <- function(XY) { + XY <- as.matrix(XY) + centroid <- colMeans(XY) + centered <- sweep(XY, 2L, centroid, "-") + Xi <- centered[, 1L] + Yi <- centered[, 2L] + Zi <- Xi * Xi + Yi * Yi + + Mxx <- mean(Xi * Xi) + Myy <- mean(Yi * Yi) + Mxy <- mean(Xi * Yi) + Mxz <- mean(Xi * Zi) + Myz <- mean(Yi * Zi) + Mzz <- mean(Zi * Zi) + Mz <- Mxx + Myy + Cov_xy <- Mxx * Myy - Mxy * Mxy + Mxz2 <- Mxz * Mxz + Myz2 <- Myz * Myz + A2 <- 4 * Cov_xy - 3 * Mz * Mz - Mzz + A1 <- Mzz * Mz + 4 * Cov_xy * Mz - Mxz2 - Myz2 - Mz * Mz * Mz + A0 <- Mxz2 * Myy + Myz2 * Mxx - Mzz * Cov_xy - 2 * Mxz * Myz * Mxy + Mz * Mz * Cov_xy + A22 <- A2 + A2 + epsilon <- 1e-12 + ynew <- 1e20 + xnew <- 0 + + for (iter in seq_len(20L)) { + yold <- ynew + ynew <- A0 + xnew * (A1 + xnew * (A2 + xnew * xnew * 4)) + if (abs(ynew) > abs(yold)) { + xnew <- 0 + break + } + Dy <- A1 + xnew * (A22 + 16 * xnew * xnew) + xold <- xnew + xnew <- xold - ynew / Dy + if (abs((xnew - xold) / xnew) < epsilon) break + if (iter >= 20L || xnew < 0) { + xnew <- 0 + break + } + } + + DET <- xnew * xnew - xnew * Mz + Cov_xy + center <- c( + Mxz * (Myy - xnew) - Myz * Mxy, + Myz * (Mxx - xnew) - Mxz * Mxy + ) / DET / 2 + matrix(c(center + centroid, sqrt(sum(center * center) + Mz + 2 * xnew)), nrow = 1L) +} + #' RANSAC circle fitting algorithm specially adapted for tree DBH estimation #' #' This function fits a circle to a set of points using the RANSAC algorithm it maximizes the points that are in the circle and the number of filled 36 degree angle segments @@ -87,9 +143,7 @@ ransac_circle_fit <- function(data,n_iterations = 1000L,distance_threshold = 0.0 # Fit circle; keep tryCatch very tight and avoid pipe circle <- tryCatch( - { - CspStandSegmentation::suppress_cat(conicfit::CircleFitByPratt, sample_points) - }, + .circle_fit_by_pratt_quiet(sample_points), warning = function(w) NULL, error = function(e) NULL ) @@ -536,4 +590,3 @@ plot_inventory <- function(plot, inventory, col = NA, cex = 1.5, label_col = "wh rgl::lines3d(c(inventory$X[i] - plot[1], inventory$X[i] - plot[1]), c(inventory$Y[i] - plot[2], inventory$Y[i] - plot[2]), c(inventory$Z[i], inventory$Height[i]), col = ifelse(length(col) >= i, col[i], col), lwd = 2) } } - diff --git a/README.md b/README.md index d00b4db..f17b3f3 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,114 @@ Authors: Julian Frey and Zoe Schindler, University of Freiburg, Chair of Forest Growth and Dendroecology -[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.17294732.svg)](https://doi.org/10.5281/zenodo.17294732) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) [![R-CMD-check](https://github.com/JulFrey/CspStandSegmentation/actions/workflows/r.yml/badge.svg)](https://github.com/JulFrey/CspStandSegmentation/actions/workflows/r.yml) +[![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.17294732.svg)](https://doi.org/10.5281/zenodo.17294732) [![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0) + +## 3Dtrees command-line workflow + +This fork retains the upstream segmentation implementation and adds a headless +CLI and container for the 3Dtrees Galaxy tool. It always creates a DTM and can +inventory any number of existing instance dimensions. CSP segmentation is +optional; when enabled, the original points and dimensions are preserved and +the output gains exactly one `PredInstance_CSP` dimension. + +Inventory-only runs avoid loading dimensions that cannot affect the DTM or the +requested inventories. They retain XYZ, Classification, requested +instance/species fields, and optional ForestMamba score/semantic fields. CSP +runs continue to load every dimension because the emitted point cloud must +preserve them. LASlib can select only the first nine extra-byte records by +position; when a requested field occurs later, the reader safely falls back to +all extra bytes, then immediately projects the in-memory cloud back to the +required fields. The container uses the same pinned `rlas` patch as 3Dtrees +standardization, so that fallback correctly loads every declared extra byte +rather than stopping after nine. + +```bash +Rscript exec/run.R \ + --input input.laz \ + --output-dir results \ + --segmentation-spec PredInstance_SAT,species_id_SAT,species_prob_SAT \ + --segmentation-spec PredInstance_FM,species_id_FM,species_prob_FM \ + --enable-csp false \ + --dtm-resolution 0.2 \ + --random-seed 42 +``` + +`--segmentation-spec INSTANCE[,SPECIES,SPECIES_PROB]` is repeatable. Species +dimensions are optional per segmentation but must be supplied as a pair. When +none are supplied, the combined inventory omits species columns and no species +composition file is created. Run `Rscript exec/run.R --help` for all controls. + +Common controls are the input, repeatable segmentation specs, non-tree IDs, +optional CSP and seed source, optional native-CRS AOI GeoJSON, DTM resolution +(default 0.2 m), and random seed. Fine-tuning controls retain upstream defaults, +including a 0.3 m CSP voxel and one routing worker. CSP geometry features are +computed only when a non-zero geometry weight requires them. + +Outputs include: + +- `dtm_full.tif` and optional `dtm_aoi.tif`; +- one TSV per instance dimension plus `inventory_combined.tsv`; +- `stand_summary.tsv` and optional `species_composition.tsv`; +- `effective_seeds.tsv` and `segmented_csp.laz` only when CSP is enabled; +- `run_metadata.json` and `resource_summary.json`. + +Inventory rows include point count, position, height, DBH, crown convex-hull +area, and an explicit measurement-quality field. ForestMamba inventories also +include median `PredScore_FM`, a mixed-score flag, and wood/leaf point counts and +shares when those source dimensions exist. Processing is fail-atomic: the final +output directory is published only after all requested products succeed. + +Existing-instance inventory keeps upstream's 500 RANSAC iterations but uses a +quiet, vectorized implementation of the same Pratt circle equations. It also +projects only the point attributes required by each requested segmentation and +skips the full preservation copy when CSP output is disabled. On the local +5.9-million-point GFZ benchmark these changes reduced tool time from 59.4 to +41.4 seconds and process peak RSS from 4.85 to 3.19 GiB. + +With selective reading and the standardization `rlas` patch, the same GFZ file +successfully loaded 14 extra-byte attributes and inventoried SAT and FM together +in 53.7 seconds at 2.14 GiB process peak RSS. Sampled CPU averaged 100.1%, +confirming that the one-thread default consumes approximately one core. The run +produced 64 tree rows plus DTM, stand, and species products without creating a +point-cloud output. + +Inventory-only execution uses two bounded passes. DTM generation defaults to +300 m tiles with a 5 m buffer and up to 10 workers. Below 50 million points the +tiles read the source directly. Larger inputs are scanned once in parallel by +point range, retaining the minimum Z at each deterministic 0.1 m cell centre; +CSF and TIN rasterization then run on that reduced surface in parallel tiles. +Each spatial worker uses one lidR thread, avoiding nested oversubscription. The +result raster retains the aligned input extent, with unsupported edge areas as +NoData. + +The second pass reads only selected inventory fields into disk-backed hash +partitions by instance ID. Each tree remains complete even when its points span +spatial chunks. CSP continues to use the full-cloud path to preserve upstream +global voxel routing and the optional point-cloud output. + +On dataset 2056 (1,141,911,324 points, 3.7 GB compact LAZ), the memory-focused +streaming candidate stage completed in 56.8 seconds with 10 workers, averaged +9.51 CPU cores, and retained 1,267,661 candidates. Parent peak RSS was 0.78 GiB +and the conservative sum of all worker peaks was 6.72 GiB under a 50 GB Docker +limit. The 300 m + 5 m CSF/TIN stage completed in 76.4 seconds; its conservative +aggregate worker peak was 2.98 GiB. An earlier exact-coordinate prototype +scaled from 543.4 seconds with one worker to 306.6 seconds with two, 206.6 +seconds with four, and 127.9 seconds with ten. Replacing per-worker XY grids +with deterministic cell centres produced the final 56.8-second result and cut +the conservative 10-worker peak sum from 16.6 GiB to 6.72 GiB. + +On the 5.9-million-point GFZ reference, reducing to deterministic 0.1 m cell +centres before CSF/TIN changed the DTM relative to direct full-cloud CSF/TIN by +about 9.1 cm RMSE (4.0 cm median absolute difference). Automatic mode therefore +keeps the direct spatial method for clouds below the 50-million-point threshold. + +Build and run the pinned container with: + +```bash +docker build -t 3dtrees-csp . +docker run --rm -v "$PWD:/work" -w /work 3dtrees-csp \ + Rscript /opt/CspStandSegmentation/exec/run.R --help +``` @@ -161,4 +268,3 @@ BibTex: file = {Full Text PDF:O\:\\Research\\Projects\\Confobi_IWW\\Literatur\\lit_database\\storage\\R7Q8BFU5\\Larysch et al. - 2025 - Quantifying and mapping the ready-to-use veneer volume of European beech trees based on terrestrial.pdf:application/pdf}, } ``` - diff --git a/exec/run.R b/exec/run.R new file mode 100644 index 0000000..f8ec5ba --- /dev/null +++ b/exec/run.R @@ -0,0 +1,32 @@ +#!/usr/bin/env Rscript + +script_path <- function() { + arguments <- commandArgs(trailingOnly = FALSE) + file_argument <- grep("^--file=", arguments, value = TRUE) + if (!length(file_argument)) return(getwd()) + dirname(normalizePath(sub("^--file=", "", file_argument[[1L]]), mustWork = TRUE)) +} + +repository_root <- normalizePath(file.path(script_path(), ".."), mustWork = TRUE) +source(file.path(repository_root, "workflow", "cli.R")) +source(file.path(repository_root, "workflow", "tool.R")) + +if ("--version" %in% commandArgs(trailingOnly = TRUE)) { + description <- read.dcf(file.path(repository_root, "DESCRIPTION")) + cat(sprintf("CspStandSegmentation %s\n", description[1L, "Version"])) + quit(status = 0L) +} + +required_packages <- c("CspStandSegmentation", "data.table", "jsonlite", "lidR", "sf", "terra") +missing_packages <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)] +if (length(missing_packages)) { + abort(sprintf("Missing required R packages: %s", paste(missing_packages, collapse = ", "))) +} + +tryCatch( + run_tool(parse_cli_args()), + error = function(error) { + message(sprintf("ERROR: %s", conditionMessage(error))) + quit(status = 1L) + } +) diff --git a/patches/rlas-read-all-extrabytes.patch b/patches/rlas-read-all-extrabytes.patch new file mode 100644 index 0000000..4be87ed --- /dev/null +++ b/patches/rlas-read-all-extrabytes.patch @@ -0,0 +1,64 @@ +--- a/src/rlasstreamer.cpp ++++ b/src/rlasstreamer.cpp +@@ -160,37 +160,42 @@ + if (unselect.find("-C") != std::string::npos) read_cha(false); + if (unselect.find("-W") != std::string::npos) read_W(false); + +- std::vector select_eb(9); ++ std::vector select_eb(header->number_attributes); + std::fill(select_eb.begin(), select_eb.end(), false); ++ auto set_select_eb = [&](size_t index, bool value) ++ { ++ if (index < select_eb.size()) ++ select_eb[index] = value; ++ }; + + if (select.find("0") != std::string::npos) + std::fill(select_eb.begin(), select_eb.end(), true); + else + { +- if (select.find("1") != std::string::npos) select_eb[0] = true; +- if (select.find("2") != std::string::npos) select_eb[1] = true; +- if (select.find("3") != std::string::npos) select_eb[2] = true; +- if (select.find("4") != std::string::npos) select_eb[3] = true; +- if (select.find("5") != std::string::npos) select_eb[4] = true; +- if (select.find("6") != std::string::npos) select_eb[5] = true; +- if (select.find("7") != std::string::npos) select_eb[6] = true; +- if (select.find("8") != std::string::npos) select_eb[7] = true; +- if (select.find("9") != std::string::npos) select_eb[8] = true; ++ if (select.find("1") != std::string::npos) set_select_eb(0, true); ++ if (select.find("2") != std::string::npos) set_select_eb(1, true); ++ if (select.find("3") != std::string::npos) set_select_eb(2, true); ++ if (select.find("4") != std::string::npos) set_select_eb(3, true); ++ if (select.find("5") != std::string::npos) set_select_eb(4, true); ++ if (select.find("6") != std::string::npos) set_select_eb(5, true); ++ if (select.find("7") != std::string::npos) set_select_eb(6, true); ++ if (select.find("8") != std::string::npos) set_select_eb(7, true); ++ if (select.find("9") != std::string::npos) set_select_eb(8, true); + } + + if (unselect.find("-0") != std::string::npos) + std::fill(select_eb.begin(), select_eb.end(), false); + else + { +- if (unselect.find("-1") != std::string::npos) select_eb[0] = false; +- if (unselect.find("-2") != std::string::npos) select_eb[1] = false; +- if (unselect.find("-3") != std::string::npos) select_eb[2] = false; +- if (unselect.find("-4") != std::string::npos) select_eb[3] = false; +- if (unselect.find("-5") != std::string::npos) select_eb[4] = false; +- if (unselect.find("-6") != std::string::npos) select_eb[5] = false; +- if (unselect.find("-7") != std::string::npos) select_eb[6] = false; +- if (unselect.find("-8") != std::string::npos) select_eb[7] = false; +- if (unselect.find("-9") != std::string::npos) select_eb[8] = false; ++ if (unselect.find("-1") != std::string::npos) set_select_eb(0, false); ++ if (unselect.find("-2") != std::string::npos) set_select_eb(1, false); ++ if (unselect.find("-3") != std::string::npos) set_select_eb(2, false); ++ if (unselect.find("-4") != std::string::npos) set_select_eb(3, false); ++ if (unselect.find("-5") != std::string::npos) set_select_eb(4, false); ++ if (unselect.find("-6") != std::string::npos) set_select_eb(5, false); ++ if (unselect.find("-7") != std::string::npos) set_select_eb(6, false); ++ if (unselect.find("-8") != std::string::npos) set_select_eb(7, false); ++ if (unselect.find("-9") != std::string::npos) set_select_eb(8, false); + } + + IntegerVector pos_eb; diff --git a/tests/testthat/test-ransac-circle-fit.R b/tests/testthat/test-ransac-circle-fit.R new file mode 100644 index 0000000..150f7c2 --- /dev/null +++ b/tests/testthat/test-ransac-circle-fit.R @@ -0,0 +1,24 @@ +test_that("RANSAC circle fitting avoids per-iteration output redirection", { + body_text <- paste(deparse(body(ransac_circle_fit)), collapse = "\n") + expect_false(grepl("suppress_cat", body_text, fixed = TRUE)) +}) + +test_that("RANSAC circle fitting preserves partial-arc results without loop overhead", { + theta <- seq(0, pi, length.out = 50L) + points <- cbind( + 10 + 0.3 * cos(theta) + sin(seq_along(theta)) * 0.001, + 20 + 0.3 * sin(theta) + cos(seq_along(theta)) * 0.001 + ) + + set.seed(42L) + result <- ransac_circle_fit( + points, + n_iterations = 100L, + distance_threshold = 0.01, + min_inliers = 3L + ) + + expect_equal(as.numeric(result$circle), c(10.00015, 20.00363, 0.2983673), tolerance = 1e-5) + expect_equal(result$inliers, 50L) + expect_equal(result$angle_segs, 20L) +}) diff --git a/tests/tool/create_enriched_fixture.py b/tests/tool/create_enriched_fixture.py new file mode 100644 index 0000000..14b6cb9 --- /dev/null +++ b/tests/tool/create_enriched_fixture.py @@ -0,0 +1,51 @@ +#!/usr/bin/env python3 +"""Create a deterministic multi-attribute LAZ fixture from mikro_segmented.laz.""" + +from pathlib import Path +import sys + +import laspy +import numpy as np + + +def add_dimension(las: laspy.LasData, name: str, type_: str, description: str) -> None: + las.add_extra_dim(laspy.ExtraBytesParams(name=name, type=type_, description=description)) + + +def main() -> None: + if len(sys.argv) != 3: + raise SystemExit("usage: create_enriched_fixture.py INPUT.laz OUTPUT.laz") + input_path, output_path = map(Path, sys.argv[1:]) + las = laspy.read(input_path) + source_ids = np.asarray(las.PredInstance, dtype=np.int32) + valid = source_ids > 0 + + add_dimension(las, "PredInstance_FM", "int32", "Synthetic FM instance") + add_dimension(las, "species_id_FM", "int32", "Synthetic species ID") + add_dimension(las, "species_prob_FM", "float32", "Synthetic species probability") + add_dimension(las, "PredScore_FM", "float32", "Synthetic FM score") + add_dimension(las, "PredSemantic_FM", "int8", "Synthetic FM semantic class") + + las.PredInstance_FM = source_ids + species = np.where(valid, source_ids % 3 + 1, -1).astype(np.int32) + probability = np.where(valid, 0.65 + (source_ids % 4) * 0.05, -1).astype(np.float32) + score = np.where(valid, 0.70 + (source_ids % 5) * 0.04, -1).astype(np.float32) + semantic = np.where(valid, np.where(np.arange(len(las.points)) % 3 == 0, 1, 2), 0).astype(np.int8) + + first_id = next((value for value in np.unique(source_ids) if value > 0), None) + if first_id is not None: + indices = np.flatnonzero(source_ids == first_id) + if len(indices) >= 3: + species[indices[-1]] = species[indices[0]] + 10 + probability[indices[-1]] = 0.95 + score[indices[-1]] = 0.15 + + las.species_id_FM = species + las.species_prob_FM = probability + las.PredScore_FM = score + las.PredSemantic_FM = semantic + las.write(output_path) + + +if __name__ == "__main__": + main() diff --git a/tests/tool/test_cli.R b/tests/tool/test_cli.R new file mode 100644 index 0000000..b56c14d --- /dev/null +++ b/tests/tool/test_cli.R @@ -0,0 +1,136 @@ +#!/usr/bin/env Rscript + +all_arguments <- commandArgs(trailingOnly = FALSE) +test_file <- sub("^--file=", "", grep("^--file=", all_arguments, value = TRUE)[[1L]]) +root <- normalizePath(file.path(dirname(test_file), "..", ".."), mustWork = TRUE) +source(file.path(root, "workflow", "cli.R")) +source(file.path(root, "workflow", "tool.R")) + +expect_error <- function(expression, pattern) { + message <- tryCatch({ force(expression); NULL }, error = conditionMessage) + stopifnot(!is.null(message), grepl(pattern, message)) +} + +spec <- parse_segmentation_spec("PredInstance_FM,species_id_FM,species_prob_FM") +stopifnot( + identical(spec$source, "FM"), + identical(spec$species, "species_id_FM"), + identical(spec$species_prob, "species_prob_FM") +) +expect_error(parse_segmentation_spec("PredInstance_FM,species_id_FM"), "supplied together") +expect_error(parse_segmentation_spec("bad-name"), "valid LAS dimension") + +config <- parse_cli_args(c( + "--input", "input.laz", + "--output-dir", "output", + "--segmentation-spec", "PredInstance_SAT", + "--segmentation-spec", "PredInstance_FM,species_id_FM,species_prob_FM", + "--dtm-resolution", "0.2", + "--routing-workers", "1" +)) +stopifnot( + length(config$segmentation_specs) == 2L, + identical(config$dtm_resolution, 0.2), + identical(config$routing_workers, 1L), + identical(config$read_chunk_size, 300), + identical(config$chunk_buffer, 5), + identical(config$inventory_partitions, 64L), + identical(config$dtm_workers, 10L), + identical(config$dtm_strategy, "auto"), + identical(config$dtm_candidate_resolution, 0.1), + identical(config$dtm_streaming_threshold, 50000000L), + !config$enable_csp +) + +worker_pids <- unlist(parallel_chunk_map( + as.list(seq_len(4L)), + workers = 2L, + fun = function(value) { + Sys.sleep(0.05) + Sys.getpid() + } +)) +stopifnot(length(unique(worker_pids)) == 2L) + +edge_chunks <- spatial_chunks(list(xmin = 0, xmax = 25, ymin = 0, ymax = 25), 25) +stopifnot( + length(edge_chunks) == 4L, + any(vapply(edge_chunks, function(chunk) chunk[["xmin"]] == 25, logical(1))), + any(vapply(edge_chunks, function(chunk) chunk[["ymin"]] == 25, logical(1))) +) + +dimensions <- inventory_read_dimensions(config$segmentation_specs) +stopifnot( + identical( + dimensions, + c( + "PredInstance_SAT", "PredInstance_FM", "species_id_FM", + "species_prob_FM", "PredScore_FM", "PredSemantic_FM" + ) + ), + identical( + build_read_selector( + dimensions, + c( + "unused", "PredInstance_FM", "species_id_FM", "species_prob_FM", + "PredScore_FM", "PredSemantic_FM", "PredInstance_SAT" + ) + ), + "c234567" + ), + identical(build_read_selector("late_dimension", c(rep("unused", 9L), "late_dimension")), "c0") +) + +expect_error( + parse_cli_args(c("--input", "input.laz", "--output-dir", "output")), + "segmentation-spec" +) +expect_error( + parse_cli_args(c( + "--input", "input.laz", "--output-dir", "output", + "--segmentation-spec", "PredInstance_SAT", "--dtm-strategy", "invalid" + )), + "dtm-strategy" +) +expect_error( + parse_cli_args(c( + "--input", "input.laz", "--output-dir", "output", + "--segmentation-spec", "PredInstance_SAT", "--dtm-resolution", "0.2", + "--dtm-candidate-resolution", "0.3" + )), + "candidate-resolution" +) +expect_error( + parse_cli_args(c( + "--input", "input.laz", "--output-dir", "output", + "--enable-csp", "true", "--seed-mode", "supplied" + )), + "seed-file" +) + +points <- data.table::data.table( + PredInstance_FM = c(1, 1, 1, 2), + species_id_FM = c(4, 4, 7, -1), + species_prob_FM = c(0.8, 0.6, 0.9, -1), + PredScore_FM = c(0.9, 0.9, 0.2, -1), + PredSemantic_FM = c(1, 2, 2, 0) +) +species <- species_by_instance(points, spec) +tree_one <- species[PredInstance_FM == 1] +stopifnot( + tree_one$species_id == 4, + tree_one$species_prob == 0.7, + abs(tree_one$species_conflict_fraction - 1 / 3) < 1e-12 +) +fm <- fm_by_instance(points, "PredInstance_FM") +fm_one <- fm[PredInstance_FM == 1] +stopifnot( + fm_one$pred_score_fm == 0.9, + isTRUE(fm_one$pred_score_mixed), + fm_one$wood_point_count == 1, + fm_one$leaf_point_count == 2, + abs(fm_one$wood_share - 1 / 3) < 1e-12, + abs(fm_one$leaf_share - 2 / 3) < 1e-12 +) + +cat("tool helper tests passed\n") diff --git a/tests/tool/test_rlas_extra_bytes.R b/tests/tool/test_rlas_extra_bytes.R new file mode 100644 index 0000000..21ecd50 --- /dev/null +++ b/tests/tool/test_rlas_extra_bytes.R @@ -0,0 +1,27 @@ +#!/usr/bin/env Rscript + +point_count <- 20L +las <- lidR::LAS(data.frame( + X = as.numeric(seq_len(point_count)), + Y = as.numeric(seq_len(point_count) %% 3L), + Z = as.numeric(seq_len(point_count)) / 10, + Classification = rep(2L, point_count) +)) + +expected <- paste0("extra_", seq_len(12L)) +for (index in seq_along(expected)) { + las <- lidR::add_lasattribute( + las, + rep(as.numeric(index), point_count), + expected[[index]], + expected[[index]] + ) +} + +path <- tempfile(fileext = ".las") +on.exit(unlink(path), add = TRUE) +lidR::writeLAS(las, path) +loaded <- lidR::readLAS(path, select = "c0") + +stopifnot(all(expected %in% names(loaded@data))) +cat(sprintf("patched rlas loaded all %d extra-byte dimensions\n", length(expected))) diff --git a/workflow/cli.R b/workflow/cli.R new file mode 100644 index 0000000..62f7458 --- /dev/null +++ b/workflow/cli.R @@ -0,0 +1,217 @@ +abort <- function(message) { + stop(message, call. = FALSE) +} + +parse_bool <- function(value, name) { + normalized <- tolower(value) + if (normalized %in% c("true", "1", "yes")) return(TRUE) + if (normalized %in% c("false", "0", "no")) return(FALSE) + abort(sprintf("%s must be true or false", name)) +} + +parse_number <- function(value, name, minimum = -Inf, maximum = Inf, integer = FALSE) { + parsed <- suppressWarnings(as.numeric(value)) + if (length(parsed) != 1L || is.na(parsed) || parsed < minimum || parsed > maximum) { + abort(sprintf("%s must be between %s and %s", name, minimum, maximum)) + } + if (integer && parsed != as.integer(parsed)) abort(sprintf("%s must be an integer", name)) + if (integer) as.integer(parsed) else parsed +} + +parse_non_tree_ids <- function(value) { + values <- strsplit(value, ",", fixed = TRUE)[[1L]] + parsed <- suppressWarnings(as.numeric(trimws(values))) + if (anyNA(parsed)) abort("--non-tree-ids must be a comma-separated list of numbers") + unique(parsed) +} + +validate_dimension <- function(value, label) { + if (!grepl("^[A-Za-z][A-Za-z0-9_]*$", value)) { + abort(sprintf("%s is not a valid LAS dimension name: %s", label, value)) + } + value +} + +parse_segmentation_spec <- function(value) { + parts <- strsplit(value, ",", fixed = TRUE)[[1L]] + parts <- trimws(parts) + if (length(parts) > 3L || !length(parts) || !nzchar(parts[[1L]])) { + abort("--segmentation-spec must be INSTANCE[,SPECIES,SPECIES_PROB]") + } + parts <- c(parts, rep("", 3L - length(parts))) + has_species <- nzchar(parts[[2L]]) + has_probability <- nzchar(parts[[3L]]) + if (xor(has_species, has_probability)) { + abort(sprintf("Species and species-probability dimensions must be supplied together for %s", parts[[1L]])) + } + instance <- validate_dimension(parts[[1L]], "Instance dimension") + source <- sub("^PredInstance_", "", instance) + if (identical(source, instance)) source <- instance + list( + instance = instance, + source = source, + species = if (has_species) validate_dimension(parts[[2L]], "Species dimension") else NULL, + species_prob = if (has_probability) validate_dimension(parts[[3L]], "Species probability dimension") else NULL + ) +} + +usage <- function() { + paste( + "Usage: run.R --input FILE_OR_DIRECTORY --output-dir DIR [options]", + "", + "Common options:", + " --segmentation-spec INSTANCE[,SPECIES,SPECIES_PROB] Repeat for each existing segmentation", + " --non-tree-ids IDS Comma-separated IDs to exclude (default: 0)", + " --enable-csp true|false Create PredInstance_CSP (default: false)", + " --seed-mode automatic|supplied CSP seed source (default: automatic)", + " --seed-file TSV X/Y/Z/TreeID table for supplied seeds", + " --aoi-json GEOJSON Optional AOI in the point-cloud CRS", + " --dtm-resolution METRES DTM cell size (default: 0.2)", + " --read-chunk-size METRES Spatial DTM/inventory tile size (default: 300)", + " --chunk-buffer METRES DTM chunk buffer (default: 5)", + " --dtm-workers INTEGER Parallel DTM workers (default: 10)", + " --dtm-strategy auto|spatial|streaming DTM read strategy (default: auto)", + " --dtm-candidate-resolution METRES Streaming low-surface grid (default: 0.1)", + " --dtm-streaming-threshold POINTS Auto threshold (default: 50000000)", + " --inventory-partitions COUNT Disk-backed instance partitions (default: 64)", + " --random-seed INTEGER RANSAC seed (default: 42)", + "", + "Fine-tuning options:", + " --routing-workers INTEGER CSP routing workers (default: 1)", + " --geometry-threads INTEGER Geometry threads (default: 1)", + " --voxel-resolution METRES CSP routing voxel (default: 0.3)", + " --geometry-k INTEGER Geometry neighbours (default: 10)", + " --verticality-weight NUMBER CSP V_w (default: 0)", + " --linearity-weight NUMBER CSP L_w (default: 0)", + " --sphericity-weight NUMBER CSP S_w (default: 0)", + " --seed-resolution METRES Raster seed cell size (default: 0.1)", + " --seed-z-min METRES Seed slice minimum (default: 0.5)", + " --seed-z-max METRES Seed slice maximum (default: 2)", + " --seed-quantile NUMBER Seed density quantile (default: 0.975)", + " --seed-eps METRES Seed clustering radius (default: 0.2)", + " --slice-min METRES Inventory slice minimum (default: 0.3)", + " --slice-max METRES Inventory slice maximum (default: 4)", + " --slice-increment METRES Inventory slice step (default: 0.2)", + " --slice-width METRES Inventory slice width (default: 0.1)", + " --max-dbh METRES Maximum accepted DBH (default: 1)", + sep = "\n" + ) +} + +parse_cli_args <- function(args = commandArgs(trailingOnly = TRUE)) { + config <- list( + input = NULL, + output_dir = NULL, + segmentation_specs = list(), + non_tree_ids = 0, + enable_csp = FALSE, + seed_mode = "automatic", + seed_file = NULL, + aoi_json = NULL, + dtm_resolution = 0.2, + read_chunk_size = 300, + chunk_buffer = 5, + dtm_workers = 10L, + dtm_strategy = "auto", + dtm_candidate_resolution = 0.1, + dtm_streaming_threshold = 50000000L, + inventory_partitions = 64L, + random_seed = 42L, + routing_workers = 1L, + geometry_threads = 1L, + voxel_resolution = 0.3, + geometry_k = 10L, + verticality_weight = 0, + linearity_weight = 0, + sphericity_weight = 0, + seed_resolution = 0.1, + seed_z_min = 0.5, + seed_z_max = 2, + seed_quantile = 0.975, + seed_eps = 0.2, + slice_min = 0.3, + slice_max = 4, + slice_increment = 0.2, + slice_width = 0.1, + max_dbh = 1 + ) + + if (any(args %in% c("--help", "-h"))) { + cat(usage(), "\n") + quit(status = 0L) + } + + repeated <- character() + i <- 1L + while (i <= length(args)) { + key <- args[[i]] + if (!startsWith(key, "--")) abort(sprintf("Unexpected argument: %s", key)) + if (i == length(args)) abort(sprintf("Missing value for %s", key)) + value <- args[[i + 1L]] + i <- i + 2L + + if (key == "--input") config$input <- value + else if (key == "--output-dir") config$output_dir <- value + else if (key == "--segmentation-spec") repeated <- c(repeated, value) + else if (key == "--non-tree-ids") config$non_tree_ids <- parse_non_tree_ids(value) + else if (key == "--enable-csp") config$enable_csp <- parse_bool(value, key) + else if (key == "--seed-mode") config$seed_mode <- value + else if (key == "--seed-file") config$seed_file <- value + else if (key == "--aoi-json") config$aoi_json <- value + else if (key == "--dtm-resolution") config$dtm_resolution <- parse_number(value, key, 0.001) + else if (key == "--read-chunk-size") config$read_chunk_size <- parse_number(value, key, 1) + else if (key == "--chunk-buffer") config$chunk_buffer <- parse_number(value, key, 0) + else if (key == "--dtm-workers") config$dtm_workers <- parse_number(value, key, 1, 256, TRUE) + else if (key == "--dtm-strategy") config$dtm_strategy <- value + else if (key == "--dtm-candidate-resolution") config$dtm_candidate_resolution <- parse_number(value, key, 0.001) + else if (key == "--dtm-streaming-threshold") config$dtm_streaming_threshold <- parse_number(value, key, 1, Inf, TRUE) + else if (key == "--inventory-partitions") config$inventory_partitions <- parse_number(value, key, 1, 4096, TRUE) + else if (key == "--random-seed") config$random_seed <- parse_number(value, key, 0, .Machine$integer.max, TRUE) + else if (key == "--routing-workers") config$routing_workers <- parse_number(value, key, 1, 256, TRUE) + else if (key == "--geometry-threads") config$geometry_threads <- parse_number(value, key, 1, 256, TRUE) + else if (key == "--voxel-resolution") config$voxel_resolution <- parse_number(value, key, 0.001) + else if (key == "--geometry-k") config$geometry_k <- parse_number(value, key, 1, 1000, TRUE) + else if (key == "--verticality-weight") config$verticality_weight <- parse_number(value, key, 0, 1) + else if (key == "--linearity-weight") config$linearity_weight <- parse_number(value, key, 0, 1) + else if (key == "--sphericity-weight") config$sphericity_weight <- parse_number(value, key, 0, 1) + else if (key == "--seed-resolution") config$seed_resolution <- parse_number(value, key, 0.001) + else if (key == "--seed-z-min") config$seed_z_min <- parse_number(value, key) + else if (key == "--seed-z-max") config$seed_z_max <- parse_number(value, key) + else if (key == "--seed-quantile") config$seed_quantile <- parse_number(value, key, 0, 1) + else if (key == "--seed-eps") config$seed_eps <- parse_number(value, key, 0.001) + else if (key == "--slice-min") config$slice_min <- parse_number(value, key, 0) + else if (key == "--slice-max") config$slice_max <- parse_number(value, key, 0) + else if (key == "--slice-increment") config$slice_increment <- parse_number(value, key, 0.001) + else if (key == "--slice-width") config$slice_width <- parse_number(value, key, 0.001) + else if (key == "--max-dbh") config$max_dbh <- parse_number(value, key, 0.001) + else abort(sprintf("Unknown option: %s", key)) + } + + if (is.null(config$input)) abort("--input is required") + if (is.null(config$output_dir)) abort("--output-dir is required") + if (!config$dtm_strategy %in% c("auto", "spatial", "streaming")) { + abort("--dtm-strategy must be auto, spatial, or streaming") + } + if (config$dtm_candidate_resolution > config$dtm_resolution) { + abort("--dtm-candidate-resolution may not exceed --dtm-resolution") + } + if (!config$seed_mode %in% c("automatic", "supplied")) abort("--seed-mode must be automatic or supplied") + if (config$enable_csp && config$seed_mode == "supplied" && is.null(config$seed_file)) { + abort("--seed-file is required when --seed-mode supplied") + } + if (config$seed_z_min >= config$seed_z_max) abort("--seed-z-min must be lower than --seed-z-max") + if (config$slice_min >= config$slice_max) abort("--slice-min must be lower than --slice-max") + if (config$slice_width > config$slice_increment / 2) { + abort("--slice-width may not exceed half of --slice-increment") + } + config$segmentation_specs <- lapply(repeated, parse_segmentation_spec) + dimensions <- vapply(config$segmentation_specs, `[[`, character(1), "instance") + if (anyDuplicated(dimensions)) abort("Each instance dimension may only be supplied once") + if (config$enable_csp && "PredInstance_CSP" %in% dimensions) { + abort("Do not supply PredInstance_CSP when CSP is enabled; the tool creates it") + } + if (!config$enable_csp && !length(config$segmentation_specs)) { + abort("Provide at least one --segmentation-spec or enable CSP") + } + config +} diff --git a/workflow/stream_dtm_candidates.py b/workflow/stream_dtm_candidates.py new file mode 100644 index 0000000..1b4ccd7 --- /dev/null +++ b/workflow/stream_dtm_candidates.py @@ -0,0 +1,294 @@ +#!/usr/bin/env python3 +"""Create a bounded, parallel low-surface candidate cloud for CSF/TIN.""" + +from __future__ import annotations + +import argparse +from concurrent.futures import ProcessPoolExecutor, as_completed +from copy import deepcopy +import json +import math +from pathlib import Path +import resource +import shutil +import time + +import laspy +import numpy as np + + +def projection_records(header: laspy.LasHeader) -> list: + """Return CRS VLRs without requiring the optional pyproj dependency.""" + records = list(header.vlrs) + if header.evlrs is not None: + records.extend(header.evlrs) + return [record for record in records if record.user_id == "LASF_Projection"] + + +def projection_signature(header: laspy.LasHeader) -> tuple: + return tuple( + (record.user_id, record.record_id, record.record_data_bytes()) + for record in projection_records(header) + ) + + +def worker( + input_path: str, + parts_dir: str, + worker_id: int, + start: int, + count: int, + chunk_points: int, + grid_min_x: float, + grid_min_y: float, + rows: int, + cols: int, + resolution: float, +) -> tuple[int, str, int, float]: + started = time.monotonic() + minimum_z = np.full((rows, cols), np.inf, dtype=np.float32) + processed = 0 + + with laspy.open(input_path, laz_backend=laspy.LazBackend.Lazrs) as reader: + reader.seek(start) + remaining = count + while remaining: + points = reader.read_points(min(chunk_points, remaining)) + if not len(points): + break + x = np.asarray(points.x) + y = np.asarray(points.y) + z = np.asarray(points.z, dtype=np.float32) + col = np.floor((x - grid_min_x) / resolution).astype(np.int64) + row = np.floor((y - grid_min_y) / resolution).astype(np.int64) + np.clip(col, 0, cols - 1, out=col) + np.clip(row, 0, rows - 1, out=row) + cells = row * cols + col + flat_z = minimum_z.ravel() + np.minimum.at(flat_z, cells, z) + processed += len(points) + remaining -= len(points) + + if processed != count: + raise RuntimeError(f"DTM worker {worker_id} read {processed} of {count} points") + prefix = Path(parts_dir) / f"worker_{worker_id:03d}" + z_path = f"{prefix}_z.npy" + np.save(z_path, minimum_z) + return ( + worker_id, + z_path, + resource.getrusage(resource.RUSAGE_SELF).ru_maxrss, + time.monotonic() - started, + ) + + +def write_candidates( + output: Path, + source_header: laspy.LasHeader, + minimum_z: np.ndarray, + grid_min_x: float, + grid_min_y: float, + resolution: float, + row_block: int = 256, +) -> int: + header = laspy.LasHeader(point_format=0, version="1.4") + header.scales = np.array((0.001, 0.001, 0.001)) + header.offsets = np.array((grid_min_x, grid_min_y, 0.0)) + for record in projection_records(source_header): + header.vlrs.append(deepcopy(record)) + written = 0 + with laspy.open( + output, + mode="w", + header=header, + do_compress=True, + laz_backend=laspy.LazBackend.Lazrs, + ) as writer: + for start in range(0, minimum_z.shape[0], row_block): + stop = min(start + row_block, minimum_z.shape[0]) + valid_rows, valid_cols = np.nonzero(np.isfinite(minimum_z[start:stop])) + if not len(valid_rows): + continue + rows = valid_rows + start + points = laspy.ScaleAwarePointRecord.zeros(len(rows), header=header) + points.x = grid_min_x + (valid_cols.astype(np.float64) + 0.5) * resolution + points.y = grid_min_y + (rows.astype(np.float64) + 0.5) * resolution + points.z = minimum_z[rows, valid_cols] + writer.write_points(points) + written += len(points) + return written + + +def process_file( + input_path: Path, + parts_dir: Path, + workers: int, + chunk_points: int, + grid_min_x: float, + grid_min_y: float, + rows: int, + cols: int, + resolution: float, +) -> tuple[np.ndarray, int, int]: + with laspy.open(input_path) as reader: + point_count = int(reader.header.point_count) + worker_count = min(workers, point_count) + boundaries = np.linspace(0, point_count, worker_count + 1, dtype=np.int64) + tasks = [ + ( + str(input_path), + str(parts_dir), + index, + int(boundaries[index]), + int(boundaries[index + 1] - boundaries[index]), + chunk_points, + grid_min_x, + grid_min_y, + rows, + cols, + resolution, + ) + for index in range(worker_count) + ] + results = [] + with ProcessPoolExecutor(max_workers=worker_count) as executor: + futures = [executor.submit(worker, *task) for task in tasks] + completed = 0 + for future in as_completed(futures): + result = future.result() + results.append(result) + worker_id = result[0] + completed += tasks[worker_id][4] + print( + f"dtm_worker={worker_id} file={input_path.name} " + f"completed={completed}/{point_count} elapsed_seconds={result[-1]:.3f}", + flush=True, + ) + + minimum_z = np.full((rows, cols), np.inf, dtype=np.float32) + peak_sum_kib = 0 + for _, z_path, peak_kib, _ in sorted(results): + candidate_z = np.load(z_path, mmap_mode="r") + np.minimum(minimum_z, candidate_z, out=minimum_z) + peak_sum_kib += peak_kib + return minimum_z, peak_sum_kib, point_count + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("--input", action="append", required=True, type=Path) + parser.add_argument("--output-dir", required=True, type=Path) + parser.add_argument("--candidate-resolution", type=float, default=0.1) + parser.add_argument("--output-resolution", type=float, default=0.2) + parser.add_argument("--workers", type=int, default=10) + parser.add_argument("--chunk-points", type=int, default=1_000_000) + args = parser.parse_args() + started = time.monotonic() + + if args.workers < 1: + parser.error("--workers must be at least 1") + if args.chunk_points < 1: + parser.error("--chunk-points must be at least 1") + if args.candidate_resolution <= 0 or args.output_resolution <= 0: + parser.error("DTM resolutions must be greater than zero") + if args.candidate_resolution > args.output_resolution: + parser.error("--candidate-resolution may not exceed --output-resolution") + for path in args.input: + if not path.is_file(): + parser.error(f"input does not exist or is not a file: {path}") + + if args.output_dir.exists() and any(args.output_dir.iterdir()): + raise SystemExit(f"Output directory is not empty: {args.output_dir}") + args.output_dir.mkdir(parents=True, exist_ok=True) + parts_dir = args.output_dir / "parts" + parts_dir.mkdir() + + headers = [] + for path in args.input: + with laspy.open(path) as reader: + if reader.header.point_count < 1: + parser.error(f"input contains no points: {path}") + headers.append(reader.header) + source_crs = projection_signature(headers[0]) + for path, header in zip(args.input[1:], headers[1:]): + if projection_signature(header) != source_crs: + parser.error(f"input CRS differs from the first input: {path}") + input_min_x = min(float(header.mins[0]) for header in headers) + input_min_y = min(float(header.mins[1]) for header in headers) + input_max_x = max(float(header.maxs[0]) for header in headers) + input_max_y = max(float(header.maxs[1]) for header in headers) + grid_min_x = math.floor(input_min_x / args.output_resolution) * args.output_resolution + grid_min_y = math.floor(input_min_y / args.output_resolution) * args.output_resolution + grid_max_x = math.ceil(input_max_x / args.output_resolution) * args.output_resolution + grid_max_y = math.ceil(input_max_y / args.output_resolution) * args.output_resolution + cols = max(1, round((grid_max_x - grid_min_x) / args.candidate_resolution)) + rows = max(1, round((grid_max_y - grid_min_y) / args.candidate_resolution)) + + global_z = np.full((rows, cols), np.inf, dtype=np.float32) + total_points = 0 + worker_peak_sum_kib = 0 + for file_index, input_path in enumerate(args.input): + file_parts = parts_dir / f"file_{file_index:03d}" + file_parts.mkdir() + file_z, peak_kib, point_count = process_file( + input_path, + file_parts, + args.workers, + args.chunk_points, + grid_min_x, + grid_min_y, + rows, + cols, + args.candidate_resolution, + ) + np.minimum(global_z, file_z, out=global_z) + total_points += point_count + worker_peak_sum_kib = max(worker_peak_sum_kib, peak_kib) + shutil.rmtree(file_parts) + + candidate_path = args.output_dir / "minimum_candidates.laz" + candidate_points = write_candidates( + candidate_path, + headers[0], + global_z, + grid_min_x, + grid_min_y, + args.candidate_resolution, + ) + elapsed_seconds = time.monotonic() - started + self_usage = resource.getrusage(resource.RUSAGE_SELF) + child_usage = resource.getrusage(resource.RUSAGE_CHILDREN) + cpu_user_seconds = self_usage.ru_utime + child_usage.ru_utime + cpu_system_seconds = self_usage.ru_stime + child_usage.ru_stime + metadata = { + "input_files": [str(path) for path in args.input], + "input_points": total_points, + "candidate_points": candidate_points, + "candidate_resolution_m": args.candidate_resolution, + "candidate_xy_mode": "cell_center", + "output_resolution_m": args.output_resolution, + "min_x": grid_min_x, + "min_y": grid_min_y, + "max_x": grid_max_x, + "max_y": grid_max_y, + "rows": rows, + "cols": cols, + "workers": min(args.workers, total_points), + "chunk_points": args.chunk_points, + "worker_peak_rss_sum_upper_bound_kb": worker_peak_sum_kib, + "parent_peak_rss_kb": self_usage.ru_maxrss, + "cpu_user_seconds": cpu_user_seconds, + "cpu_system_seconds": cpu_system_seconds, + "average_cpu_cores": (cpu_user_seconds + cpu_system_seconds) / elapsed_seconds, + "elapsed_seconds": elapsed_seconds, + } + (args.output_dir / "candidate_metadata.json").write_text( + json.dumps(metadata, indent=2) + "\n", + encoding="utf-8", + ) + print(json.dumps(metadata, sort_keys=True), flush=True) + shutil.rmtree(parts_dir) + + +if __name__ == "__main__": + main() diff --git a/workflow/tool.R b/workflow/tool.R new file mode 100644 index 0000000..47885fb --- /dev/null +++ b/workflow/tool.R @@ -0,0 +1,997 @@ +resolve_input_files <- function(path) { + if (file.exists(path) && !dir.exists(path)) return(normalizePath(path, mustWork = TRUE)) + if (!dir.exists(path)) abort(sprintf("Input does not exist: %s", path)) + files <- list.files(path, pattern = "\\.(las|laz)$", full.names = TRUE, ignore.case = TRUE) + if (!length(files)) abort(sprintf("Input directory contains no LAS/LAZ files: %s", path)) + sort(normalizePath(files, mustWork = TRUE)) +} + +inventory_read_dimensions <- function(specs) { + dimensions <- unlist(lapply(specs, function(spec) { + result <- c(spec$instance, spec$species, spec$species_prob) + if (identical(spec$source, "FM")) { + result <- c(result, "PredScore_FM", "PredSemantic_FM") + } + result + }), use.names = FALSE) + unique(dimensions[!is.na(dimensions) & nzchar(dimensions)]) +} + +extra_byte_names <- function(header) { + extra_bytes <- header@VLR$Extra_Bytes[["Extra Bytes Description"]] + if (is.null(extra_bytes)) character() else names(extra_bytes) +} + +build_read_selector <- function(dimensions, extra_dimensions) { + standard_codes <- c( + gpstime = "t", ScanAngle = "a", Intensity = "i", + NumberOfReturns = "n", ReturnNumber = "r", Classification = "c", + Synthetic_flag = "s", Keypoint_flag = "k", Withheld_flag = "w", + Overlap_flag = "o", UserData = "u", PointSourceID = "p", + EdgeOfFlightline = "e", ScanDirectionFlag = "d", R = "R", G = "G", + B = "B", NIR = "N", ScannerChannel = "C", Waveform = "W" + ) + standard <- unname(standard_codes[intersect(dimensions, names(standard_codes))]) + positions <- sort(match(intersect(dimensions, extra_dimensions), extra_dimensions)) + + # LASlib can address only the first nine extra-byte records individually. + # If a requested dimension is later, load all extra bytes rather than risk + # silently omitting a requested inventory field. + extra <- if (any(positions > 9L)) "0" else as.character(positions) + paste0(unique(c("c", standard, extra)), collapse = "") +} + +read_selector <- function(file, specs, preserve_all_dimensions) { + if (preserve_all_dimensions) return("*") + header <- lidR::readLASheader(file) + build_read_selector(inventory_read_dimensions(specs), extra_byte_names(header)) +} + +project_inventory_dimensions <- function(las, specs) { + available <- names(las@data) + keep <- unique(c( + "X", "Y", "Z", intersect("Classification", available), + intersect(inventory_read_dimensions(specs), available) + )) + dropped <- setdiff(available, keep) + las@data <- las@data[, ..keep] + + extra_bytes <- las@header@VLR$Extra_Bytes[["Extra Bytes Description"]] + if (!is.null(extra_bytes)) { + for (dimension in setdiff(names(extra_bytes), keep)) { + las@header@VLR$Extra_Bytes[["Extra Bytes Description"]][[dimension]] <- NULL + } + } + attr(las, "dropped_dimensions") <- dropped + las +} + +read_point_cloud <- function(files, specs, preserve_all_dimensions = FALSE) { + selectors <- vapply( + files, + read_selector, + character(1), + specs = specs, + preserve_all_dimensions = preserve_all_dimensions + ) + clouds <- Map(function(file, select) { + cloud <- suppressWarnings(lidR::readLAS(file, select = select)) + if (is.null(cloud) || lidR::is.empty(cloud)) return(cloud) + if (!preserve_all_dimensions) cloud <- project_inventory_dimensions(cloud, specs) + cloud + }, files, selectors) + if (any(vapply(clouds, function(cloud) is.null(cloud) || lidR::is.empty(cloud), logical(1)))) { + abort("At least one input point-cloud file is empty") + } + dropped_dimensions <- unique(unlist(lapply(clouds, attr, which = "dropped_dimensions"))) + if (is.null(dropped_dimensions)) dropped_dimensions <- character() + if (length(files) == 1L) { + las <- clouds[[1L]] + } else { + las <- CspStandSegmentation::las_merge(clouds, fill = TRUE) + } + if (is.null(las) || lidR::is.empty(las)) abort("The input point cloud is empty") + attr(las, "read_selectors") <- stats::setNames(unname(selectors), basename(files)) + attr(las, "dropped_dimensions") <- dropped_dimensions + las +} + +has_valid_ground <- function(las) { + if (!"Classification" %in% names(las@data)) return(FALSE) + ground <- las@data[Classification == 2 & is.finite(X) & is.finite(Y)] + nrow(ground) >= 3L && length(unique(ground$X)) >= 2L && length(unique(ground$Y)) >= 2L +} + +make_dtm <- function(las, resolution) { + if (has_valid_ground(las)) { + grounded <- las + method <- "classification_2" + } else { + grounded <- lidR::classify_ground(las, lidR::csf(), last_returns = FALSE) + if (!has_valid_ground(grounded)) abort("Ground classification did not produce enough ground points for a DTM") + method <- "csf" + } + dtm <- lidR::rasterize_terrain(grounded, res = resolution, algorithm = lidR::tin()) + if (is.null(dtm) || all(is.na(terra::values(dtm)))) abort("DTM generation produced no valid cells") + list(las = grounded, dtm = dtm, method = method) +} + +read_aoi <- function(path, las_crs) { + if (is.null(path)) return(NULL) + if (!file.exists(path)) abort(sprintf("AOI JSON does not exist: %s", path)) + aoi <- suppressWarnings(sf::st_read(path, quiet = TRUE)) + if (!nrow(aoi)) abort("AOI JSON contains no features") + types <- unique(as.character(sf::st_geometry_type(aoi))) + if (!all(types %in% c("POLYGON", "MULTIPOLYGON"))) { + abort("AOI JSON must contain only Polygon or MultiPolygon geometries") + } + # AOI coordinates are defined by the CLI contract to already be in the + # point-cloud CRS. GeoJSON readers otherwise label them as WGS84 even when + # they contain projected native coordinates, so assign rather than transform. + sf::st_crs(aoi) <- las_crs + aoi <- sf::st_make_valid(aoi) + union <- sf::st_union(sf::st_geometry(aoi)) + if (length(union) != 1L || sf::st_is_empty(union)) abort("AOI geometry is empty after normalization") + union +} + +polygon_area <- function(x, y) { + if (length(x) < 3L) return(NA_real_) + indices <- chull(x, y) + xx <- x[indices] + yy <- y[indices] + 0.5 * abs(sum(xx * c(yy[-1L], yy[1L]) - yy * c(xx[-1L], xx[1L]))) +} + +determine_area <- function(las, aoi) { + if (!is.null(aoi)) { + area <- as.numeric(sf::st_area(aoi)) + return(list(area_m2 = area, area_source = "aoi")) + } + list( + area_m2 = polygon_area(las@data$X, las@data$Y), + area_source = "point_cloud_convex_hull" + ) +} + +validate_specs <- function(las, specs) { + available <- names(las@data) + for (spec in specs) { + required <- c(spec$instance, spec$species, spec$species_prob) + missing <- setdiff(required[!vapply(required, is.null, logical(1))], available) + if (length(missing)) { + abort(sprintf("Requested dimensions are missing for %s: %s", spec$instance, paste(missing, collapse = ", "))) + } + } +} + +quality_label <- function(code, dbh) { + labels <- c( + `1` = "too_few_stem_slice_points", + `2` = "fallback_circle_estimate", + `3` = "rejected_or_unstable_fit", + `4` = "spline_fit" + ) + result <- unname(labels[as.character(code)]) + result[is.na(result)] <- "measurement_unavailable" + result[is.na(dbh)] <- "dbh_unavailable" + result +} + +modal_value <- function(values) { + values <- values[!is.na(values) & is.finite(values) & values >= 0] + if (!length(values)) return(NA_real_) + counts <- table(values) + as.numeric(sort(names(counts)[counts == max(counts)])[1L]) +} + +species_by_instance <- function(points, spec) { + id <- spec$instance + species <- spec$species + probability <- spec$species_prob + points[, { + chosen <- modal_value(get(species)) + valid_species <- get(species)[!is.na(get(species)) & is.finite(get(species)) & get(species) >= 0] + matching_prob <- get(probability)[get(species) == chosen & is.finite(get(probability)) & get(probability) >= 0] + list( + species_id = chosen, + species_prob = if (length(matching_prob)) stats::median(matching_prob) else NA_real_, + species_conflict_fraction = if (length(valid_species)) 1 - sum(valid_species == chosen) / length(valid_species) else NA_real_ + ) + }, by = id] +} + +fm_by_instance <- function(points, instance) { + has_score <- "PredScore_FM" %in% names(points) + has_semantic <- "PredSemantic_FM" %in% names(points) + points[, { + scores <- if (has_score) get("PredScore_FM") else numeric() + scores <- scores[is.finite(scores) & scores >= 0] + semantics <- if (has_semantic) get("PredSemantic_FM") else numeric() + wood <- if (has_semantic) sum(semantics == 1, na.rm = TRUE) else NA_integer_ + leaf <- if (has_semantic) sum(semantics == 2, na.rm = TRUE) else NA_integer_ + classified <- if (has_semantic) wood + leaf else 0L + list( + pred_score_fm = if (length(scores)) stats::median(scores) else NA_real_, + pred_score_mixed = if (length(scores)) data.table::uniqueN(scores) > 1L else NA, + wood_point_count = wood, + leaf_point_count = leaf, + wood_share = if (classified > 0L) wood / classified else NA_real_, + leaf_share = if (classified > 0L) leaf / classified else NA_real_ + ) + }, by = instance] +} + +filter_inventory_to_aoi <- function(inventory, aoi, crs) { + if (is.null(aoi) || !nrow(inventory)) return(inventory) + points <- sf::st_as_sf(inventory, coords = c("x", "y"), crs = crs, remove = FALSE) + inventory[lengths(sf::st_intersects(points, aoi)) > 0L, ] +} + +inventory_for_spec <- function(las, spec, config, aoi, crs) { + instance <- spec$instance + available <- names(las@data) + point_columns <- unique(c( + "X", "Y", "Z", intersect("Zref", available), + instance, spec$species, spec$species_prob, + if (identical(spec$source, "FM")) intersect(c("PredScore_FM", "PredSemantic_FM"), available) + )) + points <- las@data[ + !is.na(get(instance)) & is.finite(get(instance)) & + get(instance) >= 0 & !(get(instance) %in% config$non_tree_ids), + ..point_columns + ] + if (!nrow(points)) abort(sprintf("No valid instances remain for requested dimension %s", instance)) + + z_original <- if ("Zref" %in% names(points)) "Zref" else "Z" + base <- points[, list( + point_count = .N, + fallback_x = stats::median(X), + fallback_y = stats::median(Y), + fallback_z = min(get(z_original), na.rm = TRUE), + fallback_height = max(Z, na.rm = TRUE) - min(Z, na.rm = TRUE), + fallback_hull = polygon_area(X, Y) + ), by = instance] + data.table::setnames(base, instance, "instance_id") + + measured <- tryCatch( + CspStandSegmentation::forest_inventory( + las, + slice_min = config$slice_min, + slice_max = config$slice_max, + increment = config$slice_increment, + width = config$slice_width, + max_dbh = config$max_dbh, + n_cores = config$geometry_threads, + tree_id_col = instance, + non_tree_id = config$non_tree_ids + ), + error = function(error) { + warning(sprintf("Inventory measurements failed for %s: %s", instance, conditionMessage(error))) + NULL + } + ) + + # Upstream simplifies a one-tree inventory to an 8x1 matrix/data frame whose + # field names are row names. Restore the same one-row schema returned for + # multi-tree calls before applying the normal merge path. + if (!is.null(measured) && nrow(measured) && !instance %in% names(measured) && ncol(measured) == 1L) { + fields <- rownames(measured) + required_fields <- c(instance, "X", "Y", "Z", "DBH", "quality_flag", "Height", "ConvexHullArea") + if (all(required_fields %in% fields)) { + values <- as.numeric(measured[[1L]]) + names(values) <- fields + measured <- data.table::as.data.table(as.list(values)) + } + } + + if (is.null(measured) || !nrow(measured) || !instance %in% names(measured)) { + if (!is.null(measured) && nrow(measured) && !instance %in% names(measured)) { + warning(sprintf( + "Inventory measurements for %s omitted the instance column; using fallback measurements for this partition", + instance + )) + } + measured <- data.table::data.table( + instance_id = numeric(), X = numeric(), Y = numeric(), Z = numeric(), + DBH = numeric(), quality_flag = integer(), Height = numeric(), ConvexHullArea = numeric() + ) + } else { + measured <- data.table::as.data.table(measured) + data.table::setnames(measured, instance, "instance_id") + } + + result <- merge(base, measured, by = "instance_id", all.x = TRUE, sort = FALSE) + result[, `:=`( + segmentation_source = spec$source, + instance_dimension = instance, + x = data.table::fcoalesce(X, fallback_x), + y = data.table::fcoalesce(Y, fallback_y), + z = data.table::fcoalesce(Z, fallback_z), + dbh_m = DBH, + height_m = data.table::fcoalesce(Height, fallback_height), + convex_hull_area_m2 = data.table::fcoalesce(ConvexHullArea, fallback_hull), + inventory_quality_code = data.table::fcoalesce(as.integer(quality_flag), 9L) + )] + result[, inventory_quality_label := quality_label(inventory_quality_code, dbh_m)] + + if (!is.null(spec$species)) { + species <- species_by_instance(points, spec) + data.table::setnames(species, instance, "instance_id") + result <- merge(result, species, by = "instance_id", all.x = TRUE, sort = FALSE) + } + + if (identical(instance, "PredInstance_FM") || identical(spec$source, "FM")) { + fm <- fm_by_instance(points, instance) + data.table::setnames(fm, instance, "instance_id") + result <- merge(result, fm, by = "instance_id", all.x = TRUE, sort = FALSE) + } + + keep <- c( + "segmentation_source", "instance_dimension", "instance_id", "x", "y", "z", + "dbh_m", "height_m", "convex_hull_area_m2", "point_count", + "species_id", "species_prob", "species_conflict_fraction", + "pred_score_fm", "pred_score_mixed", "wood_point_count", "leaf_point_count", + "wood_share", "leaf_share", "inventory_quality_code", "inventory_quality_label" + ) + result <- result[, intersect(keep, names(result)), with = FALSE] + result <- filter_inventory_to_aoi(result, aoi, crs) + attr(result, "segmentation_source") <- spec$source + attr(result, "instance_dimension") <- instance + result +} + +safe_stat <- function(values, function_) { + values <- values[is.finite(values)] + if (length(values)) function_(values) else NA_real_ +} + +safe_quantile <- function(values, probability) { + safe_stat(values, function(x) as.numeric(stats::quantile(x, probability, names = FALSE))) +} + +stand_summary_for_inventory <- function(inventory, area) { + area_ha <- area$area_m2 / 10000 + source <- if (nrow(inventory)) unique(inventory$segmentation_source) else attr(inventory, "segmentation_source") + dimension <- if (nrow(inventory)) unique(inventory$instance_dimension) else attr(inventory, "instance_dimension") + data.table::data.table( + segmentation_source = source, + instance_dimension = dimension, + area_source = area$area_source, + area_m2 = area$area_m2, + tree_count = nrow(inventory), + trees_per_ha = if (is.finite(area_ha) && area_ha > 0) nrow(inventory) / area_ha else NA_real_, + basal_area_m2_per_ha = if (is.finite(area_ha) && area_ha > 0 && any(is.finite(inventory$dbh_m))) sum(pi * (inventory$dbh_m / 2)^2, na.rm = TRUE) / area_ha else NA_real_, + mean_dbh_m = safe_stat(inventory$dbh_m, mean), + median_dbh_m = safe_stat(inventory$dbh_m, stats::median), + dbh_p05_m = safe_quantile(inventory$dbh_m, 0.05), + dbh_p25_m = safe_quantile(inventory$dbh_m, 0.25), + dbh_p75_m = safe_quantile(inventory$dbh_m, 0.75), + dbh_p95_m = safe_quantile(inventory$dbh_m, 0.95), + mean_height_m = safe_stat(inventory$height_m, mean), + median_height_m = safe_stat(inventory$height_m, stats::median), + height_p05_m = safe_quantile(inventory$height_m, 0.05), + height_p25_m = safe_quantile(inventory$height_m, 0.25), + height_p75_m = safe_quantile(inventory$height_m, 0.75), + height_p95_m = safe_quantile(inventory$height_m, 0.95) + ) +} + +species_composition <- function(inventories) { + eligible <- Filter(function(inventory) "species_id" %in% names(inventory), inventories) + if (!length(eligible)) return(NULL) + combined <- data.table::rbindlist(eligible, fill = TRUE) + combined <- combined[!is.na(species_id)] + if (!nrow(combined)) { + return(data.table::data.table( + segmentation_source = character(), + instance_dimension = character(), + species_id = numeric(), + tree_count = integer(), + tree_proportion = numeric() + )) + } + result <- combined[, .(tree_count = .N), by = .(segmentation_source, instance_dimension, species_id)] + result[, tree_proportion := tree_count / sum(tree_count), by = .(segmentation_source, instance_dimension)] + result[] +} + +read_supplied_seeds <- function(path) { + if (!file.exists(path)) abort(sprintf("Seed file does not exist: %s", path)) + seeds <- data.table::fread(path) + required <- c("X", "Y", "Z", "TreeID") + missing <- setdiff(required, names(seeds)) + if (length(missing)) abort(sprintf("Seed file is missing columns: %s", paste(missing, collapse = ", "))) + seeds <- as.data.frame(seeds[, ..required]) + if (!nrow(seeds) || anyNA(seeds) || !all(vapply(seeds, is.numeric, logical(1))) || + any(!is.finite(as.matrix(seeds)))) { + abort("Seed columns must contain finite, non-missing numeric values") + } + if (any(seeds$TreeID <= 0) || any(seeds$TreeID != as.integer(seeds$TreeID)) || anyDuplicated(seeds$TreeID)) { + abort("Seed TreeID values must be unique positive integers") + } + seeds +} + +run_csp <- function(las, config, stage_dir) { + seeds <- if (config$seed_mode == "supplied") { + read_supplied_seeds(config$seed_file) + } else { + CspStandSegmentation::find_base_coordinates_raster( + las, + res = config$seed_resolution, + zmin = config$seed_z_min, + zmax = config$seed_z_max, + q = config$seed_quantile, + eps = config$seed_eps + ) + } + if (!nrow(seeds)) abort("CSP seed detection produced no seeds") + if (nrow(seeds) < 2L) abort("CSP requires at least two effective seeds with upstream version 0.2.0") + data.table::fwrite(seeds, file.path(stage_dir, "effective_seeds.tsv"), sep = "\t", na = "") + + weighted <- any(c(config$verticality_weight, config$linearity_weight, config$sphericity_weight) > 0) + working <- las + if (weighted) { + working <- CspStandSegmentation::add_geometry( + working, + k = config$geometry_k, + n_cores = config$geometry_threads + ) + } + working <- lidR::add_lasattribute(working, seq_len(nrow(working@data)), "CSPPointIndex", "Temporary point index") + segmented <- CspStandSegmentation::csp_cost_segmentation( + working, + seeds, + Voxel_size = config$voxel_resolution, + V_w = config$verticality_weight, + L_w = config$linearity_weight, + S_w = config$sphericity_weight, + N_cores = config$routing_workers, + N_trees = 1 + ) + if (!all(c("TreeID", "CSPPointIndex") %in% names(segmented@data))) abort("CSP did not return point assignments") + point_indices <- as.integer(round(segmented@data$CSPPointIndex)) + valid <- !is.na(point_indices) & point_indices >= 1L & point_indices <= nrow(las@data) + if (!any(valid)) abort("CSP returned no assignments that could be mapped to input points") + ids <- integer(nrow(las@data)) + ids[point_indices[valid]] <- as.integer(segmented@data$TreeID[valid]) + list( + seeds = seeds, + ids = ids, + assigned_point_count = sum(ids > 0L), + unassigned_point_count = sum(ids == 0L) + ) +} + +write_csp_cloud <- function(original_las, csp_ids, path) { + if ("PredInstance_CSP" %in% names(original_las@data)) abort("Input already contains PredInstance_CSP") + output <- lidR::add_lasattribute(original_las, csp_ids, "PredInstance_CSP", "CSP tree instance ID") + lidR::writeLAS(output, path) + invisible(path) +} + +process_peak_rss_kb <- function() { + status <- tryCatch(readLines("/proc/self/status", warn = FALSE), error = function(error) character()) + line <- grep("^VmHWM:", status, value = TRUE) + if (!length(line)) return(NA_real_) + as.numeric(sub("^VmHWM:\\s+([0-9]+).*", "\\1", line[[1L]])) +} + +write_outputs <- function(inventories, area, stage_dir) { + inventory_dir <- file.path(stage_dir, "inventories") + dir.create(inventory_dir, recursive = TRUE) + for (name in names(inventories)) { + data.table::fwrite(inventories[[name]], file.path(inventory_dir, paste0(name, ".tsv")), sep = "\t", na = "") + } + combined <- data.table::rbindlist(inventories, use.names = TRUE, fill = TRUE) + has_species <- any(vapply(inventories, function(x) "species_id" %in% names(x), logical(1))) + if (!has_species) { + species_columns <- intersect(c("species_id", "species_prob", "species_conflict_fraction"), names(combined)) + if (length(species_columns)) combined[, (species_columns) := NULL] + } + data.table::fwrite(combined, file.path(stage_dir, "inventory_combined.tsv"), sep = "\t", na = "") + summaries <- data.table::rbindlist(lapply(inventories, stand_summary_for_inventory, area = area), fill = TRUE) + data.table::fwrite(summaries, file.path(stage_dir, "stand_summary.tsv"), sep = "\t", na = "") + composition <- species_composition(inventories) + if (!is.null(composition)) { + data.table::fwrite(composition, file.path(stage_dir, "species_composition.tsv"), sep = "\t", na = "") + } +} + +publish_stage <- function(stage_dir, output_dir) { + if (dir.exists(output_dir) && length(list.files(output_dir, all.files = TRUE, no.. = TRUE))) { + abort(sprintf("Output directory is not empty: %s", output_dir)) + } + dir.create(output_dir, recursive = TRUE, showWarnings = FALSE) + entries <- list.files(stage_dir, all.files = TRUE, no.. = TRUE, full.names = TRUE) + if (!all(file.rename(entries, file.path(output_dir, basename(entries))))) abort("Failed to publish staged outputs") + unlink(stage_dir, recursive = TRUE, force = TRUE) +} + +point_cloud_bounds <- function(files) { + headers <- lapply(files, lidR::readLASheader) + list( + xmin = min(vapply(headers, function(x) x@PHB[["Min X"]], numeric(1))), + xmax = max(vapply(headers, function(x) x@PHB[["Max X"]], numeric(1))), + ymin = min(vapply(headers, function(x) x@PHB[["Min Y"]], numeric(1))), + ymax = max(vapply(headers, function(x) x@PHB[["Max Y"]], numeric(1))), + headers = headers + ) +} + +spatial_chunks <- function(bounds, size) { + x_breaks <- seq(floor(bounds$xmin / size) * size, ceiling(bounds$xmax / size) * size, by = size) + y_breaks <- seq(floor(bounds$ymin / size) * size, ceiling(bounds$ymax / size) * size, by = size) + if (length(x_breaks) < 2L) x_breaks <- c(x_breaks, x_breaks + size) + if (length(y_breaks) < 2L) y_breaks <- c(y_breaks, y_breaks + size) + # Core reads are half-open on their upper edges to avoid double counting. + # Add a final interval when the point-cloud maximum lies exactly on a grid + # line so points on that maximum are still assigned to one chunk. + if (tail(x_breaks, 1L) <= bounds$xmax) x_breaks <- c(x_breaks, tail(x_breaks, 1L) + size) + if (tail(y_breaks, 1L) <= bounds$ymax) y_breaks <- c(y_breaks, tail(y_breaks, 1L) + size) + grid <- expand.grid(x = seq_len(length(x_breaks) - 1L), y = seq_len(length(y_breaks) - 1L)) + lapply(seq_len(nrow(grid)), function(index) { + x <- grid$x[[index]] + y <- grid$y[[index]] + c(xmin = x_breaks[[x]], xmax = x_breaks[[x + 1L]], ymin = y_breaks[[y]], ymax = y_breaks[[y + 1L]]) + }) +} + +read_spatial_chunk <- function(files, selectors, extent, specs = NULL) { + filter <- sprintf( + "-keep_xy %.10f %.10f %.10f %.10f", + extent[["xmin"]], extent[["ymin"]], extent[["xmax"]], extent[["ymax"]] + ) + clouds <- Map(function(file, select) { + cloud <- suppressWarnings(lidR::readLAS(file, select = select, filter = filter)) + if (!is.null(cloud) && !lidR::is.empty(cloud) && !is.null(specs)) { + cloud <- project_inventory_dimensions(cloud, specs) + } + cloud + }, files, selectors) + clouds <- Filter(function(x) !is.null(x) && !lidR::is.empty(x), clouds) + if (!length(clouds)) return(NULL) + if (length(clouds) == 1L) clouds[[1L]] else CspStandSegmentation::las_merge(clouds, fill = TRUE) +} + +parallel_chunk_map <- function(values, workers, fun) { + workers <- min(as.integer(workers), length(values)) + if (workers <= 1L) return(lapply(values, fun)) + results <- parallel::mclapply( + values, + fun, + mc.cores = workers, + mc.preschedule = TRUE, + mc.set.seed = FALSE + ) + failed <- vapply(results, inherits, logical(1), what = "try-error") + if (any(failed)) { + abort(sprintf("Parallel DTM worker failed: %s", as.character(results[[which(failed)[[1L]]]]))) + } + results +} + +chunked_dtm <- function(files, chunks, buffer, resolution, work_dir, stage_dir, workers = 1L) { + partial_dir <- file.path(work_dir, "dtm") + dir.create(partial_dir, recursive = TRUE) + selectors <- rep("c", length(files)) + process_chunk <- function(index) { + lidR::set_lidr_threads(1L) + core <- chunks[[index]] + buffered <- core + c(xmin = -buffer, xmax = buffer, ymin = -buffer, ymax = buffer) + las <- read_spatial_chunk(files, selectors, buffered) + if (is.null(las)) return(NULL) + if (nrow(las@data) < 3L || length(unique(las@data$X)) < 2L || length(unique(las@data$Y)) < 2L) { + return(NULL) + } + result <- make_dtm(las, resolution) + tile <- terra::crop( + result$dtm, + terra::ext(core[["xmin"]], core[["xmax"]], core[["ymin"]], core[["ymax"]]) + ) + if (terra::ncell(tile) == 0L || all(is.na(terra::values(tile)))) return(NULL) + path <- file.path(partial_dir, sprintf("dtm_%05d.tif", index)) + terra::writeRaster(tile, path, overwrite = TRUE) + list(path = path, method = result$method, peak_rss_kb = process_peak_rss_kb()) + } + results <- parallel_chunk_map(as.list(seq_along(chunks)), workers, process_chunk) + results <- Filter(Negate(is.null), results) + partials <- vapply(results, `[[`, character(1), "path") + methods <- vapply(results, `[[`, character(1), "method") + if (!length(partials)) abort("Chunked DTM generation produced no valid cells") + dtm <- terra::vrt(partials) + bounds <- point_cloud_bounds(files) + target <- terra::ext( + floor(bounds$xmin / resolution) * resolution, + ceiling(bounds$xmax / resolution) * resolution, + floor(bounds$ymin / resolution) * resolution, + ceiling(bounds$ymax / resolution) * resolution + ) + dtm <- terra::extend(dtm, target) + output <- file.path(stage_dir, "dtm_full.tif") + terra::writeRaster(dtm, output, overwrite = TRUE) + list( + dtm = terra::rast(output), + method = if (all(methods == "classification_2")) "chunked_classification_2" else "chunked_csf", + chunk_count = length(chunks), + valid_chunk_count = length(partials), + workers = min(as.integer(workers), length(chunks)), + worker_peak_rss_sum_kb = sum(vapply(results, `[[`, numeric(1), "peak_rss_kb")), + worker_peak_rss_max_kb = max(vapply(results, `[[`, numeric(1), "peak_rss_kb")) + ) +} + +header_point_count <- function(header) { + value <- header@PHB[["Number of point records"]] + if (is.null(value)) abort("LAS header does not declare a point count") + as.numeric(value) +} + +streaming_dtm_candidates <- function(files, config, work_dir) { + candidate_dir <- file.path(work_dir, "dtm_candidates") + root <- if (exists("repository_root", inherits = TRUE)) { + get("repository_root", inherits = TRUE) + } else { + normalizePath(getwd(), mustWork = TRUE) + } + script <- file.path(root, "workflow", "stream_dtm_candidates.py") + if (!file.exists(script)) abort(sprintf("Streaming DTM helper is missing: %s", script)) + log_path <- file.path(work_dir, "dtm_candidates.log") + arguments <- c( + unlist(lapply(files, function(path) c("--input", path)), use.names = FALSE), + "--output-dir", candidate_dir, + "--candidate-resolution", as.character(config$dtm_candidate_resolution), + "--output-resolution", as.character(config$dtm_resolution), + "--workers", as.character(config$dtm_workers) + ) + status <- system2("python3", c(script, arguments), stdout = log_path, stderr = log_path) + if (!identical(status, 0L)) { + detail <- if (file.exists(log_path)) paste(tail(readLines(log_path, warn = FALSE), 20L), collapse = "\n") else "" + abort(sprintf("Streaming DTM candidate generation failed%s", if (nzchar(detail)) paste0(":\n", detail) else "")) + } + metadata_path <- file.path(candidate_dir, "candidate_metadata.json") + candidate_path <- file.path(candidate_dir, "minimum_candidates.laz") + if (!file.exists(metadata_path) || !file.exists(candidate_path)) { + abort("Streaming DTM candidate generation did not produce its declared outputs") + } + list( + path = candidate_path, + metadata = jsonlite::fromJSON(metadata_path, simplifyVector = TRUE) + ) +} + +make_inventory_dtm <- function(files, bounds, chunks, config, work_dir, stage_dir) { + input_points <- sum(vapply(bounds$headers, header_point_count, numeric(1))) + strategy <- config$dtm_strategy + if (identical(strategy, "auto")) { + strategy <- if (input_points >= config$dtm_streaming_threshold) "streaming" else "spatial" + } + if (identical(strategy, "spatial")) { + result <- chunked_dtm( + files, chunks, config$chunk_buffer, config$dtm_resolution, + work_dir, stage_dir, config$dtm_workers + ) + result$strategy <- "spatial" + result$candidate_metadata <- NULL + return(result) + } + + candidates <- streaming_dtm_candidates(files, config, work_dir) + candidate_bounds <- point_cloud_bounds(candidates$path) + candidate_chunks <- spatial_chunks(candidate_bounds, config$read_chunk_size) + result <- chunked_dtm( + candidates$path, candidate_chunks, config$chunk_buffer, + config$dtm_resolution, work_dir, stage_dir, config$dtm_workers + ) + result$method <- paste0("streaming_candidates_", result$method) + result$strategy <- "streaming" + result$candidate_metadata <- candidates$metadata + result +} + +update_streaming_hull <- function(hull, x, y) { + candidate <- data.table::data.table(x = x, y = y) + if (nrow(candidate) > 3L) candidate <- candidate[chull(x, y)] + if (!is.null(hull)) candidate <- data.table::rbindlist(list(hull, candidate)) + if (nrow(candidate) > 3L) candidate <- candidate[chull(x, y)] + candidate +} + +append_partition <- function(path, values) { + connection <- file(path, open = "ab") + on.exit(close(connection)) + writeBin(as.double(t(as.matrix(values))), connection, size = 8L) +} + +stream_inventory_partitions <- function(files, chunks, specs, config, work_dir) { + selectors <- vapply(files, read_selector, character(1), specs = specs, preserve_all_dimensions = FALSE) + partition_root <- file.path(work_dir, "partitions") + dir.create(partition_root, recursive = TRUE) + layouts <- list() + for (spec in specs) { + slug <- gsub("[^A-Za-z0-9_]+", "_", spec$source) + columns <- unique(c( + "X", "Y", "Z", spec$instance, spec$species, spec$species_prob, + if (identical(spec$source, "FM")) c("PredScore_FM", "PredSemantic_FM") + )) + columns <- columns[!is.na(columns) & nzchar(columns)] + directory <- file.path(partition_root, slug) + dir.create(directory) + layouts[[slug]] <- list(spec = spec, columns = columns, directory = directory) + } + + point_count <- 0 + hull <- NULL + loaded_dimensions <- character() + dropped_dimensions <- character() + for (core in chunks) { + las <- read_spatial_chunk(files, selectors, core, specs) + if (is.null(las)) next + points <- las@data[ + X >= core[["xmin"]] & X < core[["xmax"]] & + Y >= core[["ymin"]] & Y < core[["ymax"]] + ] + if (!nrow(points)) next + point_count <- point_count + nrow(points) + hull <- update_streaming_hull(hull, points$X, points$Y) + loaded_dimensions <- union(loaded_dimensions, names(points)) + dropped_dimensions <- union(dropped_dimensions, attr(las, "dropped_dimensions")) + + for (layout in layouts) { + spec <- layout$spec + instance <- spec$instance + available_columns <- intersect(layout$columns, names(points)) + selected <- points[ + !is.na(get(instance)) & is.finite(get(instance)) & + get(instance) >= 0 & !(get(instance) %in% config$non_tree_ids), + ..available_columns + ] + if (!nrow(selected)) next + missing <- setdiff(layout$columns, names(selected)) + for (dimension in missing) data.table::set(selected, j = dimension, value = NA_real_) + layout_columns <- layout$columns + selected <- selected[, ..layout_columns] + data.table::set( + selected, + j = "partition__", + value = as.integer(abs(selected[[instance]]) %% config$inventory_partitions) + 1L + ) + for (partition in unique(selected$partition__)) { + path <- file.path(layout$directory, sprintf("%05d.bin", partition)) + partition_columns <- setdiff(names(selected), "partition__") + append_partition(path, selected[partition__ == partition, ..partition_columns]) + } + } + rm(las, points) + gc(FALSE) + } + list( + layouts = layouts, + point_count = point_count, + hull = hull, + loaded_dimensions = loaded_dimensions, + dropped_dimensions = dropped_dimensions, + selectors = selectors + ) +} + +read_binary_partition <- function(path, columns) { + count <- file.info(path)$size / 8 + values <- readBin(path, numeric(), n = count, size = 8L) + if (length(values) %% length(columns) != 0L) abort(sprintf("Invalid inventory partition: %s", path)) + data.table::as.data.table(matrix(values, ncol = length(columns), byrow = TRUE, dimnames = list(NULL, columns))) +} + +inventory_from_partitions <- function(stream, dtm, config, aoi, crs) { + inventories <- list() + for (slug in names(stream$layouts)) { + layout <- stream$layouts[[slug]] + paths <- sort(list.files(layout$directory, pattern = "\\.bin$", full.names = TRUE)) + partials <- list() + for (path in paths) { + points <- read_binary_partition(path, layout$columns) + terrain <- terra::extract(dtm, data.frame(X = points$X, Y = points$Y), ID = FALSE)[[1L]] + data.table::set(points, j = "Zref", value = points$Z) + data.table::set(points, j = "Z", value = points$Z - terrain) + points <- points[is.finite(Z)] + if (!nrow(points)) next + las <- lidR::LAS(as.data.frame(points)) + partials[[length(partials) + 1L]] <- inventory_for_spec(las, layout$spec, config, aoi, crs) + rm(points, las) + gc(FALSE) + } + if (!length(partials)) abort(sprintf("No valid inventory partitions for %s", layout$spec$instance)) + inventory <- data.table::rbindlist(partials, use.names = TRUE, fill = TRUE) + data.table::setorder(inventory, instance_id) + attr(inventory, "segmentation_source") <- layout$spec$source + attr(inventory, "instance_dimension") <- layout$spec$instance + inventories[[slug]] <- inventory + } + inventories +} + +run_chunked_inventory <- function(config, files, stage_dir, output_parent, started) { + work_dir <- tempfile("csp-chunks-", tmpdir = output_parent) + dir.create(work_dir) + on.exit(unlink(work_dir, recursive = TRUE, force = TRUE), add = TRUE) + bounds <- point_cloud_bounds(files) + chunks <- spatial_chunks(bounds, config$read_chunk_size) + available <- unique(c( + "X", "Y", "Z", "Classification", + unlist(lapply(bounds$headers, extra_byte_names), use.names = FALSE) + )) + for (spec in config$segmentation_specs) { + required <- c(spec$instance, spec$species, spec$species_prob) + missing <- setdiff(required[!is.na(required)], available) + if (length(missing)) abort(sprintf("Requested dimensions are missing for %s: %s", spec$instance, paste(missing, collapse = ", "))) + } + las_crs <- sf::st_crs(bounds$headers[[1L]]) + aoi <- read_aoi(config$aoi_json, las_crs) + dtm_result <- make_inventory_dtm(files, bounds, chunks, config, work_dir, stage_dir) + if (!is.null(aoi)) { + aoi_dtm <- terra::mask(terra::crop(dtm_result$dtm, terra::vect(aoi)), terra::vect(aoi)) + terra::writeRaster(aoi_dtm, file.path(stage_dir, "dtm_aoi.tif"), overwrite = TRUE) + } + stream <- stream_inventory_partitions(files, chunks, config$segmentation_specs, config, work_dir) + area <- if (!is.null(aoi)) { + list(area_m2 = as.numeric(sf::st_area(aoi)), area_source = "aoi") + } else { + list(area_m2 = polygon_area(stream$hull$x, stream$hull$y), area_source = "point_cloud_convex_hull") + } + inventories <- inventory_from_partitions(stream, dtm_result$dtm, config, aoi, las_crs) + write_outputs(inventories, area, stage_dir) + + metadata <- list( + tool = "3Dtrees: CSP StandSegmentation", + package_version = as.character(utils::packageVersion("CspStandSegmentation")), + input_files = unname(files), input_bytes = unname(sum(file.info(files)$size)), + input_points = stream$point_count, input_dimensions = stream$loaded_dimensions, + input_read_selectors = unname(stream$selectors), + input_dimensions_dropped = stream$dropped_dimensions, + segmentation_dimensions = vapply(config$segmentation_specs, `[[`, character(1), "instance"), + species_dimensions_supplied = any(vapply(config$segmentation_specs, function(spec) !is.null(spec$species), logical(1))), + csp_enabled = FALSE, csp_seed_count = 0L, csp_unassigned_point_count = 0L, + read_mode = "spatial_chunks_and_instance_partitions", + read_chunk_size_m = config$read_chunk_size, + chunk_buffer_m = config$chunk_buffer, + spatial_chunk_count = length(chunks), + dtm_workers = dtm_result$workers, + dtm_valid_chunk_count = dtm_result$valid_chunk_count, + dtm_strategy = dtm_result$strategy, + dtm_candidate_resolution_m = if (is.null(dtm_result$candidate_metadata)) NULL else dtm_result$candidate_metadata$candidate_resolution_m, + dtm_candidate_xy_mode = if (is.null(dtm_result$candidate_metadata)) NULL else dtm_result$candidate_metadata$candidate_xy_mode, + dtm_candidate_points = if (is.null(dtm_result$candidate_metadata)) NULL else dtm_result$candidate_metadata$candidate_points, + dtm_candidate_elapsed_seconds = if (is.null(dtm_result$candidate_metadata)) NULL else dtm_result$candidate_metadata$elapsed_seconds, + dtm_candidate_average_cpu_cores = if (is.null(dtm_result$candidate_metadata)) NULL else dtm_result$candidate_metadata$average_cpu_cores, + dtm_candidate_worker_peak_rss_sum_upper_bound_kb = if (is.null(dtm_result$candidate_metadata)) NULL else dtm_result$candidate_metadata$worker_peak_rss_sum_upper_bound_kb, + dtm_candidate_parent_peak_rss_kb = if (is.null(dtm_result$candidate_metadata)) NULL else dtm_result$candidate_metadata$parent_peak_rss_kb, + dtm_worker_peak_rss_sum_kb = dtm_result$worker_peak_rss_sum_kb, + dtm_worker_peak_rss_max_kb = dtm_result$worker_peak_rss_max_kb, + inventory_partitions = config$inventory_partitions, + ground_method = dtm_result$method, dtm_resolution_m = config$dtm_resolution, + area_source = area$area_source, area_m2 = area$area_m2, + random_seed = config$random_seed, routing_workers = config$routing_workers, + geometry_threads = config$geometry_threads, voxel_resolution_m = config$voxel_resolution + ) + jsonlite::write_json(metadata, file.path(stage_dir, "run_metadata.json"), auto_unbox = TRUE, pretty = TRUE, na = "null") + completed <- Sys.time() + resources <- list( + elapsed_seconds = as.numeric(difftime(completed, started, units = "secs")), + process_peak_rss_kb = process_peak_rss_kb(), input_bytes = unname(sum(file.info(files)$size)), + input_points = stream$point_count, routing_workers = config$routing_workers, + geometry_threads = config$geometry_threads, + note = "Process-local VmHWM only; container/cgroup telemetry is recorded by benchmark runs." + ) + jsonlite::write_json(resources, file.path(stage_dir, "resource_summary.json"), auto_unbox = TRUE, pretty = TRUE, na = "null") + resources +} + +run_tool <- function(config) { + started <- Sys.time() + set.seed(config$random_seed) + lidR::set_lidr_threads(config$geometry_threads) + files <- resolve_input_files(config$input) + output_parent <- dirname(config$output_dir) + dir.create(output_parent, recursive = TRUE, showWarnings = FALSE) + stage_dir <- tempfile("csp-stage-", tmpdir = output_parent) + dir.create(stage_dir) + published <- FALSE + on.exit(if (!published) unlink(stage_dir, recursive = TRUE, force = TRUE), add = TRUE) + + if (!config$enable_csp) { + resources <- run_chunked_inventory(config, files, stage_dir, output_parent, started) + publish_stage(stage_dir, config$output_dir) + published <- TRUE + message(sprintf("Completed CSP stand inventory in %.1f seconds", resources$elapsed_seconds)) + return(invisible(resources)) + } + + original_las <- read_point_cloud( + files, + config$segmentation_specs, + preserve_all_dimensions = config$enable_csp + ) + input_point_count <- nrow(original_las@data) + input_dimensions <- names(original_las@data) + input_read_selectors <- attr(original_las, "read_selectors") + input_dimensions_dropped <- attr(original_las, "dropped_dimensions") + original_data <- if (config$enable_csp) data.table::copy(original_las@data) else NULL + validate_specs(original_las, config$segmentation_specs) + las_crs <- sf::st_crs(original_las) + aoi <- read_aoi(config$aoi_json, las_crs) + area <- determine_area(original_las, aoi) + + dtm_result <- make_dtm(original_las, config$dtm_resolution) + terra::writeRaster(dtm_result$dtm, file.path(stage_dir, "dtm_full.tif"), overwrite = TRUE) + if (!is.null(aoi)) { + aoi_dtm <- terra::mask(terra::crop(dtm_result$dtm, terra::vect(aoi)), terra::vect(aoi)) + terra::writeRaster(aoi_dtm, file.path(stage_dir, "dtm_aoi.tif"), overwrite = TRUE) + } + normalized <- lidR::normalize_height(dtm_result$las, lidR::tin(), dtm = dtm_result$dtm) + + specs <- config$segmentation_specs + inventory_las <- normalized + csp_seed_count <- 0L + csp_unassigned_point_count <- 0L + if (config$enable_csp) { + csp <- run_csp(normalized, config, stage_dir) + csp_seed_count <- nrow(csp$seeds) + csp_unassigned_point_count <- csp$unassigned_point_count + original_las@data <- original_data + write_csp_cloud(original_las, csp$ids, file.path(stage_dir, "segmented_csp.laz")) + inventory_las <- normalized + inventory_las <- lidR::add_lasattribute(inventory_las, csp$ids, "PredInstance_CSP", "CSP tree instance ID") + specs <- c(specs, list(list(instance = "PredInstance_CSP", source = "CSP", species = NULL, species_prob = NULL))) + } + + validate_specs(inventory_las, specs) + inventories <- list() + for (spec in specs) { + result <- inventory_for_spec(inventory_las, spec, config, aoi, las_crs) + slug <- gsub("[^A-Za-z0-9_]+", "_", spec$source) + inventories[[slug]] <- result + } + write_outputs(inventories, area, stage_dir) + + metadata <- list( + tool = "3Dtrees: CSP StandSegmentation", + package_version = as.character(utils::packageVersion("CspStandSegmentation")), + input_files = unname(files), + input_bytes = unname(sum(file.info(files)$size)), + input_points = input_point_count, + input_dimensions = input_dimensions, + input_read_selectors = unname(input_read_selectors), + input_dimensions_dropped = unname(input_dimensions_dropped), + segmentation_dimensions = vapply(specs, `[[`, character(1), "instance"), + species_dimensions_supplied = any(vapply(specs, function(spec) !is.null(spec$species), logical(1))), + csp_enabled = config$enable_csp, + csp_seed_count = csp_seed_count, + csp_unassigned_point_count = csp_unassigned_point_count, + ground_method = dtm_result$method, + dtm_resolution_m = config$dtm_resolution, + area_source = area$area_source, + area_m2 = area$area_m2, + random_seed = config$random_seed, + routing_workers = config$routing_workers, + geometry_threads = config$geometry_threads, + voxel_resolution_m = config$voxel_resolution + ) + jsonlite::write_json(metadata, file.path(stage_dir, "run_metadata.json"), auto_unbox = TRUE, pretty = TRUE, na = "null") + completed <- Sys.time() + resources <- list( + elapsed_seconds = as.numeric(difftime(completed, started, units = "secs")), + process_peak_rss_kb = process_peak_rss_kb(), + input_bytes = unname(sum(file.info(files)$size)), + input_points = input_point_count, + routing_workers = config$routing_workers, + geometry_threads = config$geometry_threads, + note = "Process-local VmHWM only; container/cgroup telemetry is recorded by benchmark runs." + ) + jsonlite::write_json(resources, file.path(stage_dir, "resource_summary.json"), auto_unbox = TRUE, pretty = TRUE, na = "null") + + publish_stage(stage_dir, config$output_dir) + published <- TRUE + message(sprintf("Completed CSP stand inventory in %.1f seconds", resources$elapsed_seconds)) + invisible(config$output_dir) +}