---
title: "Single_Farm_Analysis_cowfootR"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Single_Farm_Analysis_cowfootR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


`
```{r 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)
```

# Comprehensive Single Farm Carbon Footprint Analysis

This vignette demonstrates how to conduct a detailed carbon footprint assessment for an individual dairy farm using all available functions in cowfootR. We'll work through a realistic example with comprehensive data collection and analysis.

## Farm Profile: Estancia Las Flores

For this analysis, we'll assess a medium-sized dairy farm in Uruguay with the following characteristics:

- **Location**: Temperate climate, well-drained soils
- **System**: Mixed grazing with supplementation
- **Herd**: 150 dairy cows plus young stock
- **Area**: 200 hectares total
- **Production**: 950,000 litres annually

## Step 1: System Boundaries Definition

First, we define what emission sources to include in our assessment:

```{r}
# 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)
```

For this analysis, we'll use farm-gate boundaries as they represent the most common assessment scope.

## Step 2: Detailed Farm Data Collection

### Herd Composition and Management

```{r}
# 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 and Nutrition Management

```{r}
# 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)
```

### Land Use and Soil Management

```{r}
# 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 Consumption

```{r}
# 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)
```

### Additional Inputs

```{r}
# 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)
```

## Step 3: Emission Calculations by Source

Now we calculate emissions from each source using the detailed farm data.

### Enteric Fermentation Emissions

```{r}
# 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
```

### Manure Management Emissions

```{r}
# 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)
```

### Soil N2O Emissions

```{r}
# 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)
```

### Energy-Related Emissions

```{r}
# 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)
```

### Purchased Input Emissions

```{r}
# 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)
```

## Step 4: Total Emissions and Analysis

### Aggregate All Emission Sources

```{r}
# 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
```

### Visualize Emission Sources

```{r 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))
```

## Step 5: Intensity Calculations

### Milk Intensity

```{r}
# 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)
```

### Area Intensity

```{r}
# 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)
```

### Benchmarking Against Regional Standards

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

print(area_benchmark$benchmarking)
```

## Step 6: Detailed Results Analysis

### Per-Animal Emissions

```{r}
# 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")
```

### Feed Efficiency Analysis

```{r}
# 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")
```

## Step 7: Mitigation Scenario Analysis

Let's analyze potential mitigation strategies:

### Scenario 1: Improved Feed Efficiency

```{r}
# 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 2: Enhanced Manure Management

```{r}
# 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")
```

## Step 8: Comprehensive Results Visualization

### Multi-Source Emissions Chart

```{r 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")
```

### Intensity Metrics Dashboard

```{r}
# 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")
```

## Conclusion

This detailed single-farm analysis demonstrates the comprehensive capabilities of cowfootR for dairy farm carbon footprint assessment. The modular approach allows for detailed investigation of each emission source while maintaining consistency with international standards.

Key takeaways from this analysis:

- **Systematic approach**: Following the step-by-step methodology ensures completeness and accuracy
- **Data quality matters**: More detailed farm-specific data leads to more accurate results
- **Benchmarking provides context**: Comparing results against regional standards helps identify improvement opportunities
- **Scenario analysis**: Testing mitigation strategies helps prioritize actions
- **Uncertainty awareness**: Understanding data limitations guides decision-making

For processing multiple farms simultaneously, see the "Batch Farm Assessment" vignette. For methodology details, consult the "IPCC Tier Comparison" vignette.

---

*This analysis used cowfootR version `r packageVersion("cowfootR")` following IDF 2022 and IPCC 2019 standards.*