## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  fig.alt = "Figure generated by this vignette; see the surrounding text for details.",
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  warning = FALSE,
  message = FALSE
)

library(cowfootR)
library(ggplot2)
library(dplyr)
library(knitr)

## -----------------------------------------------------------------------------
# Most common scope: farm-gate boundaries
boundaries <- set_system_boundaries("farm_gate")
boundaries

## -----------------------------------------------------------------------------
# Includes upstream emissions (e.g., purchased inputs)
boundaries_extended <- set_system_boundaries("cradle_to_farm_gate")
boundaries_extended

# Custom selection of sources
boundaries_partial <- set_system_boundaries(
  scope = "partial",
  include = c("enteric", "manure", "soil")
)
boundaries_partial

## -----------------------------------------------------------------------------
herd_data <- list(
  dairy_cows_milking = 120,
  dairy_cows_dry = 30,
  heifers_total = 45,
  calves_total = 35,
  bulls_total = 3,
  body_weight_cows = 580,
  body_weight_heifers = 380,
  body_weight_calves = 180,
  body_weight_bulls = 750,
  milk_yield_per_cow = 6300,
  annual_milk_litres = 950000,
  fat_percent = 3.7,
  protein_percent = 3.3,
  milk_density = 1.032
)

herd_data

## -----------------------------------------------------------------------------
feed_data <- list(
  concentrate_kg = 220000,
  grain_dry_kg = 80000,
  grain_wet_kg = 45000,
  ration_kg = 60000,
  byproducts_kg = 25000,
  proteins_kg = 35000,
  dry_matter_intake_cows = 19.5,
  dry_matter_intake_heifers = 11.0,
  dry_matter_intake_calves = 6.0,
  dry_matter_intake_bulls = 14.0,
  ym_percent = 6.3
)

feed_data

## -----------------------------------------------------------------------------
land_data <- list(
  area_total = 200,
  area_productive = 185,
  area_fertilized = 160,
  pasture_permanent = 140,
  pasture_temporary = 30,
  crops_feed = 12,
  crops_cash = 3,
  infrastructure = 8,
  woodland = 7,
  soil_type = "well_drained",
  climate_zone = "temperate",
  n_fertilizer_synthetic = 2400,
  n_fertilizer_organic = 500,
  n_excreta_pasture = 15000,
  n_crop_residues = 800
)

land_data

## -----------------------------------------------------------------------------
energy_data <- list(
  diesel_litres = 12000,
  petrol_litres = 1800,
  lpg_kg = 600,
  natural_gas_m3 = 0,
  electricity_kwh = 48000,
  country = "UY"
)

energy_data

## -----------------------------------------------------------------------------
other_inputs <- list(
  plastic_kg = 450,
  transport_km = 120,
  fert_type = "mixed",
  plastic_type = "mixed",
  region = "global"
)

other_inputs

## -----------------------------------------------------------------------------
enteric_cows <- calc_emissions_enteric(
  n_animals = herd_data$dairy_cows_milking + herd_data$dairy_cows_dry,
  cattle_category = "dairy_cows",
  avg_milk_yield = herd_data$milk_yield_per_cow,
  avg_body_weight = herd_data$body_weight_cows,
  dry_matter_intake = feed_data$dry_matter_intake_cows,
  ym_percent = feed_data$ym_percent,
  tier = 2,
  boundaries = boundaries
)

enteric_heifers <- calc_emissions_enteric(
  n_animals = herd_data$heifers_total,
  cattle_category = "heifers",
  avg_body_weight = herd_data$body_weight_heifers,
  dry_matter_intake = feed_data$dry_matter_intake_heifers,
  ym_percent = feed_data$ym_percent,
  tier = 2,
  boundaries = boundaries
)

enteric_calves <- calc_emissions_enteric(
  n_animals = herd_data$calves_total,
  cattle_category = "calves",
  avg_body_weight = herd_data$body_weight_calves,
  dry_matter_intake = feed_data$dry_matter_intake_calves,
  tier = 2,
  boundaries = boundaries
)

enteric_bulls <- calc_emissions_enteric(
  n_animals = herd_data$bulls_total,
  cattle_category = "bulls",
  avg_body_weight = herd_data$body_weight_bulls,
  dry_matter_intake = feed_data$dry_matter_intake_bulls,
  tier = 2,
  boundaries = boundaries
)

enteric_summary <- data.frame(
  Category = c("Dairy Cows", "Heifers", "Calves", "Bulls"),
  Animals = c(150, herd_data$heifers_total, herd_data$calves_total, herd_data$bulls_total),
  CH4_kg = c(enteric_cows$ch4_kg, enteric_heifers$ch4_kg, enteric_calves$ch4_kg, enteric_bulls$ch4_kg),
  CO2eq_kg = c(enteric_cows$co2eq_kg, enteric_heifers$co2eq_kg, enteric_calves$co2eq_kg, enteric_bulls$co2eq_kg)
)

kable(enteric_summary, caption = "Enteric emissions by animal category")

total_enteric <- enteric_summary$CO2eq_kg |> sum()

## -----------------------------------------------------------------------------
total_animals <- sum(
  herd_data$dairy_cows_milking, herd_data$dairy_cows_dry,
  herd_data$heifers_total, herd_data$calves_total, herd_data$bulls_total
)

manure_emissions <- calc_emissions_manure(
  n_cows = total_animals,
  manure_system = "pasture",
  tier = 2,
  avg_body_weight = 500,
  diet_digestibility = 0.67,
  climate = "temperate",
  include_indirect = TRUE,
  boundaries = boundaries
)

manure_emissions

## -----------------------------------------------------------------------------
soil_emissions <- calc_emissions_soil(
  n_fertilizer_synthetic = land_data$n_fertilizer_synthetic,
  n_fertilizer_organic = land_data$n_fertilizer_organic,
  n_excreta_pasture = land_data$n_excreta_pasture,
  n_crop_residues = land_data$n_crop_residues,
  area_ha = land_data$area_total,
  soil_type = land_data$soil_type,
  climate = land_data$climate_zone,
  include_indirect = TRUE,
  boundaries = boundaries
)

soil_emissions

## -----------------------------------------------------------------------------
energy_emissions <- calc_emissions_energy(
  diesel_l = energy_data$diesel_litres,
  petrol_l = energy_data$petrol_litres,
  lpg_kg = energy_data$lpg_kg,
  natural_gas_m3 = energy_data$natural_gas_m3,
  electricity_kwh = energy_data$electricity_kwh,
  country = energy_data$country,
  include_upstream = FALSE,
  boundaries = boundaries
)

energy_emissions

## -----------------------------------------------------------------------------
input_emissions <- calc_emissions_inputs(
  conc_kg = feed_data$concentrate_kg,
  fert_n_kg = land_data$n_fertilizer_synthetic,
  plastic_kg = other_inputs$plastic_kg,
  feed_grain_dry_kg = feed_data$grain_dry_kg,
  feed_grain_wet_kg = feed_data$grain_wet_kg,
  feed_ration_kg = feed_data$ration_kg,
  feed_byproducts_kg = feed_data$byproducts_kg,
  feed_proteins_kg = feed_data$proteins_kg,
  region = other_inputs$region,
  fert_type = other_inputs$fert_type,
  plastic_type = other_inputs$plastic_type,
  transport_km = other_inputs$transport_km,
  boundaries = boundaries
)

input_emissions

## -----------------------------------------------------------------------------
enteric_combined <- list(source = "enteric", co2eq_kg = total_enteric)

total_emissions <- calc_total_emissions(
  enteric_combined,
  manure_emissions,
  soil_emissions,
  energy_emissions,
  input_emissions
)

total_emissions

## ----fig.width=8, fig.height=6, fig.alt="Horizontal bar chart showing total farm emissions by source (enteric, manure, soil, energy, inputs) and each source's percentage of total CO2eq."----
emission_breakdown <- data.frame(
  Source = names(total_emissions$breakdown),
  Emissions = as.numeric(total_emissions$breakdown),
  Percentage = round(as.numeric(total_emissions$breakdown) / total_emissions$total_co2eq * 100, 1)
)

ggplot(emission_breakdown, aes(x = reorder(Source, Emissions), y = Emissions)) +
  geom_col(alpha = 0.8) +
  geom_text(aes(label = paste0(Percentage, "%")), hjust = -0.1, size = 3) +
  coord_flip() +
  labs(
    title = "Farm emissions by source",
    subtitle = paste("Total:", format(round(total_emissions$total_co2eq), big.mark = ","), "kg CO₂eq/year"),
    x = "Emission source",
    y = "Emissions (kg CO₂eq/year)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  )

## -----------------------------------------------------------------------------
milk_intensity <- calc_intensity_litre(
  total_emissions = total_emissions,
  milk_litres = herd_data$annual_milk_litres,
  fat = herd_data$fat_percent,
  protein = herd_data$protein_percent,
  milk_density = herd_data$milk_density
)

milk_intensity

## -----------------------------------------------------------------------------
area_breakdown <- list(
  pasture_permanent = land_data$pasture_permanent,
  pasture_temporary = land_data$pasture_temporary,
  crops_feed = land_data$crops_feed,
  crops_cash = land_data$crops_cash,
  infrastructure = land_data$infrastructure,
  woodland = land_data$woodland
)

area_intensity <- calc_intensity_area(
  total_emissions = total_emissions,
  area_total_ha = land_data$area_total,
  area_productive_ha = land_data$area_productive,
  area_breakdown = area_breakdown,
  validate_area_sum = TRUE
)

area_intensity

## -----------------------------------------------------------------------------
area_benchmark <- benchmark_area_intensity(
  cf_area_intensity = area_intensity,
  region = "uruguay"
)

area_benchmark$benchmarking

## -----------------------------------------------------------------------------
improved_inputs <- calc_emissions_inputs(
  conc_kg = feed_data$concentrate_kg * 0.9,
  fert_n_kg = land_data$n_fertilizer_synthetic,
  plastic_kg = other_inputs$plastic_kg,
  feed_grain_dry_kg = feed_data$grain_dry_kg,
  feed_grain_wet_kg = feed_data$grain_wet_kg,
  feed_ration_kg = feed_data$ration_kg,
  feed_byproducts_kg = feed_data$byproducts_kg,
  feed_proteins_kg = feed_data$proteins_kg,
  region = other_inputs$region,
  fert_type = other_inputs$fert_type,
  plastic_type = other_inputs$plastic_type,
  transport_km = other_inputs$transport_km,
  boundaries = boundaries
)

total_improved <- calc_total_emissions(
  enteric_combined,
  manure_emissions,
  soil_emissions,
  energy_emissions,
  improved_inputs
)

scenario_comparison <- data.frame(
  Scenario = c("Baseline", "Improved Feed Efficiency"),
  Total_Emissions = c(total_emissions$total_co2eq, total_improved$total_co2eq),
  Reduction_kg = c(0, total_emissions$total_co2eq - total_improved$total_co2eq),
  Reduction_percent = c(
    0,
    round((total_emissions$total_co2eq - total_improved$total_co2eq) / total_emissions$total_co2eq * 100, 1)
  )
)

kable(scenario_comparison, caption = "Mitigation scenario analysis (illustrative)")

## -----------------------------------------------------------------------------
# Scenario: Switch from pasture to anaerobic digester
improved_manure <- calc_emissions_manure(
  n_cows = total_animals,
  manure_system = "anaerobic_digester",
  tier = 2,
  avg_body_weight = 500,
  diet_digestibility = 0.67,
  climate = "temperate",
  retention_days = 45,
  system_temperature = 35,
  include_indirect = TRUE,
  boundaries = boundaries
)

# Calculate total with improved manure management
total_improved_manure <- calc_total_emissions(
  enteric_combined,
  improved_manure,
  soil_emissions,
  energy_emissions,
  input_emissions
)

manure_comparison <- data.frame(
  System = c("Pasture", "Anaerobic Digester"),
  CH4_kg = c(manure_emissions$ch4_kg, improved_manure$ch4_kg),
  N2O_kg = c(manure_emissions$n2o_total_kg, improved_manure$n2o_total_kg),
  Total_CO2eq = c(manure_emissions$co2eq_kg, improved_manure$co2eq_kg),
  Reduction_kg = c(0, manure_emissions$co2eq_kg - improved_manure$co2eq_kg)
)

kable(manure_comparison, caption = "Manure Management Comparison")

## ----fig.width=10, fig.height=6, fig.alt="Stacked bar chart with a detailed breakdown of farm emissions, splitting enteric emissions into cows vs young stock, plus manure, soil, energy, and purchased inputs."----
detailed_emissions <- data.frame(
  Source = c(
    "Enteric - Cows", "Enteric - Young Stock", "Manure Management",
    "Soil N2O", "Energy Use", "Purchased Inputs"
  ),
  Emissions = c(
    enteric_cows$co2eq_kg,
    enteric_heifers$co2eq_kg + enteric_calves$co2eq_kg + enteric_bulls$co2eq_kg,
    manure_emissions$co2eq_kg,
    soil_emissions$co2eq_kg,
    energy_emissions$co2eq_kg,
    input_emissions$total_co2eq_kg
  )
)

ggplot(detailed_emissions, aes(x = "Farm emissions", y = Emissions, fill = Source)) +
  geom_col() +
  labs(
    title = "Single farm carbon footprint breakdown (detailed view)",
    subtitle = paste("Total:", format(round(total_emissions$total_co2eq), big.mark = ","), "kg CO₂eq/year"),
    x = "",
    y = "Emissions (kg CO₂eq/year)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 14, hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5),
    axis.text.x = element_blank(),
    legend.position = "right"
  )

## -----------------------------------------------------------------------------
kpi_summary <- data.frame(
  Metric = c(
    "Milk intensity (kg CO₂eq/kg FPCM)",
    "Area intensity - total (kg CO₂eq/ha)",
    "Area intensity - productive (kg CO₂eq/ha)",
    "Land use efficiency (%)",
    "Milk yield (L/cow/year)",
    "Stocking rate (cows/ha)"
  ),
  Value = c(
    round(milk_intensity$intensity_co2eq_per_kg_fpcm, 3),
    round(area_intensity$intensity_per_total_ha, 0),
    round(area_intensity$intensity_per_productive_ha, 0),
    round(area_intensity$land_use_efficiency * 100, 1),
    round(herd_data$annual_milk_litres / 150, 0),
    round(150 / land_data$area_total, 2)
  )
)

kable(kpi_summary, caption = "Key performance indicators (illustrative)")