Fix add_measures join: key divide area by flowpath_id, not divide_id

add_measures joined incremental catchment area onto flowpaths with
`by = c("flowpath_id" = "divide_id")`, assuming the legacy 1:1 convention
where divide_id == flowpath_id. In the current schema these differ
(e.g. cat-* vs fp-*), so the join matched nothing: areasqkm came back NA and
was forced to 0, which then propagated through accumulate_downstream to make
total_dasqkm uniformly 0 in flowpath/network outputs.

Prefer the explicit divide->flowpath link (`flowpath_id` on divides) when
present; fall back to the legacy divide_id join otherwise. Verified on
ngen_vpu09: 22977/22989 flowpaths now receive correct areas (the remaining 12
are genuinely divide-less). Adds regression tests for both join paths.

Author: Mike Johnson

R/utils.R

@@ -47,10 +47,22 @@ add_measures <- function(flowpaths, divides) {
   flowpaths$lengthkm <- add_lengthkm(flowpaths)
   divides$areasqkm <- add_areasqkm(divides)
   flowpaths$areasqkm <- NULL
-  flowpaths <- dplyr::left_join(flowpaths,
-                                dplyr::select(sf::st_drop_geometry(divides), divide_id, areasqkm),
-                                by = c("flowpath_id" = "divide_id")
-  )
+  div_tab <- sf::st_drop_geometry(divides)
+  # Join incremental catchment area onto flowpaths. Prefer the explicit
+  # divide->flowpath link (`flowpath_id` on divides), required by the current
+  # schema where divide_id != flowpath_id (e.g. cat-* vs fp-*); joining on
+  # divide_id there matches nothing and silently zeroes areasqkm (and any
+  # downstream-accumulated total area). Fall back to the legacy 1:1 convention
+  # (divide_id == flowpath_id) when divides carry no flowpath_id.
+  if ("flowpath_id" %in% names(div_tab)) {
+    flowpaths <- dplyr::left_join(
+      flowpaths, dplyr::select(div_tab, flowpath_id, areasqkm),
+      by = "flowpath_id")
+  } else {
+    flowpaths <- dplyr::left_join(
+      flowpaths, dplyr::select(div_tab, divide_id, areasqkm),
+      by = c("flowpath_id" = "divide_id"))
+  }
   list(
     flowpaths = rename_geometry(flowpaths, "geometry"),
     divides   = rename_geometry(divides, "geometry")

tests/testthat/test-add-measures.R

@@ -0,0 +1,38 @@
+test_that("add_measures joins divide area by flowpath_id when divide_id != flowpath_id", {
+  skip_if_not_installed("sf")
+  fp <- sf::st_sf(
+    flowpath_id = c("fp-1", "fp-2"),
+    geometry = sf::st_sfc(
+      sf::st_linestring(matrix(c(0, 0, 1000, 1000), ncol = 2, byrow = TRUE)),
+      sf::st_linestring(matrix(c(1000, 1000, 2000, 2000), ncol = 2, byrow = TRUE))),
+    crs = 5070)
+  # cat-* divides carry an explicit flowpath_id link; divide_id != flowpath_id
+  div <- sf::st_sf(
+    divide_id   = c("cat-1", "cat-2"),
+    flowpath_id = c("fp-1", "fp-2"),
+    geometry = sf::st_sfc(
+      sf::st_polygon(list(matrix(c(0,0, 1000,0, 1000,1000, 0,1000, 0,0), ncol = 2, byrow = TRUE))),
+      sf::st_polygon(list(matrix(c(0,0, 2000,0, 2000,2000, 0,2000, 0,0), ncol = 2, byrow = TRUE)))),
+    crs = 5070)
+  res <- add_measures(fp, div)
+  expect_true(all(res$flowpaths$areasqkm > 0))          # not zeroed by a bad join
+  expect_equal(res$flowpaths$areasqkm[2] / res$flowpaths$areasqkm[1], 4, tolerance = 1e-6)
+})
+
+test_that("add_measures falls back to legacy divide_id == flowpath_id", {
+  skip_if_not_installed("sf")
+  fp <- sf::st_sf(
+    flowpath_id = c("1", "2"),
+    geometry = sf::st_sfc(
+      sf::st_linestring(matrix(c(0, 0, 1000, 1000), ncol = 2, byrow = TRUE)),
+      sf::st_linestring(matrix(c(1000, 1000, 2000, 2000), ncol = 2, byrow = TRUE))),
+    crs = 5070)
+  div <- sf::st_sf(   # no flowpath_id column -> legacy join on divide_id
+    divide_id = c("1", "2"),
+    geometry = sf::st_sfc(
+      sf::st_polygon(list(matrix(c(0,0, 1000,0, 1000,1000, 0,1000, 0,0), ncol = 2, byrow = TRUE))),
+      sf::st_polygon(list(matrix(c(0,0, 2000,0, 2000,2000, 0,2000, 0,0), ncol = 2, byrow = TRUE)))),
+    crs = 5070)
+  res <- add_measures(fp, div)
+  expect_true(all(res$flowpaths$areasqkm > 0))
+})