project_3.R

# =============================================================================
# Project 3: Socioeconomic Analysis & Intervention Planning
# =============================================================================
#
# This script demonstrates how bitfield-encoded outputs from independent
# analyses (projects 1 and 2) can be decoded and recombined into new derived
# metrics for decision support. It is a SIMULATION designed to showcase the
# bitfield approach and produce plausible visuals — not a validated
# socioeconomic model. The heuristics, weights, and thresholds are chosen to
# be realistic enough to illustrate cross-workflow metadata preservation, but
# should not be used for actual policy decisions.
#
# Workflow:
#   1. Load geophysical basis data (human footprint, land cover)
#   2. Decode bitfields from project 1 (livestock) and project 2 (ecology)
#   3. Simulate market distortion types and magnitudes
#   4. Derive GDP adjustment estimates from ecological and economic signals
#   5. Calculate an economic-ecological misalignment index
#   6. Classify system trajectories (temporal dimension)
#   7. Prioritize interventions and identify cross-sectoral synergies
#   8. Encode everything into a 32-bit bitfield
#
# Study area: Black Forest region, SW Germany (~7.6-9.3°E, 47.5-48.8°N)
# =============================================================================

library(terra)
library(geodata)
library(tidyverse)
library(bitfield)

source("functions.R")

set.seed(42)

xmin <- 7.6
xmax <- 9.3
ymin <- 47.5
ymax <- 48.8
area_extent <- ext(xmin, xmax, ymin, ymax)


# 1. get (bio)geophysical basis ----
#
## get footprint
footprint <- footprint(year = 2009, path = tempdir()) |>
  crop(area_extent)
names(footprint) <- "human_footprint"

## get landcover data (not encoded in bitfields, so load fresh)
lc_meta <- data.frame(id = c(1, 2, 4:9),
                      var = c("trees", "grassland", "shrubs", "cropland", "built",
                              "bare", "water", "wetland"))

landcover <- list()
landcover$grassland <- landcover(var = "grassland", path = tempdir()) |>
  crop(area_extent) |> unlist()
landcover$cropland <- landcover(var = "cropland", path = tempdir()) |>
  crop(area_extent) |> unlist()
landcover$trees <- landcover(var = "trees", path = tempdir()) |>
  crop(area_extent) |> unlist()
landcover$water <- landcover(var = "water", path = tempdir()) |>
  crop(area_extent) |> unlist()


# 2. decode bitfield from project_1 to retrieve layers ----
#
lstRegistry <- readRDS("lstRegistry.rds")
lstBitfield <- rast("lstBitfield.tif")

decoded_lst <- bf_decode(x = lstBitfield,
                         registry = lstRegistry,
                         verbose = FALSE)

animals_median <- decoded_lst$`Median livestock density`
animals_sd <- decoded_lst$`Standard deviation`

animals_mean <- rast("animals_mean.tif")


# 3. decode bitfield from project_2 to retrieve layers ----
#
ecoRegistry <- readRDS("eco_reg.rds")
ecoBitfield <- rast("ecoBitfield.tif")

decoded_eco <- bf_decode(x = ecoBitfield,
                         registry = ecoRegistry,
                         verbose = FALSE)

carrying_capacity <- decoded_eco$`Ecological carrying capacity`
exceedance_risk <- decoded_eco$`Exceedance risk`
exceedance <- decoded_eco$`Exceedance magnitude`
deficit_type <- decoded_eco$`Resource limitation type`
deficit_magnitude <- decoded_eco$`Resource limitation magnitude`

plot(c(carrying_capacity, exceedance_risk, exceedance, deficit_type, deficit_magnitude))

# 4. simulate market distortions ----
## Distortion type ----
# 0=none, 1=subsidy, 2=tax, 3=regulation, 4=externality, 5=monopoly, 6=information, 7=public good
distortion_type <- footprint

# determine landscape context
is_agriculture <- landcover$cropland > 0.3
is_forest <- landcover$trees > 0.3
is_grassland <- landcover$grassland > 0.3
is_urban_proximate <- footprint > 15
is_protected_area <- footprint < 8 & (landcover$trees > 0.5 | landcover$water > 0.2)
is_remote <- footprint < 12 & !is_urban_proximate

# assign values based on heuristics
values(distortion_type) <- 0                                  # no distortion (base state)
distortion_type[is_agriculture & !is_protected_area] <- 1     # subsidy, common in productive agricultural areas
distortion_type[is_protected_area & is_urban_proximate] <- 2  # tax, sensitive transition zones
distortion_type[is_forest & is_urban_proximate] <- 3          # regulation, natural areas with moderate human activity
distortion_type[is_agriculture & is_urban_proximate &
                  !is_protected_area] <- 4                    # externality, intensive use areas without compensating mechanisms
distortion_type[is_remote & is_agriculture] <- 5              # monopoly, isolated areas with concentrated economic activity
distortion_type[is_grassland & is_forest] <- 6                # information failure, transition areas between different land use types
distortion_type[is_protected_area & !is_urban_proximate] <- 7 # public good, providing ecosystem services

names(distortion_type) <- "distortion_type"
levels(distortion_type) <- data.frame(id = 0:7,
                                      distortion = c("none", "subsidy", "tax",
                                                     "regulation", "externality",
                                                     "monopoly", "information",
                                                     "public good"))

## distortion magnitude ----
distortion_magnitude <- footprint
values(distortion_magnitude) <- 0

# normalize to 0-1 scale
footprint_norm <- (footprint - min(values(footprint), na.rm=TRUE)) /
  (max(values(footprint), na.rm=TRUE) - min(values(footprint), na.rm=TRUE))

distortion_magnitude <- ifel(distortion_type == 1,
                             round(2 + 2 * footprint_norm),
                             distortion_magnitude) # subsidy; at least distortion of 2 for subsidies and scaling factor of 2 as well
distortion_magnitude <- ifel(distortion_type == 2,
                             round(1 + 2 * footprint_norm),
                             distortion_magnitude) # tax
distortion_magnitude <- ifel(distortion_type == 3,
                             round(1 + 3 * footprint_norm),
                             distortion_magnitude) # regulation
distortion_magnitude <- ifel(distortion_type == 4,
                             round(3 + 3 * footprint_norm),
                             distortion_magnitude) # externality
distortion_magnitude <- ifel(distortion_type == 5,
                             round(4 + 3 * footprint_norm),
                             distortion_magnitude) # monopoly
distortion_magnitude <- ifel(distortion_type == 6,
                             round(2 + 3 * footprint_norm),
                             distortion_magnitude) # information
distortion_magnitude <- ifel(distortion_type == 7,
                             round(2 + 4 * footprint_norm),
                             distortion_magnitude) # public good; very high scaling as undersupply worsens dramatically with intensity

names(distortion_magnitude) <- "distortion_magnitude"
levels(distortion_magnitude) <- data.frame(id = 0:7,
                                           magnitude = c("none", "minimal", "low",
                                                         "moderate", "substantial",
                                                         "high", "severe",
                                                         "extreme"))

# plot
plot(c(footprint, distortion_type, distortion_magnitude))


# 5. calculate GDP adjustment estimates ----
# estimate hidden ecological costs as % of conventional GDP
exceedance_finite <- ifel(is.finite(exceedance), exceedance, NA)
exceedance_max <- global(exceedance_finite, "max", na.rm = TRUE)[1, 1]
exceedance_norm <- clamp(exceedance_finite / exceedance_max, 0, 1)

# water and soil limitations typically have higher remediation costs
resource_weights <- classify(deficit_type,
                             rcl = matrix(c(0, 0.5,    # no limitation: moderate impact
                                            1, 1.3,    # water limitation: high cost
                                            2, 0.8,    # forage limitation: moderate cost
                                            3, 1.2),   # soil limitation: high cost
                                          ncol=2, byrow=TRUE))

# resource adjustment component
resource_adjustment <- resource_weights * (deficit_magnitude / 7)

# market distortion component
distortion_weights <- classify(distortion_type,
                               rcl = matrix(c(0, 0.2,    # no distortion: minimal impact
                                              1, 1.1,    # subsidy: often hides true costs
                                              2, 0.8,    # tax: may already internalize costs
                                              3, 0.9,    # regulation: partial internalization
                                              4, 1.3,    # externality: significant hidden costs
                                              5, 1.2,    # monopoly: market inefficiency costs
                                              6, 1.0,    # information failure: moderate impact
                                              7, 0.7),   # public good: often undervalued
                                            ncol=2, byrow=TRUE))

# distortion effect
distortion_effect <- distortion_weights * (distortion_magnitude / 7)

# calculate combined GDP adjustment percentage with weighted components
gdp_adjustment <-
  (0.4 * exceedance_norm) +           # highest weight: carrying capacity exceedance
  (0.3 * resource_adjustment) +       # medium weight: resource limitations
  (0.3 * distortion_effect)           # medium weight: market distortions

# ecological systems often show non-linear impacts beyond certain thresholds
gdp_adjustment <- ifel(exceedance > 1.2,
                       gdp_adjustment * 1.5,  # 50% amplification when >20% over capacity
                       gdp_adjustment)

# Constrain values between 0-1 (0-100%)
gdp_adjustment <- clamp(gdp_adjustment, 0, 1)

# Scale to 0-63 range for 6-bit encoding
gdp_adjustment <- round(gdp_adjustment * 63)
names(gdp_adjustment) <- "gdp_adjustment"


# 6. calculate economic-ecological misalignment index ----
# component 1: ecological pressure - how much current use exceeds sustainable capacity
ecological_pressure <- clamp(exceedance / exceedance_max, 0, 1)

# component 2: resource constraint severity - limitations in critical resources
resource_constraint <- deficit_magnitude / 7

# component 3: market failure - economic signals that fail to reflect ecological realities
market_failure <- distortion_magnitude / 7

# component 4: hidden costs - ecological costs not captured in economic accounting
hidden_costs <- gdp_adjustment / 63

# base misalignment components
base_misalignment <-
  (0.35 * ecological_pressure) +    # highest weight: direct ecological overshoot
  (0.25 * resource_constraint) +    # medium weight: resource limitations
  (0.25 * market_failure) +         # medium weight: market distortions
  (0.15 * hidden_costs)             # lower weight: derived metric

# calculate neighborhood effect to capture spatial interdependencies that
# reflect how sustainability issues cross property boundaries and require
# landscape-level solutions
neighborhood_effect <- focal(base_misalignment,
                             w = matrix(1, 3, 3),  # 3x3 moving window
                             fun = "mean",
                             na.rm = TRUE)

# normalize neighborhood effect
ne_min <- global(neighborhood_effect, "min", na.rm = TRUE)[1, 1]
ne_max <- global(neighborhood_effect, "max", na.rm = TRUE)[1, 1]
neighborhood_norm <- (neighborhood_effect - ne_min) / (ne_max - ne_min)

# calculate final misalignment with spatial effects
misalignment <- (0.7 * base_misalignment) + (0.3 * neighborhood_norm)

# constrain to 0-1 range
misalignment <- clamp(misalignment, 0, 1)

# scale to 0-255 range for 8-bit encoding
misalignment <- round(misalignment * 255)
names(misalignment) <- "misalignment"


# 7. generate system trajectory classifications ----
# this adds a critical temporal dimension to the analysis

# normalize inputs for clarity
misalign_norm <- misalignment/255
deficit_norm <- deficit_magnitude/7

# generate individual condition masks to improve readability
is_sustainable <- exceedance < 1
is_moderate_exceedance <- exceedance >= 1 & exceedance < 1.5
is_high_exceedance <- exceedance >= 1.5
is_approaching_threshold <- exceedance >= 0.9 & exceedance < 1

is_low_misalignment <- misalign_norm < 0.3
is_moderate_misalignment <- misalign_norm >= 0.3 & misalign_norm < 0.6
is_high_misalignment <- misalign_norm >= 0.6

is_resource_limited <- deficit_norm > 0.5

system_trajectory <- carrying_capacity
values(system_trajectory) <- NA

# Classifications are applied in order of increasing severity. Later (worse)
# conditions deliberately overwrite earlier ones where ranges overlap, so that
# a cell matching both e.g. "sustainable intensification" and "cyclic
# fluctuation" is assigned the more severe trajectory.

# stable systems: low misalignment, well below carrying capacity
system_trajectory <- ifel(misalign_norm < 0.3 & exceedance < 0.7, 0, system_trajectory)

# sustainable intensification: low misalignment, approaching but not exceeding capacity
system_trajectory <- ifel(misalign_norm < 0.3 & exceedance >= 0.7 & exceedance < 1, 1, system_trajectory)

# approaching capacity: moderate misalignment, close to capacity
system_trajectory <- ifel(misalign_norm >= 0.3 & misalign_norm < 0.4 & exceedance >= 0.8 & exceedance < 1, 2, system_trajectory)

# cyclic fluctuation: significant resource limitations driving oscillations
system_trajectory <- ifel(deficit_magnitude/7 > 0.6 & exceedance >= 0.7 & exceedance <= 1.3, 3, system_trajectory)

# gradual degradation: moderate misalignment, moderately exceeding capacity
system_trajectory <- ifel(misalign_norm >= 0.3 & misalign_norm < 0.6 & exceedance >= 1 & exceedance < 1.5, 4, system_trajectory)

# rapid degradation: high misalignment, significantly exceeding capacity
system_trajectory <- ifel(misalign_norm >= 0.6 & exceedance >= 1.5, 5, system_trajectory)

# approaching threshold: moderate to high misalignment, near critical threshold
system_trajectory <- ifel(misalign_norm >= 0.4 & exceedance >= 0.9 & exceedance < 1.1, 6, system_trajectory)

# post-threshold reorganization: very high misalignment but reduced pressure
system_trajectory <- ifel(misalign_norm > 0.7 & exceedance < 0.9, 7, system_trajectory)

# default to gradual degradation if exceeding capacity
system_trajectory <- ifel(is.na(system_trajectory) & exceedance >= 1, 4, system_trajectory)

# default to stable for any remaining areas
system_trajectory <- ifel(is.na(system_trajectory), 0, system_trajectory)

names(system_trajectory) <- "system_trajectory"
levels(system_trajectory) <- data.frame(
  id = 0:7,
  trajectory = c("stable", "sustainable intensification", "approaching capacity",
    "cyclic fluctuation", "gradual degradation", "rapid degradation",
    "approaching threshold", "post-threshold reorganization")
)

# plot
plot(c(gdp_adjustment, misalignment, system_trajectory))


# 8. calculate intervention priorities ----
#
# different trajectories have different implications for intervention urgency
trajectory_weights <- classify(system_trajectory,
                               rcl = matrix(c(
                                 0, 0.0,  # stable: lowest priority
                                 1, 0.2,  # sustainable intensification: low priority
                                 2, 0.3,  # approaching capacity: moderate priority
                                 3, 0.5,  # cyclic fluctuation: moderate-high priority
                                 4, 0.7,  # gradual degradation: high priority
                                 5, 0.9,  # rapid degradation: very high priority
                                 6, 1.0,  # approaching threshold: highest priority
                                 7, 0.8   # post-threshold reorganization: high priority (window of opportunity)
                               ), ncol=2, byrow=TRUE))

# higher values indicate areas where intervention is more critical
intervention_priority <-
  (0.35 * (misalignment/255)) +                      # economic-ecological misalignment (35%)
  (0.30 * clamp(exceedance / exceedance_max, 0, 1)) + # carrying capacity exceedance (30%), capped at 1
  (0.15 * (deficit_magnitude/7)) +                   # resource limitation severity (15%)
  (0.20 * trajectory_weights)                        # system trajectory (increased to 20%)

# Scale to 0-31 range for 5-bit encoding
intervention_priority <- round(intervention_priority * 31)
names(intervention_priority) <- "intervention_priority"


# 9. calculate cross-sectoral synergy potential ----
# identifies opportunities where food security interventions
# simultaneously advance other development goals
water_synergy <- deficit_type == 1
ecosystem_synergy <- exceedance > 1.5 & system_trajectory >= 4
economic_synergy <- misalignment/255 > 0.6
climate_synergy <- deficit_magnitude/7 > 0.5

# combine into a single bitmap (0-15)
cross_sectoral_synergy <-
  (water_synergy * 8) +          # bit 3: water security
  (ecosystem_synergy * 4) +      # bit 2: ecosystem conservation
  (economic_synergy * 2) +       # bit 1: economic development
  (climate_synergy * 1)          # bit 0: climate resilience
names(cross_sectoral_synergy) <- "synergies"
levels(cross_sectoral_synergy) <- data.frame(
  id = 0:15,
  synergies = c("none", "climate", "economic", "economic + climate", "ecosystem",
                 "ecosystem + climate", "ecosystem + economic",
                 "ecosystem + economic + climate", "water", "water + climate",
                 "water + economic", "water + economic + climate",
                 "water + ecosystem", "water + ecosystem + climate",
                 "water + ecosystem + economic", "all sectors")
)

# plot ----
plot(c(intervention_priority, cross_sectoral_synergy))


# 10. Create bitfield registry ----
#
socio_registry <- bf_registry(
  name = "socioeconomic_planning",
  description = "Socioeconomic analysis and intervention planning metrics",
  template = ecoBitfield)

socio_registry <- bf_map(
  protocol = "category",
  data = distortion_type,
  registry = socio_registry,
  name = "Market distortion type",
  x = distortion,
  na.val = 0)

socio_registry <- bf_map(
  protocol = "category",
  data = distortion_magnitude,
  registry = socio_registry,
  name = "Market distortion magnitude",
  x = magnitude,
  na.val = 0)

socio_registry <- bf_map(
  protocol = "integer",
  data = gdp_adjustment,
  registry = socio_registry,
  name = "GDP adjustment estimates",
  x = gdp_adjustment,
  fields = list(significand = 6),
  na.val = 0) # bf_analyze(gdp_adjustment) shows that range doesn't include 0, so setting na.val to this to reduce bit requirements by 1.

socio_registry <- bf_map(
  protocol = "integer",
  data = misalignment,
  registry = socio_registry,
  name = "Economic-Ecological Misalignment Index",
  x = misalignment,
  fields = list(significand = 8),
  na.val = 0)

socio_registry <- bf_map(
  protocol = "category",
  data = system_trajectory,
  registry = socio_registry,
  name = "System trajectory",
  x = trajectory,
  na.val = 0)

socio_registry <- bf_map(
  protocol = "integer",
  data = intervention_priority,
  registry = socio_registry,
  name = "Intervention priority",
  x = intervention_priority,
  fields = list(significand = 5),
  na.val = 0)

socio_registry <- bf_map(
  protocol = "category",
  data = cross_sectoral_synergy,
  registry = socio_registry,
  name = "Cross-sectoral synergy",
  x = synergies,
  na.val = 0)

socio_field <- bf_encode(registry = socio_registry)

# 11. export items ----
#
writeRaster(x = socio_field, filename = "intervBitfield.tif", datatype = "INT4U", overwrite = TRUE)

bf_export(registry = socio_registry, format = "yaml", file = paste0(getwd(), "/meta/intervMeta.yml"))
saveRDS(object = socio_registry, file = "intervReg.rds")