project_2.R

# =============================================================================
# Project 2: Ecological Economics & Carrying Capacity Assessment
# =============================================================================
#
# This script demonstrates how a bitfield from an upstream analysis (project 1)
# can be decoded and combined with new data to produce derived ecological
# metrics. It is a SIMULATION designed to showcase the bitfield approach and
# produce plausible visuals — not a validated ecological model. Climate, soil,
# and topography data are real; carrying capacity estimates and resource
# limitation heuristics are simplified approximations to illustrate
# cross-workflow metadata preservation in a 24-bit encoding.
#
# Workflow:
#   1. Load geophysical basis data (climate, soil, topography)
#   2. Decode the livestock bitfield from project 1
#   3. Estimate ecological carrying capacity (Miami model NPP)
#   4. Determine carrying capacity exceedance
#   5. Calculate resource limitations (water, forage, soil)
#   6. Encode everything into a 24-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 ----
#
## Download and crop climate, soil and elevation data ----
tavg <- worldclim_global(var = "tavg", res = 0.5, path = tempdir()) |>
  crop(area_extent)
tmax <- worldclim_global(var = "tmax", res = 0.5, path = tempdir()) |>
  crop(area_extent)
prec <- worldclim_global(var = "prec", res = 0.5, path = tempdir()) |>
  crop(area_extent)
soil_pH <- soil_world(var = "phh2o", depth = 5, path = tempdir()) |>
  crop(area_extent)
soil_carbon <- soil_world(var = "soc", depth = 5, path = tempdir()) |>
  crop(area_extent)
soil_clay <- soil_world(var = "clay", depth = 15, path = tempdir()) |> # for some reason the depth = 5 layer doesn't work currently
  crop(area_extent)

elevation <- elevation_30s(country = "DEU", mask = FALSE, path = tempdir()) |>
  crop(area_extent)
names(elevation) <- "elevation"
slope <- terrain(elevation, "slope", unit = "degrees")
aspect <- terrain(elevation, "aspect", unit = "degrees")

# Summarize climate variables
mean_temp <- mean(tavg)
names(mean_temp) <- "yearly_average_temperature"
max_month <- max(tmax)
names(max_month) <- "warmest_month"
annual_precip <- sum(prec)
names(annual_precip) <- "yearly_total_precipitation"

# Plot
plot(c(mean_temp, annual_precip, soil_pH, soil_carbon, soil_clay))


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

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

animals_median <- decoded$`Median livestock density`
animals_sd <- decoded$`Standard deviation`
skewness <- decoded$Skewness
kurtosis <- decoded$Kurtosis
dist_values <- decoded$`Distribution type`

animals_mean <- rast("animals_mean.tif")

# Load landcover data from project_1 (needed for carrying capacity calculation)
# These were used in project_1 but not encoded in the bitfield
landcover <- list()
lc_meta <- data.frame(id = c(1, 2, 4:9),
                     var = c("trees", "grassland", "shrubs", "cropland", "built",
                             "bare", "water", "wetland"),
                     name = c("Forest", "Grassland", "Shrubland", "Cropland",
                              "Urban", "Bare", "Water", "Wetland"))

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$shrubs <- landcover(var = "shrubs", path = tempdir()) |>
  crop(area_extent) |> unlist()
landcover$bare <- landcover(var = "bare", path = tempdir()) |>
  crop(area_extent) |> unlist()


# 3. calculate carrying capacity ----
#
# calculate NPP using Miami model with climate data
npp_temp <- 3000 / (1 + exp(1.315 - 0.119 * mean_temp))
npp_precip <- 3000 * (1 - exp(-0.000664 * annual_precip))
potential_npp <- min(npp_temp, npp_precip)
names(potential_npp) <- "potential_npp"

# determine landcover factor from simplified classifications
land_factor <-
  1.5 * landcover$grassland +                       # grassland most productive for livestock
  1.0 * landcover$cropland +                        # cropland provides some grazing
  0.7 * landcover$trees +                           # forest has limited grazing
  0.5 * landcover$shrubs +                          # shrubland has minimal forage
  0.2 * landcover$bare                              # bare land can support minimal vegetation
names(land_factor) <- "land_support"

# determine soil factor
pH_factor <- 1 - (0.2 * abs(soil_pH - 6.5) / 2)     # optimal pH for grassland is around 6-7, with decreasing productivity away from this range
carbon_factor <- 0.7 + (0.3 * (soil_carbon / 200))  # organic carbon improves soil fertility and moisture retention
clay_factor <- 1 - (0.2 * abs(soil_clay - 25) / 20) # optimal clay is around 20-30%, with decreased productivity at extremes
soil_factor <- (pH_factor + carbon_factor + clay_factor) / 3
names(soil_factor) <- "soil_support"

# determine slope factor (Otieno, 2019)
slope_factor <- 1.0 - (min(slope, 45) / 45)^1.5
names(slope_factor) <- "slope_limitation"

# define sustainable utilization rate
proper_use_factor <- 0.25

# Calculate actual carrying capacity
usable_npp <- potential_npp *
  land_factor *
  soil_factor *
  slope_factor *
  proper_use_factor

# Convert NPP (g/m²/yr) to livestock units per hectare assuming 2300 kg dry matter per livestock unit per year
carrying_capacity <- usable_npp * 10 / 2300
names(carrying_capacity) <- "carrying_capacity"

# plot
plot(c(potential_npp, land_factor, soil_factor, slope_factor, carrying_capacity))


# 4. Determine carrying capacity exceedance ----
exceedance_prob <- .exceedance_probability(animals_mean, animals_sd, dist_values,
                                           skewness, kurtosis, carrying_capacity)

exceedance <- animals_mean / carrying_capacity
names(exceedance) <- "exceedance_magnitude"

# Plot
plot(c(animals_mean, carrying_capacity, exceedance_prob, exceedance))


# 5. Calculate resource limitations ----
## Water limitation ----
water_demand <- animals_mean * 45                   # water need per livestock unit (liters/day)

# south-facing slopes have increased evaporation, reducing available water
south_factor <- (abs(aspect - 180) / 180) * 0.3 + 0.7  # 0.7-1.0 range (0.7 = south-facing, 1.0 = north-facing)

# get the driest month from the precipitation data, reduced by 30% for runoff
driest_month <- min(prec) * 0.7
driest_daily <- driest_month / 30                    # convert to daily precipitation

# increased animal water requirements during heat periods
heat_stress_factor <- 1 + ((max_month - 15) / 4)

# apply to demand
water_demand <- water_demand * heat_stress_factor

# calculate water availability - baseline calculation
available_water <- driest_daily * 0.4 * 10000 *      # precipitation with retention factor
  0.05                                                # ~5% becomes available for drinking

# areas with cropland have reduced water availability due to competition
crop_intensity <- (landcover$cropland)^1.5           # non-linear relationship with intensity
crop_extraction <- 1 - (crop_intensity * 0.3)        # up to 30% reduction in high-intensity cropland

# apply all modifying factors to animal-available water
adjusted_water <- available_water *
  (1 - 0.8 * (slope / 45)) *                         # topographic slope effect
  south_factor *                                      # south-facing slope effect
  crop_extraction                                     # agricultural water extraction

# calculate deficit ratio using heat-adjusted demand
water_deficit <- water_demand / adjusted_water
water_deficit[water_deficit < 1] <- 0                # no deficit where supply exceeds demand
names(water_deficit) <- "water_deficit"


## Forage limitation ----
forage_demand <- animals_mean * 10                  # kg dry matter per livestock unit per day

forage_availability <- potential_npp * 0.0005 *     # kg/m²/yr
  (0.7 * landcover$grassland +
   0.5 * landcover$cropland +
   0.2 * landcover$trees)
forage_availability <- forage_availability * 10000  # annual production per ha
forage_deficit <- (forage_demand * 365) / (forage_availability + 0.001)  # add small constant to avoid division by zero
forage_deficit[forage_deficit < 1] <- 0
names(forage_deficit) <- "forage_deficit"

## Soil fertility limitation ----
soil_ph_factor <- 1 - abs(soil_pH - 6.5) / 3.5      # optimal pH around 6.5
soil_carbon_factor <- soil_carbon / 200             # higher carbon is better for soil resilience
soil_clay_factor <- 1 - abs(soil_clay - 25) / 75    # optimal clay around 25%

soil_fertility <- min(c(soil_ph_factor, soil_carbon_factor, soil_clay_factor))
soil_fertility <- ifel(soil_fertility < 0, 0, soil_fertility)  # ensure non-negative

base_impact <- 0.2
slope_multiplier <- 1 + (slope / 45)^1.5            # steeper slope means more erosion risk
soil_demand <- animals_mean * base_impact * slope_multiplier

soil_deficit <- soil_demand / (soil_fertility + 0.001)  # add small constant to avoid division by zero
soil_deficit[soil_deficit < 1] <- 0
names(soil_deficit) <- "soil_deficit"

## Resource limitation type ----
# 0=no limitation, 1=water, 2=forage, 3=soil
deficit_type <- ifel(water_deficit > forage_deficit & water_deficit > soil_deficit & water_deficit > 0, 1,
                      ifel(forage_deficit > soil_deficit & forage_deficit > 0, 2,
                           ifel(soil_deficit > 0, 3, 0)))
names(deficit_type) <- "deficit_type"
levels(deficit_type) <- data.frame(id = 0:3,
                                    limitation = c("none", "water", "forage", "soil"))

## Resource limitation magnitude ----
deficit_magnitude <- ifel(deficit_type == 1, log10(water_deficit + 1),
                           ifel(deficit_type == 2, log10(forage_deficit + 1),
                                ifel(deficit_type == 3, log10(soil_deficit + 1), 0)))

# Scale to 0-7 range for 3-bit encoding
deficit_magnitude <- (deficit_magnitude / max(values(deficit_magnitude), na.rm = TRUE)) * 7
names(deficit_magnitude) <- "deficit_magnitude"

# Plot
plot(c(water_deficit, forage_deficit, soil_deficit, deficit_type, deficit_magnitude))


# 6. Create bitfield according to Table A2 specifications ----
#
# Create registry for the ecological analysis
eco_registry <- bf_registry(
  name = "ecological_overshoot_analysis",
  description = "Ecological carrying capacity analysis with 24-bit encoding",
  template = lstBitfield)

eco_registry <- bf_map(
  protocol = "numeric",
  data = carrying_capacity,
  registry = eco_registry,
  name = "Ecological carrying capacity",
  x = carrying_capacity,
  fields = list(exponent = 3, significand = 4))

eco_registry <- bf_map(
  protocol = "numeric",
  data = exceedance_prob,
  registry = eco_registry,
  name = "Exceedance risk",
  x = exceedance_probability,
  fields = list(exponent = 4, significand = 2))

eco_registry <- bf_map(
  protocol = "numeric",
  data = exceedance,
  registry = eco_registry,
  name = "Exceedance magnitude",
  x = exceedance_magnitude,
  fields = list(exponent = 3, significand = 3))

eco_registry <- bf_map(
  protocol = "category",
  data = deficit_type,
  registry = eco_registry,
  name = "Resource limitation type",
  x = limitation,
  na.val = 0)

eco_registry <- bf_map(
  protocol = "integer",
  data = deficit_magnitude,
  registry = eco_registry,
  name = "Resource limitation magnitude",
  x = deficit_magnitude,
  range = c(0, 7),
  fields = list(significand = 3),
  na.val= 0)

# Encode the bitfield
eco_rast <- bf_encode(registry = eco_registry)

# 7. export items ----
#
writeRaster(x = eco_rast, filename = "ecoBitfield.tif", datatype = "INT4U", overwrite = TRUE)

bf_export(registry = eco_registry, format = "yaml", file = paste0(getwd(), "/meta/ecoMeta.yml"))
saveRDS(object = eco_registry, file = "eco_reg.rds")