## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  fig.alt = "Figura generada por la viñeta; ver texto para detalles.",
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  warning = FALSE,
  message = FALSE
)

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

## -----------------------------------------------------------------------------
# Define comprehensive farm-gate boundaries
boundaries <- set_system_boundaries("farm_gate")
print(boundaries)

# Alternative: cradle-to-farm-gate (includes upstream emissions)
boundaries_extended <- set_system_boundaries("cradle_to_farm_gate")
print(boundaries_extended)

# Custom boundaries example
boundaries_partial <- set_system_boundaries(
  scope = "partial",
  include = c("enteric", "manure", "soil")
)
print(boundaries_partial)

## -----------------------------------------------------------------------------
# Detailed herd data
herd_data <- list(
  # Main herd
  dairy_cows_milking = 120,
  dairy_cows_dry = 30,
  
  # Young stock
  heifers_total = 45,
  calves_total = 35,
  bulls_total = 3,
  
  # Animal characteristics (kg live weight)
  body_weight_cows = 580,
  body_weight_heifers = 380,
  body_weight_calves = 180,
  body_weight_bulls = 750,
  
  # Production parameters
  milk_yield_per_cow = 6300,  # kg/cow/year
  annual_milk_litres = 950000,
  fat_percent = 3.7,
  protein_percent = 3.3,
  milk_density = 1.032
)

print(herd_data)

## -----------------------------------------------------------------------------
# Feed inputs (all in kg dry matter per year)
feed_data <- list(
  # Purchased feeds
  concentrate_kg = 220000,
  grain_dry_kg = 80000,
  grain_wet_kg = 45000,
  ration_kg = 60000,
  byproducts_kg = 25000,
  proteins_kg = 35000,
  
  # Nutritional parameters
  dry_matter_intake_cows = 19.5,  # kg DM/cow/day
  dry_matter_intake_heifers = 11.0,
  dry_matter_intake_calves = 6.0,
  dry_matter_intake_bulls = 14.0,
  ym_percent = 6.3  # Methane conversion factor
)

print(feed_data)

## -----------------------------------------------------------------------------
# Detailed land use breakdown
land_data <- list(
  # Total areas (hectares)
  area_total = 200,
  area_productive = 185,
  area_fertilized = 160,
  
  # Area breakdown
  pasture_permanent = 140,
  pasture_temporary = 30,
  crops_feed = 12,
  crops_cash = 3,
  infrastructure = 8,
  woodland = 7,
  
  # Soil and climate
  soil_type = "well_drained",
  climate_zone = "temperate",
  
  # Nitrogen inputs (kg N/year)
  n_fertilizer_synthetic = 2400,
  n_fertilizer_organic = 500,
  n_excreta_pasture = 15000,  # Estimated from grazing
  n_crop_residues = 800
)

print(land_data)

## -----------------------------------------------------------------------------
# Energy use breakdown
energy_data <- list(
  # Fuel consumption (litres/year)
  diesel_litres = 12000,
  petrol_litres = 1800,
  
  # Other energy sources
  lpg_kg = 600,
  natural_gas_m3 = 0,
  electricity_kwh = 48000,
  
  # Country for electricity factors
  country = "UY"
)

print(energy_data)

## -----------------------------------------------------------------------------
# Other purchased inputs
other_inputs <- list(
  # Materials (kg/year)
  plastic_kg = 450,
  
  # Transport (optional)
  transport_km = 120,  # Average transport distance for feeds
  
  # Fertilizer types
  fert_type = "mixed",
  plastic_type = "mixed",
  
  # Regional factors
  region = "global"  # Can be "EU", "US", "Brazil", etc.
)

print(other_inputs)

## -----------------------------------------------------------------------------
# Calculate enteric emissions for each animal category
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
)

# Summary of enteric emissions
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 emissions
total_enteric <- enteric_cows$co2eq_kg + enteric_heifers$co2eq_kg + 
                enteric_calves$co2eq_kg + enteric_bulls$co2eq_kg

## -----------------------------------------------------------------------------
# Calculate manure emissions for the entire herd
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",  # Extensive grazing system
  tier = 2,
  avg_body_weight = 500,  # Weighted average
  diet_digestibility = 0.67,
  climate = "temperate",
  include_indirect = TRUE,
  boundaries = boundaries
)

print(manure_emissions)

## -----------------------------------------------------------------------------
# Calculate soil emissions from all N sources
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
)

print(soil_emissions)

## -----------------------------------------------------------------------------
# Calculate emissions from energy use
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,  # Only direct emissions
  boundaries = boundaries
)

print(energy_emissions)

## -----------------------------------------------------------------------------
# Calculate emissions from purchased inputs
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
)

print(input_emissions)

## -----------------------------------------------------------------------------
# Create combined enteric emissions object
enteric_combined <- list(
  source = "enteric",
  co2eq_kg = total_enteric
)

# Calculate total emissions
total_emissions <- calc_total_emissions(
  enteric_combined,
  manure_emissions,
  soil_emissions,
  energy_emissions,
  input_emissions
)

total_emissions

## ----fig.width=8, fig.height=6------------------------------------------------
# Create detailed breakdown
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)
)

# Create bar chart
ggplot(emission_breakdown, aes(x = reorder(Source, Emissions), y = Emissions)) +
  geom_col(fill = "steelblue", 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))

## -----------------------------------------------------------------------------
# Calculate emissions per kg of FPCM
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
)

print(milk_intensity)

## -----------------------------------------------------------------------------
# Calculate emissions per hectare
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
)

print(area_intensity)

## -----------------------------------------------------------------------------
# Benchmark against Uruguayan standards
area_benchmark <- benchmark_area_intensity(
  cf_area_intensity = area_intensity,
  region = "uruguay"
)

print(area_benchmark$benchmarking)

## -----------------------------------------------------------------------------
# Calculate emissions per animal category
per_animal_analysis <- data.frame(
  Category = c("Dairy Cows", "All Animals"),
  Number = c(150, total_animals),
  Total_Emissions = c(total_emissions$total_co2eq, total_emissions$total_co2eq),
  Emissions_per_Head = c(
    total_emissions$total_co2eq / 150,
    total_emissions$total_co2eq / total_animals
  ),
  Milk_per_Head = c(herd_data$annual_milk_litres / 150, NA)
)

kable(per_animal_analysis, 
      digits = 0,
      caption = "Per-Animal Emission Analysis")

## -----------------------------------------------------------------------------
# Calculate feed-related metrics
total_purchased_feed <- sum(
  feed_data$concentrate_kg,
  feed_data$grain_dry_kg,
  feed_data$grain_wet_kg,
  feed_data$ration_kg,
  feed_data$byproducts_kg,
  feed_data$proteins_kg
)

feed_analysis <- data.frame(
  Metric = c("Total Feed Purchases", "Feed Emissions", "Feed CO2eq per kg DM",
             "Feed Efficiency", "Milk from Feed"),
  Value = c(
    total_purchased_feed,
    input_emissions$total_co2eq_kg,
    input_emissions$total_co2eq_kg / total_purchased_feed,
    herd_data$annual_milk_litres / total_purchased_feed,
    herd_data$annual_milk_litres
  ),
  Unit = c("kg DM", "kg CO₂eq", "kg CO₂eq/kg DM", "L milk/kg DM", "L")
)

kable(feed_analysis, digits = 2, caption = "Feed Efficiency Analysis")

## -----------------------------------------------------------------------------
# Scenario: 10% reduction in concentrate use with maintained production
improved_inputs <- calc_emissions_inputs(
  conc_kg = feed_data$concentrate_kg * 0.9,  # 10% reduction
  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
)

# Calculate total for improved scenario
total_improved <- calc_total_emissions(
  enteric_combined,
  manure_emissions,
  soil_emissions,
  energy_emissions,
  improved_inputs
)

# Compare scenarios
scenario_comparison <- data.frame(
  Scenario = c("Baseline", "Improved Feed Efficiency"),
  Total_Emissions = c(total_emissions$total_co2eq, total_improved$total_co2eq),
  Input_Emissions = c(input_emissions$total_co2eq_kg, improved_inputs$total_co2eq_kg),
  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")

## -----------------------------------------------------------------------------
# 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-----------------------------------------------
# Prepare data for comprehensive visualization
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
  ),
  Category = c("Enteric", "Enteric", "Manure", "Soil", "Energy", "Inputs")
)

# Create stacked bar chart
ggplot(detailed_emissions, aes(x = "Farm Emissions", y = Emissions, fill = Source)) +
  geom_col() +
  geom_text(aes(label = ifelse(Emissions > 2000, 
                              paste0(round(Emissions/1000, 1), "k"), "")),
            position = position_stack(vjust = 0.5),
            color = "white", fontweight = "bold") +
  labs(title = "Estancia Las Flores - Carbon Footprint Breakdown",
       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") +
  scale_fill_brewer(type = "qual", palette = "Set2")

## -----------------------------------------------------------------------------
# Calculate key performance indicators
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)
  ),
  Benchmark = c("< 1.2", "< 8,000", "< 8,500", "> 85%", "> 6,000", "0.5-1.2"),
  Performance = c(
    ifelse(milk_intensity$intensity_co2eq_per_kg_fpcm < 1.2, "Good", "Needs Improvement"),
    ifelse(area_intensity$intensity_per_total_ha < 8000, "Good", "Needs Improvement"),
    ifelse(area_intensity$intensity_per_productive_ha < 8500, "Good", "Needs Improvement"),
    ifelse(area_intensity$land_use_efficiency > 0.85, "Good", "Needs Improvement"),
    ifelse(herd_data$annual_milk_litres / 150 > 6000, "Good", "Needs Improvement"),
    "Within Range"
  )
)

kable(kpi_summary, caption = "Key Performance Indicators")