R/compute_lumi.R

In cnefetools: Access and Analysis of Brazilian CNEFE Address Data

#' Compute land-use mix indicators on a spatial grid
#'
#' @description
#' `compute_lumi()` reads CNEFE records for a given municipality,
#' assigns each address point to spatial units (either H3 hexagonal cells or
#' user-provided polygons), and computes the residential proportion (`p_res`) and land-use mix
#' indices, such as the Entropy Index (`ei`), the Herfindahl-Hirschman Index (`hhi`),
#' the Balance Index (`bal`), the Index of Concentration at Extremes (`ice`), the adapted HHI (`hhi_adp`),
#' and the Bidirectional Global-centered Balance Index (`bgbi`), following the methodology
#' proposed in Pedreira Junior et al. (2025, 2026). The 2026 article introduces the BGBI,
#' and the adapted HHI is documented in the 2025 preprint.
#'
#' @param code_muni Integer. Seven-digit IBGE municipality code.
#' @param year Integer. The CNEFE data year. Currently only 2022 is supported.
#'   Defaults to 2022.
#' @param polygon_type `r lifecycle::badge("deprecated")` The aggregation mode is
#'   now inferred from `polygon`: leave it `NULL` for an H3 grid, or pass an
#'   [`sf::sf`] object for user polygons. Passing `polygon_type` still works and
#'   warns.
#' @param polygon An [`sf::sf`] object with polygon geometries. Supplying it
#'   switches the output from an H3 grid to these polygons. A warning is issued reporting the percentage of
#'   CNEFE points covered by the polygon area. If no CNEFE points fall within
#'   the polygon, an error is raised.
#' @param crs_output The CRS for the output object. Only used when `polygon` is
#'   supplied. Default is `NULL`, which uses the original CRS of
#'   the `polygon` argument. Can be an EPSG code (e.g., 4326, 31983) or any CRS
#'   object accepted by [sf::st_transform()].
#' @param h3_resolution Integer. H3 grid resolution (default: 9). Only used for
#'   the H3 grid, so it is ignored when `polygon` is supplied.
#' @param verbose Logical; if `TRUE`, prints messages and timing information.
#' @param cache Logical. If `TRUE` (default), the downloaded data is stored as
#'   a gzipped CSV in the user cache directory and reused in future calls. If
#'   `FALSE`, a temporary file is used and deleted after the call.
#' @param cache_dir Character. Directory to use for cached downloads. If `NULL`
#'   (default), the `CNEFETOOLS_CACHE_DIR` environment variable is used when it
#'   is set, otherwise [tools::R_user_dir()] with `which = "cache"`. Use this to
#'   point large downloads at a secondary drive or a shared volume.
#' @param backend Character. `"duckdb"` (default) uses DuckDB with the H3
#'   extension, reading the cached gzipped CSV directly. `"r"` computes H3 in R
#'   using h3jsr instead, and needs no DuckDB extension.
#'
#'   `"r"` exists for environments where DuckDB extensions cannot be installed,
#'   such as some restricted computing clusters. It is **not** the lighter
#'   option: it materialises the filtered address table in R memory, so its
#'   footprint grows with the municipality, while DuckDB aggregates in a
#'   streaming fashion and stays nearly flat. On São Paulo (5.7 million
#'   addresses) the measured peak is about 8.4 GB under `"r"` against 0.6 GB
#'   under `"duckdb"`, alongside being roughly 15 times slower.
#'
#'   If the constraint is memory rather than installability, keep the DuckDB
#'   backend and cap it with the `cnefetools.duckdb_config` option instead. See
#'   `?cnefetools` for that option, and the benchmark article at
#'   <https://pedreirajr.github.io/cnefetools/articles/bench_duckdb.html> for
#'   the measurements.
#'
#' @return An [`sf::sf`] object containing:
#' \describe{
#'   \item{When `polygon` is `NULL` (H3 grid):}{
#'     \itemize{
#'       \item `id_hex`: H3 cell identifier
#'       \item `p_res`, `ei`, `hhi`, `bal`, `ice`, `hhi_adp`, `bgbi`: land-use
#'         mix indicators
#'       \item `geometry`: hexagon geometry (CRS 4326)
#'     }
#'   }
#'   \item{When `polygon` is supplied:}{
#'     \itemize{
#'       \item Original columns from `polygon`
#'       \item `p_res`, `ei`, `hhi`, `bal`, `ice`, `hhi_adp`, `bgbi`: land-use
#'         mix indicators
#'       \item `geometry`: polygon geometry (in the original or `crs_output` CRS)
#'     }
#'   }
#' }
#'
#' @details
#' ## Binary land-use classification
#'
#' The indices computed here rest on a binary split. An address is counted as
#' residential when `COD_ESPECIE == 1` (private household), and as
#' non-residential otherwise. This follows the formulation of the indices as
#' published in Pedreira Junior et al. (2026), where the measures are defined
#' and empirically validated on that two-category basis.
#'
#' ## Exclusion of buildings under construction
#'
#' `compute_lumi()` drops records with `COD_ESPECIE == 7` (building under
#' construction or renovation), because such records describe a transitional
#' state rather than a realised land use. Note that [cnefe_counts()] does **not**
#' apply this exclusion and reports these records as `addr_type7`.
#'
#' ## The citywide baseline P
#'
#' Two indices use a citywide residential share P: the `bgbi` index, which is
#' referenced against it, and the Balance Index (`bal`), which uses it through
#' r = P / (1 - P). The other indices are computed entirely within each spatial
#' unit. Two properties of P are worth stating.
#'
#' First, P is computed from CNEFE address-type counts rather than from census
#' population, so it describes the distribution of address types and not the
#' distribution of residents.
#'
#' Second, P is always computed over the full municipality, including when
#' `polygon` is supplied, so it does not adapt to the area the supplied
#' polygons happen to cover. This is intended, as P describes the context the
#' addresses sit in, which is the municipality, and a sub-area of a city is
#' still part of that wider context. A baseline recomputed over the sub-area
#' would measure something different, namely mix relative to the sub-area
#' itself rather than relative to the city.
#'
#' @references
#' Pedreira Junior, J. U.; Louro, T. V.; Assis, L. B. M.; Brito, P. L.;
#' Bomfim, F. G. (2026).
#' BGBI: A citywide-referenced and bidirectional land use mix index for
#' planning and policy evaluation.
#' *Land Use Policy*, 169, 108135.
#' https://doi.org/10.1016/j.landusepol.2026.108135
#'
#' Pedreira Junior, J. U.; Louro, T. V.; Assis, L. B. M.; Brito, P. L. (2025).
#' Measuring land use mix with address-level census data.
#' *engrXiv* preprint. https://engrxiv.org/preprint/view/5975
#' (where the adapted HHI (`hhi_adp`) is documented)
#'
#' Massey, D. S. (2001).
#' The prodigal paradigm returns: ecology comes back to sociology.
#' In A. Booth & A. C. Crouter (Eds.),
#' *Does It Take a Village? Community Effects on Children, Adolescents, and Families*.
#' Lawrence Erlbaum.
#'
#' Song, Y.; Merlin, L.; Rodriguez, D. (2013).
#' Comparing measures of urban land use mix.
#' *Computers, Environment and Urban Systems*, 42, 1--13.
#' https://doi.org/10.1016/j.compenvurbsys.2013.08.001
#'
#' @examples
#' \donttest{
#' # Compute land-use mix indices on H3 hexagons
#' lumi <- compute_lumi(code_muni = 2929057, cache = FALSE)
#'
#' # Compute land-use mix indices on user-provided polygons (neighborhoods of Lauro de Freitas-BA)
#' # Using geobr to download neighborhood boundaries
#' library(geobr)
#' nei_ldf <- subset(
#'   read_neighborhood(year = 2022),
#'   code_muni == 2919207
#' )
#' lumi_poly <- compute_lumi(
#'   code_muni = 2919207,
#'   polygon = nei_ldf,
#'   cache = FALSE
#' )
#' }
#'
#' @export
compute_lumi <- function(
  code_muni,
  year = 2022,
  polygon_type = lifecycle::deprecated(),
  polygon = NULL,
  crs_output = NULL,
  h3_resolution = 9,
  verbose = TRUE,
  cache = TRUE,
  cache_dir = NULL,
  backend = c("duckdb", "r")
) {
  polygon_type <- .resolve_polygon_mode(polygon, polygon_type, fn = "compute_lumi")
  backend <- match.arg(backend)
  code_muni <- .normalize_code_muni(code_muni)
  year <- .validate_year(year)

  # Validate polygon argument
  if (polygon_type == "user") {
    .validate_polygon_arg(polygon, crs_output = crs_output)
  }

  # Get the appropriate index for the requested year
  cnefe_index <- .get_cnefe_index(year)

  # Name (optional)
  info <- cnefe_index[
    cnefe_index$code_muni == code_muni,
    ,
    drop = FALSE
  ]
  city_name <- if (
    nrow(info) > 0 && "name_muni" %in% names(info) && !is.na(info$name_muni[1])
  ) {
    info$name_muni[1]
  } else {
    as.character(code_muni)
  }

  if (verbose) {
    cli::cli_alert_info("Processing municipality code {.val {code_muni}}...")
  }

  # We will return sf objects
  rlang::check_installed(
    "sf",
    reason = "to return an sf grid in `compute_lumi()`."
  )

  # Branch: H3 grid vs user-provided polygon
  if (polygon_type == "hex") {
    out <- .compute_lumi_h3(
      code_muni = code_muni,
      year = year,
      h3_resolution = h3_resolution,
      backend = backend,
      cnefe_index = cnefe_index,
      verbose = verbose,
      cache = cache,
      cache_dir = cache_dir
    )
  } else {
    out <- .compute_lumi_user_poly(
      code_muni = code_muni,
      year = year,
      polygon = polygon,
      crs_output = crs_output,
      backend = backend,
      cnefe_index = cnefe_index,
      verbose = verbose,
      cache = cache,
      cache_dir = cache_dir
    )
  }

  return(out)
}


# -----------------------------------------------------------------------------
# Internal: Shared LUMI index computation helpers
# -----------------------------------------------------------------------------

# Balance index (BAL)
.bal_fun <- function(p, P) {
  num <- abs(p - (P / (1 - P)) * (1 - p))
  den <- p + (P / (1 - P)) * (1 - p)
  ifelse(is.na(num) | is.na(den), NA_real_, 1 - num / den)
}

# BGBI function
.bgbi_fun <- function(p, P) {
  num <- (2 * p - 1) - (2 * P - 1)
  den <- 1 - (2 * p - 1) * (2 * P - 1)
  den[den == 0] <- NA_real_
  ifelse(is.na(num) | is.na(den), NA_real_, num / den)
}

# Numerically safe entropy term: treat 0*log(0) as 0
.safe_plogp <- function(x) {
  ifelse(is.na(x) | x <= 0, 0, x * log(x))
}

#' Compute LUMI indices from n_res and n_tot columns
#'
#' Adds p_res, ei, hhi, bal, ice, hhi_adp, bgbi columns to a data frame.
#'
#' @param df A data frame with `n_res` and `n_tot` columns.
#' @param P Global residential proportion for the municipality.
#' @return The data frame with LUMI columns added.
#' @keywords internal
#' @noRd
.compute_lumi_indices <- function(df, P) {
  df |>
    dplyr::mutate(
      p_res = dplyr::if_else(
        .data$n_tot > 0,
        .data$n_res / .data$n_tot,
        NA_real_
      ),
      q_rest = dplyr::if_else(!is.na(.data$p_res), 1 - .data$p_res, NA_real_),

      # EI (k=2)
      ei = dplyr::if_else(
        !is.na(.data$p_res),
        -(.safe_plogp(.data$p_res) + .safe_plogp(.data$q_rest)) / log(2),
        NA_real_
      ),

      # HHI (2 categories)
      hhi = dplyr::if_else(
        !is.na(.data$p_res),
        (.data$p_res^2 + .data$q_rest^2),
        NA_real_
      ),

      # BAL
      bal = dplyr::if_else(
        !is.na(.data$p_res),
        .bal_fun(.data$p_res, P),
        NA_real_
      ),

      # scaled HHI (min=0.5 for k=2)
      hhi_sc = dplyr::if_else(
        !is.na(.data$hhi),
        (.data$hhi - 0.5) / (1 - 0.5),
        NA_real_
      ),

      # ICE (Index of Concentration at Extremes)
      ice = dplyr::if_else(
        !is.na(.data$p_res),
        .data$p_res - .data$q_rest,
        NA_real_
      ),

      hhi_adp = dplyr::if_else(
        !is.na(.data$p_res) & !is.na(.data$hhi_sc),
        sign(.data$p_res - .data$q_rest) * .data$hhi_sc,
        NA_real_
      ),

      bgbi = dplyr::if_else(
        !is.na(.data$p_res),
        .bgbi_fun(.data$p_res, P),
        NA_real_
      )
    )
}


# -----------------------------------------------------------------------------
# Internal: H3 grid aggregation (original behavior)
# -----------------------------------------------------------------------------
.compute_lumi_h3 <- function(
  code_muni,
  year,
  h3_resolution,
  backend,
  cnefe_index,
  verbose,
  cache = TRUE,
  cache_dir = NULL
) {
  # ---------------------------------------------------------------------------
  # Step 1/3: Ensure the cached data file exists
  # ---------------------------------------------------------------------------
  if (verbose) {
    cli::cli_progress_step("Step 1/3: Ensuring the CNEFE data file...",
                           msg_done = "Step 1/3 (CNEFE data ready)")
  }

  zip_info <- .cnefe_ensure_zip(
    code_muni = code_muni,
    index = cnefe_index,
    cache = cache,
    cache_dir = cache_dir,
    year = year,
    verbose = verbose,
    retry_timeouts = c(300L, 600L, 1800L)
  )
  zip_path <- zip_info$zip_path


  if (verbose) {
    cli::cli_progress_done("Step 1/3: Ensuring the CNEFE data file...")
  }

  # ---------------------------------------------------------------------------
  # Step 2/3: Aggregate counts per hex (n_res, n_tot)
  # ---------------------------------------------------------------------------
  if (verbose) {
    cli::cli_progress_step("Step 2/3: Counting addresses per H3 cell...",
                           msg_done = "Step 2/3 (Addresses counted)")
  }

  counts_hex <- NULL

  if (identical(backend, "duckdb")) {
    rlang::check_installed(
      "DBI",
      reason = "to use backend = 'duckdb' in `compute_lumi()`."
    )
    con <- .duckdb_connect(
      extensions = "h3",
      reason = "to use backend = 'duckdb' in `compute_lumi()`.",
      fallback = "backend = \"r\"",
      verbose = verbose
    )

    {
      {
        src <- .cnefe_csv_uri(zip_path)
        if (isTRUE(src$needs_zipfs)) {
          # A cache written by an older version is still a ZIP.
          .duckdb_quiet(.duckdb_ensure_extension(con, "zipfs", verbose = verbose))
        }
        uri <- src$uri
        uri_sql <- gsub("'", "''", uri)

        # only keep the columns needed; exclude COD_ESPECIE == 7
        sql <- sprintf(
          "
        WITH src AS (
          SELECT
            CAST(LONGITUDE AS DOUBLE) AS lon,
            CAST(LATITUDE  AS DOUBLE) AS lat,
            try_cast(COD_ESPECIE AS INTEGER) AS cod
          FROM read_csv_auto('%s', delim=';', header=true, strict_mode=false)
        ),
        filtered AS (
          SELECT
            lower(hex(CAST(h3_latlng_to_cell(lat, lon, %d) AS UBIGINT))) AS id_hex,
            cod
          FROM src
          WHERE
            lon IS NOT NULL AND lat IS NOT NULL
            AND cod BETWEEN 1 AND 8
            AND cod != 7
        )
        SELECT
          id_hex,
          SUM(CASE WHEN cod = 1 THEN 1 ELSE 0 END)::BIGINT AS n_res,
          COUNT(*)::BIGINT AS n_tot
        FROM filtered
        GROUP BY 1;
      ",
          uri_sql,
          as.integer(h3_resolution)
        )

        counts_hex <- .duckdb_quiet(DBI::dbGetQuery(con, sql)) |>
          dplyr::as_tibble() |>
          dplyr::mutate(
            id_hex = as.character(.data$id_hex),
            n_res = as.integer(.data$n_res),
            n_tot = as.integer(.data$n_tot)
          )
      }
    }

  } else {
    # backend "r"
    tab <- read_cnefe(
      code_muni = code_muni,
      year = year,
      output = "arrow",
      cache = cache,
      cache_dir = cache_dir,
      verbose = FALSE
    )

    # Verbs are pushed down to the Arrow table and collected last (#80 R1.10).
    # as.data.frame() first would materialise all 34 columns as an R data frame
    # before three of them are kept. Measured on Fortaleza, 1.19M rows: peak
    # memory falls from 108.4 MB to 68.3 MB.
    df <- tab |>
      dplyr::transmute(
        LONGITUDE = as.numeric(.data$LONGITUDE),
        LATITUDE = as.numeric(.data$LATITUDE),
        COD_ESPECIE = as.integer(.data$COD_ESPECIE)
      ) |>
      dplyr::filter(
        !is.na(.data$LONGITUDE),
        !is.na(.data$LATITUDE),
        !is.na(.data$COD_ESPECIE),
        .data$COD_ESPECIE %in% 1L:8L,
        .data$COD_ESPECIE != 7L
      ) |>
      dplyr::collect()

    if (nrow(df) == 0L) {
      if (verbose) {
        cli::cli_alert_warning(
          "No valid CNEFE points after filtering (COD_ESPECIE 1:8, excluding 7). Returning an empty {.cls sf}."
        )
      }
      return(.empty_lumi_sf())
    }

    coords <- df |>
      dplyr::transmute(lon = .data$LONGITUDE, lat = .data$LATITUDE)

    id_hex <- suppressMessages(h3jsr::point_to_cell(
      coords,
      res = h3_resolution,
      simple = TRUE
    ))

    counts_hex <- df |>
      dplyr::mutate(id_hex = as.character(id_hex)) |>
      dplyr::filter(!is.na(.data$id_hex)) |>
      dplyr::group_by(.data$id_hex) |>
      dplyr::summarise(
        n_res = sum(.data$COD_ESPECIE == 1L, na.rm = TRUE),
        n_tot = dplyr::n(),
        .groups = "drop"
      ) |>
      dplyr::mutate(
        n_res = as.integer(.data$n_res),
        n_tot = as.integer(.data$n_tot)
      )
  }

  if (verbose) {
    cli::cli_progress_done("Step 2/3: Counting addresses per H3 cell...")
  }

  if (is.null(counts_hex) || nrow(counts_hex) == 0L) {
    if (verbose) {
      cli::cli_alert_warning("No hexagons found after aggregation. Returning an empty {.cls sf}.")
    }
    return(.empty_lumi_sf())
  }

  # ---------------------------------------------------------------------------
  # Step 3/3: Build H3 grid from ids + compute indices
  # ---------------------------------------------------------------------------
  if (verbose) {
    cli::cli_progress_step("Step 3/3: Building grid and computing LUMI...",
                           msg_done = "Step 3/3 (Land use mix indices computed)")
  }

  hex_grid <- build_h3_grid(
    h3_resolution = h3_resolution,
    code_muni     = code_muni,
    year          = year
  )

  # Global city residential proportion P (exclude COD_ESPECIE == 7 already)
  P <- sum(counts_hex$n_res, na.rm = TRUE) / sum(counts_hex$n_tot, na.rm = TRUE)

  out <- hex_grid |>
    dplyr::left_join(counts_hex, by = "id_hex") |>
    dplyr::mutate(
      n_res = dplyr::coalesce(as.integer(.data$n_res), 0L),
      n_tot = dplyr::coalesce(as.integer(.data$n_tot), 0L)
    ) |>
    .compute_lumi_indices(P) |>
    dplyr::select(
      "id_hex",
      "p_res",
      "ei",
      "hhi",
      "bal",
      "ice",
      "hhi_adp",
      "bgbi",
      "geometry"
    )

  if (verbose) {
    cli::cli_progress_done("Step 3/3: Building grid and computing LUMI...")
  }

  return(out)
}


# -----------------------------------------------------------------------------
# Internal: User-provided polygon aggregation
# -----------------------------------------------------------------------------
.compute_lumi_user_poly <- function(
  code_muni,
  year,
  polygon,
  crs_output,
  backend,
  cnefe_index,
  verbose,
  cache = TRUE,
  cache_dir = NULL
) {
  # ---------------------------------------------------------------------------
  # Step 1/3: Ensure the cached data file exists and prepare polygon
  # ---------------------------------------------------------------------------
  if (verbose) {
    cli::cli_progress_step("Step 1/3: Ensuring data and preparing polygon...",
                           msg_done = "Step 1/3 (Data and polygon ready)")
  }

  zip_info <- .cnefe_ensure_zip(
    code_muni = code_muni,
    index = cnefe_index,
    cache = cache,
    cache_dir = cache_dir,
    year = year,
    verbose = verbose,
    retry_timeouts = c(300L, 600L, 1800L)
  )
  zip_path <- zip_info$zip_path

  # Store original CRS for output transformation
  original_crs <- sf::st_crs(polygon)

  # Determine output CRS: use crs_output if provided, otherwise use original
  if (is.null(crs_output)) {
    output_crs <- original_crs
  } else {
    output_crs <- sf::st_crs(crs_output)
  }

  # Fix invalid geometries before any spatial operation
  polygon <- sf::st_make_valid(polygon)

  # Transform polygon to WGS84 internally for spatial join with CNEFE points
  polygon_4326 <- sf::st_transform(polygon, 4326)

  # Add row ID for joining
  polygon_4326 <- polygon_4326 |>
    dplyr::mutate(.poly_row_id = dplyr::row_number())

  if (verbose) {
    cli::cli_progress_done("Step 1/3: Ensuring data and preparing polygon...")
  }

  # ---------------------------------------------------------------------------
  # Step 2/3: Read CNEFE points and perform spatial join
  # ---------------------------------------------------------------------------
  if (verbose) {
    cli::cli_progress_step("Step 2/3: Counting addresses per polygon...",
                           msg_done = "Step 2/3 (Addresses counted)")
  }

  if (identical(backend, "duckdb")) {
    rlang::check_installed(
      "DBI",
      reason = "to use backend = 'duckdb' in `compute_lumi()`."
    )
    rlang::check_installed(
      "duckdb",
      reason = "to use backend = 'duckdb' in `compute_lumi()`."
    )

    join_result <- .compute_lumi_user_poly_duckdb(
      zip_path = zip_path,
      polygon = polygon_4326,
      verbose = verbose
    )
  } else {
    join_result <- .compute_lumi_user_poly_r(
      code_muni = code_muni,
      year = year,
      polygon = polygon_4326,
      verbose = verbose,
      cache = cache,
      cache_dir = cache_dir
    )
  }

  if (verbose) {
    cli::cli_progress_done("Step 2/3: Counting addresses per polygon...")
  }

  # ---------------------------------------------------------------------------
  # Step 3/3: Report coverage, compute LUMI indices, and join back to polygon
  # ---------------------------------------------------------------------------
  if (verbose) {
    cli::cli_progress_step("Step 3/3: Computing land use mix indices...",
                           msg_done = "Step 3/3 (Land use mix indices computed)")
  }

  # Extract coverage statistics
  total_points <- join_result$total_points
  points_matched <- join_result$points_matched
  points_outside <- join_result$points_outside
  counts <- join_result$counts
  total_n_res <- join_result$total_n_res
  total_n_tot <- join_result$total_n_tot

  # Check for zero coverage
  if (points_matched == 0L) {
    cli::cli_abort(c(
      "No CNEFE coordinates were captured within the provided polygon.",
      "i" = "This may indicate that:",
      "*" = "The municipality code {.val {code_muni}} does not correspond to the polygon's municipality, or",
      "*" = "The polygon is not located within municipality {.val {code_muni}}."
    ))
  }

  # Calculate coverage percentage
  coverage_pct <- (points_matched / total_points) * 100

  # Issue warning about coverage
  if (points_outside > 0L) {
    cli::cli_warn(c(
      "Polygon coverage: {.val {sprintf('%.1f', coverage_pct)}%} of CNEFE points captured.",
      "i" = "{.val {points_matched}} of {.val {total_points}} points are within the provided polygon.",
      "i" = "{.val {points_outside}} points fell outside the polygon and were not counted."
    ))
  } else {
    cli::cli_inform(c(
      "v" = "All {.val {total_points}} CNEFE points were captured within the provided polygon."
    ))
  }

  # Global P from all municipality CNEFE points (not just those in polygons)
  P <- total_n_res / total_n_tot

  # Join counts to polygon and compute LUMI indices
  if (nrow(counts) == 0L || all(is.na(counts$.poly_row_id))) {
    # No matched points: add empty LUMI columns
    out <- polygon_4326 |>
      dplyr::mutate(
        n_res = 0L,
        n_tot = 0L
      ) |>
      .compute_lumi_indices(P) |>
      dplyr::select(
        -"n_res", -"n_tot", -"q_rest", -"hhi_sc",
        -".poly_row_id"
      )
  } else {
    counts_per_poly <- counts |>
      dplyr::filter(!is.na(.data$.poly_row_id))

    out <- polygon_4326 |>
      dplyr::left_join(counts_per_poly, by = ".poly_row_id") |>
      dplyr::mutate(
        n_res = dplyr::coalesce(as.integer(.data$n_res), 0L),
        n_tot = dplyr::coalesce(as.integer(.data$n_tot), 0L)
      ) |>
      .compute_lumi_indices(P) |>
      dplyr::select(
        -"n_res", -"n_tot", -"q_rest", -"hhi_sc",
        -".poly_row_id"
      )
  }

  # Transform to output CRS
  out <- sf::st_transform(out, output_crs)

  if (verbose) {
    cli::cli_progress_done("Step 3/3: Computing land use mix indices...")
  }

  return(out)
}


# -----------------------------------------------------------------------------
# Internal: DuckDB backend for user polygon spatial join (compute_lumi)
# -----------------------------------------------------------------------------
.compute_lumi_user_poly_duckdb <- function(
  zip_path,
  polygon,
  verbose
) {
  res <- NULL

  con <- .duckdb_connect(
      extensions = "spatial",
    reason = "to use backend = 'duckdb' in `compute_lumi()`.",
    fallback = "backend = \"r\"",
    verbose = verbose
  )

  .duckdb_quiet({
    {
      src <- .cnefe_csv_uri(zip_path)
      if (isTRUE(src$needs_zipfs)) {
        # A cache written by an older version is still a ZIP.
        .duckdb_quiet(.duckdb_ensure_extension(con, "zipfs", verbose = verbose))
      }
      uri <- src$uri
      uri_sql <- gsub("'", "''", uri)

      DBI::dbExecute(con, sprintf(
        "
      CREATE TABLE cnefe_pts AS
      SELECT
        ROW_NUMBER() OVER () AS pt_id,
        CAST(LONGITUDE AS DOUBLE) AS lon,
        CAST(LATITUDE  AS DOUBLE) AS lat,
        try_cast(COD_ESPECIE AS INTEGER) AS COD_ESPECIE,
        ST_Point(CAST(LONGITUDE AS DOUBLE), CAST(LATITUDE AS DOUBLE)) AS geom
      FROM read_csv_auto('%s', delim=';', header=true, strict_mode=false)
      WHERE
        LONGITUDE IS NOT NULL AND LATITUDE IS NOT NULL
        AND try_cast(COD_ESPECIE AS INTEGER) BETWEEN 1 AND 8
        AND try_cast(COD_ESPECIE AS INTEGER) != 7;
      ",
        uri_sql
      ))

      total_points <- DBI::dbGetQuery(
        con,
        "SELECT COUNT(*) AS n FROM cnefe_pts;"
      )$n[1]

  if (total_points == 0L) {
    return(list(
      counts = dplyr::tibble(
        .poly_row_id = integer(0),
        n_res = integer(0),
        n_tot = integer(0)
      ),
      total_points = 0L,
      points_matched = 0L,
      points_outside = 0L,
      total_n_res = 0L,
      total_n_tot = 0L
    ))
  }

  # Global totals from all municipality points (computed inside DuckDB)
  global_totals <- DBI::dbGetQuery(con,
    "
    SELECT
      SUM(CASE WHEN COD_ESPECIE = 1 THEN 1 ELSE 0 END)::INTEGER AS total_n_res,
      COUNT(*)::INTEGER AS total_n_tot
    FROM cnefe_pts;
    "
  )
  total_n_res <- as.integer(global_totals$total_n_res[1])
  total_n_tot <- as.integer(global_totals$total_n_tot[1])

  # Write user polygon to DuckDB via duckspatial
  duckspatial::ddbs_write_table(
    conn = con,
    # Normalize geometry column to "geom"; duckspatial preserves the
    # input sf geometry name, but the SQL below hardcodes "geom".
    data = sf::st_set_geometry(polygon[, ".poly_row_id"], "geom"),
    name = "user_polygons",
    overwrite = TRUE
  )

  # duckspatial 1.0.0 (DuckDB 1.5+) writes GEOMETRY with embedded CRS metadata,
  # which RTREE does not accept. Strip CRS via WKB round-trip first.
  DBI::dbExecute(
    con,
    "ALTER TABLE user_polygons ALTER COLUMN geom SET DATA TYPE GEOMETRY
     USING ST_GeomFromWKB(ST_AsWKB(geom));"
  )

  # Spatial index on user polygons for faster joins
  DBI::dbExecute(
    con,
    "CREATE INDEX IF NOT EXISTS poly_geom_idx ON user_polygons USING RTREE (geom);"
  )

  # Spatial join in DuckDB via ST_Within (LEFT JOIN to track coverage)
  DBI::dbExecute(con,
                 "
      CREATE TABLE joined AS
      SELECT
        p.pt_id,
        p.COD_ESPECIE,
        u.\".poly_row_id\" AS poly_row_id
      FROM cnefe_pts p
      LEFT JOIN user_polygons u
        ON ST_Within(p.geom, u.geom);
      "
  )

  # Coverage stats (computed inside DuckDB)
  coverage <- DBI::dbGetQuery(con,
    "
    SELECT
      COUNT(DISTINCT CASE WHEN poly_row_id IS NOT NULL THEN pt_id END) AS matched
    FROM joined;
    "
  )
  unique_pts_matched <- as.integer(coverage$matched[1])
  points_outside <- total_points - unique_pts_matched

  # Aggregate n_res and n_tot per polygon (inside DuckDB)
  counts <- DBI::dbGetQuery(con,
    "
    SELECT
      poly_row_id AS \".poly_row_id\",
      SUM(CASE WHEN COD_ESPECIE = 1 THEN 1 ELSE 0 END)::INTEGER AS n_res,
      COUNT(*)::INTEGER AS n_tot
    FROM joined
    WHERE poly_row_id IS NOT NULL
    GROUP BY poly_row_id;
    "
  ) |>
    dplyr::as_tibble() |>
    dplyr::mutate(
      n_res = as.integer(.data$n_res),
      n_tot = as.integer(.data$n_tot)
    )

  res <- list(
    counts = counts,
    total_points = total_points,
    points_matched = unique_pts_matched,
    points_outside = points_outside,
    total_n_res = total_n_res,
    total_n_tot = total_n_tot
  )

    }
  })

  return(res)

}


# -----------------------------------------------------------------------------
# Internal: R backend for user polygon spatial join (compute_lumi)
# -----------------------------------------------------------------------------
.compute_lumi_user_poly_r <- function(
  code_muni,
  year,
  polygon,
  verbose,
  cache = TRUE,
  cache_dir = NULL
) {
  # Read CNEFE data via Arrow
  tab <- read_cnefe(
    code_muni = code_muni,
    year = year,
    output = "arrow",
    cache = cache,
    cache_dir = cache_dir,
    verbose = FALSE
  )

  # Verbs are pushed down to the Arrow table and collected last (#80 R1.10).
  # as.data.frame() first would materialise all 34 columns as an R data frame
  # before three of them are kept. Measured on Fortaleza, 1.19M rows: peak
  # memory falls from 108.4 MB to 68.3 MB.
  df <- tab |>
    dplyr::transmute(
      LONGITUDE = as.numeric(.data$LONGITUDE),
      LATITUDE = as.numeric(.data$LATITUDE),
      COD_ESPECIE = as.integer(.data$COD_ESPECIE)
    ) |>
    dplyr::filter(
      !is.na(.data$LONGITUDE),
      !is.na(.data$LATITUDE),
      !is.na(.data$COD_ESPECIE),
      .data$COD_ESPECIE %in% 1L:8L,
      .data$COD_ESPECIE != 7L
    ) |>
    dplyr::collect()

  total_points <- nrow(df)

  if (total_points == 0L) {
    return(list(
      counts = dplyr::tibble(
        .poly_row_id = integer(0),
        n_res = integer(0),
        n_tot = integer(0)
      ),
      total_points = 0L,
      points_matched = 0L,
      points_outside = 0L,
      total_n_res = 0L,
      total_n_tot = 0L
    ))
  }

  # Global totals from all municipality points
  total_n_res <- as.integer(sum(df$COD_ESPECIE == 1L, na.rm = TRUE))
  total_n_tot <- as.integer(nrow(df))

  # Convert to sf points and add point ID
  cnefe_pts <- sf::st_as_sf(
    df,
    coords = c("LONGITUDE", "LATITUDE"),
    crs = 4326
  )
  cnefe_pts$.pt_id <- seq_len(nrow(cnefe_pts))

  # Spatial join
  joined <- sf::st_join(cnefe_pts, polygon[, ".poly_row_id"], join = sf::st_within)

  # Count unique points matched
  unique_pts_matched <- length(unique(joined$.pt_id[!is.na(joined$.poly_row_id)]))
  points_outside <- total_points - unique_pts_matched

  # Aggregate n_res and n_tot per polygon
  counts <- joined |>
    sf::st_drop_geometry() |>
    dplyr::filter(!is.na(.data$.poly_row_id)) |>
    dplyr::group_by(.data$.poly_row_id) |>
    dplyr::summarise(
      n_res = as.integer(sum(.data$COD_ESPECIE == 1L, na.rm = TRUE)),
      n_tot = as.integer(dplyr::n()),
      .groups = "drop"
    )

  return(list(
    counts = counts,
    total_points = total_points,
    points_matched = unique_pts_matched,
    points_outside = points_outside,
    total_n_res = total_n_res,
    total_n_tot = total_n_tot
  ))
}