Ruminant_feed_balance_modelling.qmd

---
title: ""
format:
  pdf:
    cite-method: natbib
    include-in-header:
      - text: |
          \usepackage{afterpage}
          \def\blankpage{%
          \clearpage%
          \thispagestyle{empty}%
          \addtocounter{page}{-1}%
          \null%
          \clearpage}
    toc: true 
    number-sections: true
    number-depth: 3
    number-by: equations
crossref:
  lof-title: "List of Figures"
bibliography: references.bib
link-citations: true
header-includes:
  \usepackage{fancyhdr}
  \usepackage{graphicx}
  \usepackage{pdflscape}    
  \bibliographystyle{apalike} 
  \usepackage{hyperref}    
  \usepackage{eso-pic}
  \usepackage{threeparttable}
  \usepackage{float}
  \usepackage{sectsty}
  \usepackage{amsmath}
  \usepackage{tikz}
  \usepackage[round]{natbib}
  \usepackage{fvextra}
  \usepackage{booktabs}
  \DefineVerbatimEnvironment{Highlighting}{Verbatim}{breaklines,commandchars=\\\{\}}
  \floatplacement{table}{H}
  \renewcommand{\thesection}{\arabic{section}}
  \AtBeginDocument{
    \thispagestyle{empty}
      \begin{tikzpicture}[remember picture,overlay]
        \node at (current page.center) [yshift=4cm] [anchor=north,yshift=0cm] {\Huge\textbf{Ruminant Feed Balance Modelling in R}};
        \node at (current page.south) [yshift=8cm] [anchor=north,yshift=0cm] {\normalsize\textbf{January 2025}};
        \node at (current page.south) [yshift=3cm] [anchor=south,yshift=-1cm] {\includegraphics[width=18cm]{img/logos.png}};
      \end{tikzpicture}
      \clearpage}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(tidy = TRUE, width = 90)
```

\newpage

# Introduction

This course is designed to guide learners on ruminant feed balance modelling in R. The course gives an overview of the concept of feed balance modelling, and introduces the main modelling steps based by \citet{fravalFeedBalancesRuminant2024}. Learners are taught how to gather spatial and non-spatial data, pre-process it, and ultimately assess feed balances. The methodology has be implemented in Ethiopia, Burkina Faso and Nigeria. The codes and data in these course are for Nigeria. The course emphasizes practical learning, with learners working on hands-on sessions. Detailed explanations are provided, helping learners follow along with the material and providing a valuable resource for future reference.

\newpage

# Background

Feed balances are used to evaluate the adequacy of available livestock feed resources in meeting the dietary requirements of livestock \citep{mottetFeedBalancesSheet2024}. Feed resources typically include natural grasses, browse, crop by-products (e.g., crop residues), and agro-industrial by-products. Feed balance assessments are conducted at specific geographical scales and over defined time periods. In this course, we focus on estimating feed (energy) balance for ruminant livestock at the national scale, while ensuring relevance at sub-national levels and across different time frames.

## Objectives

The objectives of this course are to:-

-   Teach participants where to collect spatial and non-spatial data required for livestock feed balance assessment for Nigeria.
-   Teach participants how to import, edit, and export spatial and non-spatial data in preparation for feed balance modeling.
-   Guide participants in evaluating livestock feed supply and determining livestock nutritional requirements.
-   Equip participants with the skills to assess feed balances, allowing them to evaluate if livestock feed supply meets demand.

## Learning Outcomes

By the end of the course, learners will be able to:-

-   Locate and collect essential spatial and non-spatial data for livestock feed balance assessment for Nigeria
-   Import, edit, and export spatial and non-spatial data for use in livestock feed balance assessments.
-   Assess livestock feed supply and determine livestock nutritional requirements.
-   Evaluate feed supply against livestock requirements to assess feed balances.

## Approach

The approach involves estimating (i) feed availability using freely accessible geospatial data and (ii) ruminant livestock feed requirements, both assessed over time. @fig-FeedBal_model provides an overview of the feed balance model design.

![Feed balance model for ruminant livestock. A minus (−) indicates excluded biomass, while an asterisk (\*) denotes inclusion of feed concentration values in subsequent calculations. ME = Metabolizable Energy.](img/FeedBal_model.png){#fig-FeedBal_model fig-pos="H"}

The course is organized into three workflows: `1Data-download` for downloading data, `2Feed-geoprocessing` for estimating feed availability and `3Balance-estimates` for estimating feed balances.

\newpage

# Setting up R environment

Learners will need to install the following free and open source software:-

-   R programming language (4.4+) <https://cran.r-project.org>.
-   R Studio Desktop (3.6+) <https://rstudio.com/products/rstudio/download/#download>.

Create a new folder within your work space and name it `feed-balance-modelling`. Using R, we assign the path to the folder the variable name `root`.

```{r include = FALSE}
root <- "/home/s2255815/rdrive/AU_IBAR/ruminant-feed-balance"
```

For Linux/Unix systems

```{r eval=FALSE}
# Linux systems
root <- "/home/feed-balance-modelling"
```

For Windows systems

```{r eval=FALSE}
# Windows systems
root <- "C:/Documents/feed-balance-modelling"
```

Next, we set the country of interest.

```{r}
country <- "Nigeria"
```

\newpage

# Data collation

This tutorial aims to guide learners through the key spatial and non-spatial data sources required for assessing feed balance. Learners will be introduced to a variety of datasets, and by the end of the tutorial, they will understand the necessary data requirements and know how to gather these datasets for feed balance modeling.

The tutorial leverages freely available spatial and non-spatial data to perform feed balance modeling. Recent advancements in R packages have enabled streamlined access to spatial data, including administrative boundaries, climatic variables, and environmental datasets. The required data is listed in Table \ref{tbl:model-input-data}.

```{=tex}
\begin{landscape}
\begin{table}[h]
\caption{List of datasets required for feed balance modelling}
\label{tbl:model-input-data}
\centering
\begin{tabular}{p{3cm}p{7cm}p{2cm}p{2cm}p{2cm}p{3cm}}
\toprule
\textbf{Dataset} & \textbf{Description} & \textbf{Availability} & \textbf{Spatial resolution} & \textbf{Year} & \textbf{Source} \\
\midrule
Administrative boundaries & Administrative boundaries for target country & Global & NA & 2023 & \citet{gadmDatabaseGlobalAdministrative2024} \\ 
Aggregation zones & Areas with similar livestock management practices & Regional  & NA & 2018 & \citet{fmardEcologicalFeedDistribution2023} \\ 
Land use & Fractional cover of different land uses  & Global & 100 meter & 2019 & \citet{buchhornCopernicusGlobalLand2020a} \\
Above ground dry matter productivity & Vegetation's overall growth rate  & Global & 300 meter & 2020-2023 & \citet{copernicusDryMatterProductivity2024} \\
Crop type and area  & Location, extent, and patterns of feedable crops & Global & 10 kilometer & 2020 & \citet{ifpriGlobalSpatiallyDisaggregatedCrop2024} \\
Phenology & Seasonal pattern of variation in vegetated land surfaces observed from remote sensing & Global & 500 meter & 2020-2023 & \citet{grayUserGuideCollection2024} \\
Burned areas & Burnt scars & Global & 300 meter & 2020-2023 & \citet{padillaNearrealTimeMonitoring2024} \\
Protected areas & Marine and terrestrial protected areas & Global & NA & 2024 & \citet{unep-wcmcandiucnProtectedPlanetWorld2023} \\
Tree cover & Forest and non-forest tree cover  & Global, Regional  & 100 meter & 2019 & \citet{reinerMoreOneQuarter2023} \\
Crop harvest index & Ratio of harvested product dry weight to total above-ground biomass dry weight at plant maturity  & Local & NA & 2024 & \citet{ilriFeedAssessmentTool2022} \\
Feed parameters & Nutritional quality of feed items/type e.g., dry matter digestibility, crude protein  & Local & NA & 2024 & \citet{ilriSubSaharanAfricaFeeds2020} \\
Livestock population & Type and number of livestock  & Global & 10 kilometer & 2020 & \citet{gilbertGlobalDistributionData2018} \\
Livestock parameters & Animal characteristics e.g., live weight, age & Local & NA & 2024 & Experts \\
\bottomrule
\end{tabular}
\end{table}
\end{landscape}
```
We provide brief reproducible examples illustrating how to download and prepare such data for feed balance modeling.

## Administrative boundaries

We create a new folder under `feed-balance-modelling`, and name it `AdminBound` and assign it the variable name `outdir`. We can download administrative boundaries of world countries with the `geodata` R package. Here we use `geodata` to download the administrative boundaries for levels `0`, `1` and `2` for Nigeria from `Database of Global Administrative Areas` \citep{gadmDatabaseGlobalAdministrative2024}, and store the data in `AdminBound` folder.

```{r eval=FALSE, message=FALSE, warning=FALSE}
library(geodata)
library(sf)

country <- "Nigeria"

outdir <- paste0(root, "/src/1Data-download/SpatialData/inputs/AdminBound/", country); dir.create(outdir, F, T)

admin_levels <- c("0", "1", "2")
for(admin_level in admin_levels){
  aoi <- geodata::gadm(country = "NGA", level=admin_level, path = paste0(outdir), version="latest") %>% sf::st_as_sf()
  write_sf(aoi, paste0(outdir, "/gadm40_NGA_", admin_level, ".shp"), append = FALSE)
  write_sf(aoi, paste0(outdir, "/aoi", admin_level, ".shp"), append = FALSE)
}

```

## Aggregation zones

We use ecological and feed distribution zones as defined by \citet{fmardEcologialFeedDistribution2024}, the most recent version can be downloaded at <https://drive.google.com/file/d/10HsGspftDNgq-fjAjeUwB8QsmHZWUJyT/view?usp=sharing>.

```{r message=FALSE, warning=FALSE}
library(googledrive)

outdir <- paste0(root, "/src/1Data-download/SpatialData/inputs/AggregationZones"); dir.create(outdir, F, T)

drive_deauth()
drive_user()

public_file <- drive_get(as_id("10HsGspftDNgq-fjAjeUwB8QsmHZWUJyT"))
drive_download(public_file, path = paste0(outdir, "/Ecological_and_Feed_Distribution.zip"), overwrite = TRUE)

unzip(zipfile = paste0(outdir, "/Ecological_and_Feed_Distribution.zip"), exdir = paste0(outdir, "/"))
```

## Feed parameters

Several feed parameters are needed for model parameterization. In this course, we estimate the metabolizable energy (ME) concentration of feeds using species-specific estimates derived from literature reviews and databases such as the Sub-Saharan Africa Feeds composition Database (Link: <https://feedsdatabase.ilri.org>. Feed quality parameters include metabolizable Energy (ME), neutral detergent fibre (NDF), in vitro organic matter digestibility (IVOMD), and crude protein (CP). To facilitate model parameterization, a comprehensive feed parameter file has been prepared and is available for download at: <https://drive.google.com/drive/folders/1SpB1p9i4MGU1gMahF4M3Uc-HZr8FGoqd>

```{r eval=FALSE}
outdir <- paste0(root, "/src/1Data-download/Tables/inputs/", country, "/CropParams"); dir.create(outdir, F, T)

drive_deauth()
drive_user()

#folder link to id
public_folder <- "https://drive.google.com/drive/folders/1SpB1p9i4MGU1gMahF4M3Uc-HZr8FGoqd"
folder_id <- drive_get(as_id(public_folder))

#find files in folder
public_files <- drive_ls(folder_id)

for(i in 1:nrow(public_files)){
  public_file <- public_files[i, ]
  file_name <- public_file$name
  drive_download(public_file, path = paste0(outdir, "/", file_name), overwrite = TRUE)
}
```

Other useful websites for feed parameters include:

-   Feedipedia <https://www.feedipedia.org>
-   Tropical forages <http://www.tropicalforages.info>

## Livestock parameters

Several livestock parameters are needed for model parameterization. In this course, we utilize the Gridded Livestock of the World (GLW) database to spatially disaggregate herds and flocks into species and categories including sheep, goats (adult and young), cattle (bulls, steers, cows, heifers, and calves), horses, and donkeys. We compile essential livestock parameters to characterize herd dynamics and productivity. These parameters include: live weight, age, daily weight gain, fertility rate, lactation length, daily milk yield, annual milk yield, daily walking distance, and proportion of population used for work. To facilitate model parameterization, a comprehensive feed parameter file has been prepared and is available for download at: <https://drive.google.com/drive/folders/1-3N_kmMgcHr_tayylSxlMJAr-2PBGFXd>

```{r eval=FALSE}
outdir <- paste0(root, "/src/1Data-download/Tables/inputs/", country, "/LivestockParams"); dir.create(outdir, F, T)

drive_deauth()
drive_user()
  
#folder link to id
public_folder = "https://drive.google.com/drive/folders/1-3N_kmMgcHr_tayylSxlMJAr-2PBGFXd"
folder_id = drive_get(as_id(public_folder))

#find files in folder
public_files = drive_ls(folder_id)

for(i in 1:nrow(public_files)){
  public_file <- public_files[i, ]
  file_name <- public_file$name
  drive_download(public_file, path = paste0(outdir, "/", file_name), overwrite = TRUE)
}
```

Extended workflows and detailed scripts for collating other required datasets are available at <https://github.com/ilri/ruminant-feed-balance/tree/main/src/1Data-download>

# Feed geoprocessing

This tutorial aims to guide learners through the geo-processiong techniques used for making the datasets we collated in the previous session ready for feed balance modelling. The workflow is comprised of R scripts that should be executed in numeric-alphabetic sequence i.e., 0a, 0b, 1, 2a, 2b etc.

Here we provide brief reproducible examples illustrating execute some of the R scripts.

## Cropping rasters

The spatial datasets downloaded in the previous session included both vector and raster data. For example `dry matter productivity` herein referred to as `DMP` layers have a global extent, `aggregation zones` are at a national extent, the `tree cover` data is at a regional extent, and the `protected areas` data is provided as a vector layer.

To focus our analysis on the area of interest, we will crop the layers from their global or regional extents to a smaller extent. We will use the boundary of Nigeria as the cropping extent to refine the datasets for the analysis.

### Dry Matter Productivity

We begin by processing the `Dry Matter Productivity` layers. We create an outputs folder `outdir` to store the clipped files. We then read in the Nigeria boundary layer `aoi` to define the area of interest. Next, we list all `NetCDF` files downloaded in the previous session. We iterate through the files, cropping each to the specified area of interest, and save the resulting clipped files in the `outdir` folder.

```{r eval=FALSE}
# Load required packages
library(dplyr)
library(raster)
library(rgdal)
library(sf)

# output folder
outdir <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/Feed_DrySeason/DMP"); dir.create(outdir, F, T)

# read AOI
aoi <- read_sf(paste0(root, "/src/1Data-download/SpatialData/inputs/AdminBound/", country, "/aoi0.shp"))

nc_files <- list.files("/src/1Data-download/SpatialData/inputs/Feed/DMP", pattern = ".nc$", full.names = TRUE, recursive = TRUE)

for (nc_file in nc_files){
  
  nc_name <- gsub('.{3}$', '', basename(nc_file))
  
  iDMP <- raster::raster(nc_file, varname="DMP", ncdf=TRUE)
  iDMP <- crop(iDMP, extent(aoi))
  iDMP <- mask(iDMP, aoi)
  
  # save as GeoTIFF
  raster::writeRaster(iDMP, filename = paste0(outdir, "/", nc_name, ".tif"), overwrite=TRUE)
  
}
```

### Protected areas

As mentioned earlier, the `protected areas` layer is provided as a vector layer and will be converted into a raster format. A reference layer from the `DMP` layer collection is used for clipping and resampling during this process.

```{r eval=FALSE}
#Libraries
library(terra)

indir <- paste0(root, "/src/1Data-download/SpatialData/inputs/ProtectedAreas")
outdir <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/ProtectedAreas"); dir.create(outdir, F, T)

# load livelihoods vectors
wdpaNGA <- vect(paste0(indir, "/WDPA_WDOECM_Oct2024_Public_NGA.shp"))

# reference raster?
#r <- rast(ext(wdpaNGA), resolution = 0.00297619, crs = crs(wdpaNGA))
dmpTemp <- rast(paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/Feed_DrySeason/DMP/c_gls_DMP300-RT6_202301100000_GLOBE_OLCI_V1.1.2.tif"))

# rasterize the SpatVector
wdpaNGA <- rasterize(wdpaNGA, dmpTemp, field = "STATUS_YR")

# Write output
writeRaster(wdpaNGA, paste0(outdir, "/WDPAGlobal.tif"), overwrite=TRUE)
```

### Crop digestible fraction

We derive `crop digestible` fraction layers by calculating feedable versus non-feedable crops on using the `crop type and distribution` layers.

```{r eval=FALSE}
# Load libraries
library(terra)
library(dplyr)
library(readr)

#Runs with 16gb ram and 40+gb hdd space
terraOptions(tempdir = "/home/scratch/AUTemp")
terraOptions(memfrac=0.5)
terraOptions(todisk=TRUE)

pathLU <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/Feed_DrySeason/LandUse")
filesLU <- list.files(path = pathLU, pattern = ".tif$", full.names = T)

pathSPAM <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/SPAM2020")
pathSPAMInter <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/SPAM2020/intermediate"); dir.create(pathSPAMInter, F, T)

#end of file name should be physical area_cropname_a 
filesSPAM <- list.files(path = pathSPAM, pattern = "_a.tif$", full.names = T)

#stSPAM <- stack(filesSPAM)
stSPAM <- rast(filesSPAM)

###Calculate non-feed crops proportion from SPAM model.
tmpNonFeed <- read_csv(paste0(root, "/src/1Data-download/Tables/inputs/CropParams/crop_harvest index.csv")) %>% filter(Excluded != "0") %>%
  pull(codeSPAM) %>%
  unique()
tmpNonFeedIndex <- grep(pattern = paste(tmpNonFeed, collapse = "|"), names(stSPAM))
iSPAMtotalArea <- app(stSPAM, fun = sum, na.rm=TRUE)
iSPAMnonFeedArea <- app(stSPAM[[tmpNonFeedIndex]], fun = sum, na.rm=TRUE)

iSPAMnonFeedFrac <- (iSPAMnonFeedArea / iSPAMtotalArea)
iSPAMnonFeedFrac <- classify(iSPAMnonFeedFrac, cbind(NA, NA, 0), right=FALSE) # Replace missing values with 0
iSPAMAnimalDigestCropFrac <- 1 - iSPAMnonFeedFrac

iSPAMAnimalDigestCropFrac <- writeRaster(iSPAMAnimalDigestCropFrac, paste0(pathSPAMInter, "/animal_digest_frac.tif"), overwrite = T)
```

### Crop residue fraction

We derive `crop residue` fraction layers by combining data on `crop type and distribution` with `harvest index` obtained from published sources.

```{r eval=FALSE}
# Load libraries
library(dplyr)
library(raster)
library(rgdal)
library(readr)

#Runs with 16gb ram and 40+gb hdd space
rasterOptions(tmpdir="/home/scratch/AUTemp")
rasterOptions(maxmemory = 1e+60)
rasterOptions(todisk = TRUE)

#setwd("/exports/eddie/scratch/sfraval/feed-surfaces/")
cropHI <- read_csv(paste0(root, "/src/1Data-download/Tables/inputs/CropParams/crop_harvest index.csv"))

pathSPAM <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/SPAM2020")
pathSPAMInter <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/SPAM2020/intermediate")

#end of file name should be physical area_cropname_a
filesSPAM <- list.files(path = pathSPAM, pattern = "_a.tif$", full.names = T)
stSPAM <- stack(filesSPAM)
gc()

crops <-sub(".*_a_(.*?)_a\\.tif$", "\\1", filesSPAM)

tmpCropIndex <- grep(pattern = paste(crops, collapse = "|"), names(stSPAM))
iSPAMcropArea <- calc(stack(stSPAM[[tmpCropIndex]]), fun = sum, na.rm = TRUE)

print("past 1")

stHI <- stack()
stSPAMcropProp <- stack()
gc()
for(i in 1: length(crops)){
  tmpCropIndex <- grep(pattern = paste(crops[i], collapse = "|"), names(stSPAM))
  iSPAMtmpArea <- overlay(stSPAM[[tmpCropIndex]], fun = sum, na.rm=TRUE)
  icrop <- stSPAM[[tmpCropIndex]]
  icrop[icrop >0] <- (1 - cropHI$harvest_index[cropHI$codeSPAM == crops[i]])
  stHI <- stack(stHI, icrop)
  
  stSPAMcropProp <- stack(stSPAMcropProp, overlay(iSPAMtmpArea, iSPAMcropArea, fun = function(x, y){(x/y)}))
  
  print(paste("Loop", i))
}
gc()

iSPAMcropResFrac <- weighted.mean(stHI, stSPAMcropProp, na.rm = T)
iSPAMcropResFrac <- reclassify(iSPAMcropResFrac, cbind(NA, NA, 0.8), right=FALSE) #Assume that 80% is available for animals

print("past mean")

writeRaster(iSPAMcropResFrac, paste0(pathSPAMInter, "/crop_res_frac.tif"), overwrite = T)
```

We calculate the proportion of feed items i.e., cereals, roots, legumes, and oil crops using the `crop area and distribution` data. We write the new layers in the `SPAM2020` folder.

```{r eval=FALSE}
# Load libraries
library(terra)
library(readr)

terraOptions(tempdir = "/home/scratch/AUTemp")
terraOptions(memfrac=0.5)
terraOptions(todisk=TRUE)

spamPath <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/SPAM2020")

cropLookup <- read_csv(paste0(root, "/src/1Data-download/Tables/inputs/CropParams/Crop classification_feed basket.csv"))
filesSPAM <- list.files(path = spamPath, pattern = "_a.tif$", full.names = T)
stCrops <- rast(filesSPAM)

# Index major feed types
indexCere <- grep(pattern = paste(cropLookup$codeSPAM[cropLookup$codeBasket_grouped_NGA == "cere"], collapse = "|"), sub(".*_a_(.*?)_a\\.tif$", "\\1", filesSPAM))
indexRoots <- grep(pattern = paste(cropLookup$codeSPAM[cropLookup$codeBasket_grouped_NGA == "roots"], collapse = "|"), sub(".*_a_(.*?)_a\\.tif$", "\\1", filesSPAM))
indexLeg <- grep(pattern = paste(cropLookup$codeSPAM[cropLookup$codeBasket_grouped_NGA == "leg"], collapse = "|"), sub(".*_a_(.*?)_a\\.tif$", "\\1", filesSPAM))
indexOilc <- grep(pattern = paste(cropLookup$codeSPAM[cropLookup$codeBasket_grouped_NGA == "oilc"], collapse = "|"), sub(".*_a_(.*?)_a\\.tif$", "\\1", filesSPAM))

# Extraction area for major feed categories
areaTotal <- app(stCrops, fun = sum, na.rm = T)
areaCere <- app(stCrops[[indexCere]], fun = sum, na.rm = T)
areaRoots <- app(stCrops[[indexRoots]], fun = sum, na.rm = T)
areaLeg <- app(stCrops[[indexLeg]], fun = sum, na.rm = T)
areaOilc <- app(stCrops[[indexOilc]], fun = sum, na.rm = T)

# Calculate proportions
propCere <- areaCere / areaTotal
propRoots <- areaRoots / areaTotal
propLeg <- areaLeg / areaTotal
propOilc <- areaOilc / areaTotal

# Write outputs
writeRaster(propCere, paste0(spamPath, "/propCereSPAM.tif"), overwrite=TRUE)
writeRaster(propRoots, paste0(spamPath, "/propRootsSPAM.tif"), overwrite=TRUE)
writeRaster(propLeg, paste0(spamPath, "/propLegSPAM.tif"), overwrite=TRUE)
writeRaster(propOilc, paste0(spamPath, "/propOilcSPAM.tif"), overwrite=TRUE)
```

Additional workflows and detailed scripts for preparing other files are available at:

-   Landuse <https://github.com/ilri/ruminant-feed-balance/blob/main/src/2Feed-geoprocessing/1bPrepareLanduse_clip.R>
-   Tree cover <https://github.com/ilri/ruminant-feed-balance/blob/main/src/2Feed-geoprocessing/1cPrepareTreeCover_clip.R>
-   Digital Earth Africa land use <https://github.com/ilri/ruminant-feed-balance/blob/main/src/2Feed-geoprocessing/1dPrepareDEALanduse_clip.R>
-   Phenology <https://github.com/ilri/ruminant-feed-balance/blob/main/src/2Feed-geoprocessing/2aPreparePhenologyModis.R>
-   Burned area <https://github.com/ilri/ruminant-feed-balance/blob/main/src/2Feed-geoprocessing/3aBurnedDaysclip.R>
-   SPAM <https://github.com/ilri/ruminant-feed-balance/blob/main/src/2Feed-geoprocessing/4aSPAMclip.R>
-   Livestock population <https://github.com/ilri/ruminant-feed-balance/blob/main/src/2Feed-geoprocessing/6prepareLivestockPopulation.R>
-   Livestock production systems <https://github.com/ilri/ruminant-feed-balance/blob/main/src/2Feed-geoprocessing/7prepareLivestockSystems.R>

### Seasonal DM availability

We calculate the number of cropping days in a year by analyzing phenology data (e.g., green-up and senescence dates), over a time series `2020-2013`. We then define `wet` and `dry` periods. Next, we combine `crop residue`, `browse`, `natural grass` fractions, and `DMP` data to estimate feed availability for `dry` and `wet` season over the time series.

```{r eval=FALSE}
yearOffset <- (0*365) # Base year = 2020

# Load libraries

library(dplyr)
library(raster)
library(rgdal)

rasterOptions(tmpdir="/home/scratch/AUTemp")
rasterOptions(maxmemory = 5e+20) # 6e+10 ~51GB allowed
rasterOptions(todisk = TRUE)

# read AOI
aoi <- readOGR(paste0(root, "/src/1Data-download/SpatialData/inputs/AdminBound/", country, "/aoi0.shp"))

yearList <- c("2020", "2021", "2022", "2023")

lapply(yearList, function(year){
  
  cropOutdir <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/Cropping_days"); dir.create(cropOutdir, F, T)
  FeedQuantityOutdir <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/Feed_DrySeason/Feed_quantity/", year); dir.create(FeedQuantityOutdir, F, T)
  
  pathPhen <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/Feed_DrySeason/PhenologyModis/", year, "/outputTif")
  filesPhenology <- list.files(path = pathPhen,pattern=".tif$", full.names = T)
  
  pathDMP <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/Feed_DrySeason/DMP")
  filesDMP <- list.files(path = pathDMP, pattern = paste0("RT6_", year, ".*\\.tif$"), full.names = TRUE)
  stDMP <- stack(filesDMP)
  
  datesDMP <- sub(".*RT6_(.{8}).*", "\\1", filesDMP)
  datesDMP <- as.Date(datesDMP, "%Y%m%d")
  datesDMPdiff <- as.numeric(datesDMP - as.Date("1970/01/01")) #convert to same date format as Modis phenology
  pathLU <-  paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/Feed_DrySeason/LandUse")
  filesLU <- list.files(path = pathLU, pattern = "300.tif$", full.names = T)
  pathSPAM <- paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/SPAM2020")
  
  #end of file name should be physical area_cropname_a
  filesSPAM <- list.files(path = pathSPAM, pattern = "_a.tif$", full.names = T)
  iSPAMAnimalDigestCropFrac <- raster(paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/SPAM2020/animal_digest_frac.tif"))
  
  rProtectedAreas <- stack(paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/ProtectedAreas/WDPAGlobal.tif"))
  rNonProtectedAreas <- calc(rProtectedAreas, fun = function(x){ifelse(x == 0, 1, 0)})
  rm(rProtectedAreas)
  
  print("past protected")
  
  stLU <- stack(filesLU)
  LUcrops300DEA <- raster(paste0(pathLU, "/LUcrops300DEA.tif"))
  
  stPhen <- stack(raster(grep("phenoGreenup1.tif", filesPhenology, value=TRUE)), 
                  raster(grep("phenoSenescence1.tif", filesPhenology, value=TRUE)), 
                  raster(grep("phenoGreenup2.tif", filesPhenology, value=TRUE)), 
                  raster(grep("phenoSenescence2.tif", filesPhenology, value=TRUE)))
  gc()
  
  ##Crop land use to test area
  LUcrops300DEA <- extend(LUcrops300DEA, extent(stDMP[[1]]))
  LUcrops300DEA <- crop(LUcrops300DEA, extent(stDMP[[1]]))
  LUcrops300DEA <- mask(LUcrops300DEA, aoi)
  stLU <- extend(stLU, extent(stDMP[[1]]))
  stLU <- crop(stLU, extent(stDMP[[1]]))
  stLU <- mask(stLU, aoi)
  
  ##Revise grass and shrub area
  diffCrop <- LUcrops300DEA - stLU$LUcrops300
  stLU$LUgrassShrub300 <- sum(stLU$LUgrass300, stLU$LUshrub300, na.rm = T)
  stLU$LUgrassShrub300 <- stLU$LUgrassShrub300 - LUcrops300DEA
  
  stLU$LUcrops300 <- LUcrops300DEA
  
  stPhen$phenoGreenup2 <- calc(stPhen$phenoGreenup2, fun = function(x){ifelse(x > max(datesDMPdiff)+30, NA, x)})
  stPhen$phenoSenescence2 <- calc(stPhen$phenoSenescence2, fun = function(x){ifelse(x > max(datesDMPdiff)+30, NA, x)})
  stPhen$phenoGreenup1 <- calc(stPhen$phenoGreenup1, fun = function(x){ifelse(x < min(datesDMPdiff)-30, NA, x)})
  stPhen$phenoSenescence1 <- calc(stPhen$phenoSenescence1, fun = function(x){ifelse(x < min(datesDMPdiff)-30, NA, x)})
  gc()
  
  growing2 <- (stPhen$phenoSenescence2 - stPhen$phenoGreenup2)
  growing2 <- reclassify(growing2, c(365, Inf, 0))
  growingDays <- sum((stPhen$phenoSenescence1 - stPhen$phenoGreenup1), growing2, na.rm = T)
  growingDays <- reclassify(growingDays, c(300, Inf, 300))
  gc()
  
  writeRaster(growingDays, paste0(cropOutdir, "/croppingDays_", year, ".tif"), overwrite = T)
  
  print("past 0")
  
  names(stDMP) <- paste0("d", datesDMPdiff)
  
  stLU$LUtree300 <- reclassify(stLU$LUtree300, c(-Inf, 0, 0, 200, Inf, 0)) 
  stLU$LUtree300[is.na(stLU$LUtree300)] <- 0
  
  print("past 1")
  
  #####Estimate total DMP per ha
  grassFracDry <- 1# 0.33 #max 0.55
  grassFracWet <- 1# 0.55 #max 0.55
  browseShrubFrac <- 1 #0.38 #max 0.38
  browseForestFrac <- 1
  iResidueUtil <- 1 #max 0.6
  iSPAMHarvestResidueFrac <- raster(paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/SPAM2020/crop_res_frac.tif"))
  gc()
  
  print("past overlay 1")
  
  residueFrac <- iSPAMHarvestResidueFrac

  shrubFrac <- raster(paste0(root, "/src/2Feed-geoprocessing/SpatialData/inputs/", country, "/TreeCover/treecover300m.tif"))/100
  
  gc()
  
  funGrowingGrassWet <- function(dmp, crops, grassShrub, forest, shrubFrac, greenup, senesence, greenup2, senesence2, nonprotected) {ifelse((greenup <= datesDMPdiff[i] & senesence >= datesDMPdiff[i]) | (greenup2 <= datesDMPdiff[i] & senesence2 >= datesDMPdiff[i]), (dmp*9*grassShrub*grassFracWet*(1-shrubFrac))+(dmp*9*forest*grassFracWet*(1-shrubFrac)*nonprotected), NA) } #@feedFrac is the proportion of crops grown that have feedable residues - i.e. excluding coffee, tea, ect.
  funGrowingGrassDry <- function(dmp, crops, grassShrub, forest,  shrubFrac, greenup, senesence, greenup2, senesence2, nonprotected) {ifelse((greenup > datesDMPdiff[i]) | (senesence < datesDMPdiff[i] & senesence + 60 > datesDMPdiff[i]) | (senesence2 < datesDMPdiff[i]), (dmp*9*grassShrub*grassFracDry*(1-shrubFrac))+(dmp*9*forest*grassFracDry*(1-shrubFrac)*nonprotected), NA) } #@feedFrac is the proportion of crops grown that have feedable residues - i.e. excluding coffee, tea, ect.
  funGrowingBrowse <- function(dmp, crops, grassShrub, forest,  shrubFrac, nonprotected) {(dmp*9*grassShrub*shrubFrac*browseShrubFrac)+(dmp*9*forest*nonprotected*shrubFrac*browseForestFrac)} 
  funGrowingCrops <- function(dmp, crops, greenup, senesence, feedFrac, resFrac, greenup2, senesence2) {ifelse((greenup <= datesDMPdiff[i] & senesence >= datesDMPdiff[i]) | (greenup2 <= datesDMPdiff[i] & senesence2 >= datesDMPdiff[i]), (dmp*9*crops), NA) } #@feedFrac is the proportion of crops grown that have feedable residues - i.e. excluding coffee, tea, ect.
  funGrowingAftermath <- function(dmp, crops, greenup, senesence, greenup2, senesence2, nonprotected) {ifelse((greenup > datesDMPdiff[i]) | (senesence < datesDMPdiff[i] & senesence + 60 > datesDMPdiff[i]) | (senesence2 < datesDMPdiff[i]), (dmp*9*crops), NA) } #@feedFrac is the proportion of crops grown that have feedable residues - i.e. excluding coffee, tea, ect.
  
  for(i in 1:length(names(stDMP))){
    
    iDMPGrassGrowing <- overlay(stDMP[[i]], stLU$LUcrops300, stLU$LUgrassShrub300, stLU$LUtree300, shrubFrac, stPhen$phenoGreenup1, stPhen$phenoSenescence1, stPhen$phenoGreenup2, stPhen$phenoSenescence2, rNonProtectedAreas, fun = funGrowingGrassWet)
    writeRaster(iDMPGrassGrowing, paste0(FeedQuantityOutdir, "/grassWetDMP", datesDMP[i], ".tif"), overwrite = TRUE)  
    rm(iDMPGrassGrowing)
    gc()
    
    iDMPGrassDry <- overlay(stDMP[[i]], stLU$LUcrops300, stLU$LUgrassShrub300, stLU$LUtree300, shrubFrac, stPhen$phenoGreenup1, stPhen$phenoSenescence1, stPhen$phenoGreenup2, stPhen$phenoSenescence2, rNonProtectedAreas, fun = funGrowingGrassDry)
    writeRaster(iDMPGrassDry, paste0(FeedQuantityOutdir, "/grassDryDMP", datesDMP[i], ".tif"), overwrite = TRUE)  
    rm(iDMPGrassDry)
    gc()
    
    iDMPBrowse <- overlay(stDMP[[i]], stLU$LUcrops300, stLU$LUgrassShrub300, stLU$LUtree300, shrubFrac, rNonProtectedAreas, fun = funGrowingBrowse)
    writeRaster(iDMPBrowse, paste0(FeedQuantityOutdir, "/browseDMP", datesDMP[i], ".tif"), overwrite = TRUE)  
    rm(iDMPBrowse)
    gc()
    
    iDMPCropGrowing <- overlay(stDMP[[i]], stLU$LUcrops300, stPhen$phenoGreenup1, stPhen$phenoSenescence1, iSPAMAnimalDigestCropFrac, residueFrac, stPhen$phenoGreenup2, stPhen$phenoSenescence2, fun = funGrowingCrops) 
    writeRaster(iDMPCropGrowing, paste0(FeedQuantityOutdir, "/cropDMP", datesDMP[i], ".tif"), overwrite = TRUE)
    rm(iDMPCropGrowing)
    gc()
    
    iDMPAftermath <- overlay(stDMP[[i]], stLU$LUcrops300, stPhen$phenoGreenup1, stPhen$phenoSenescence1, stPhen$phenoGreenup2, stPhen$phenoSenescence2, fun = funGrowingAftermath)
    writeRaster(iDMPAftermath, paste0(FeedQuantityOutdir, "/aftermathDMP", datesDMP[i], ".tif"), overwrite = TRUE)  
    rm(iDMPAftermath)
    gc()
    print(paste("cycle", i))
    
  }
  
  gc()
  
  iDMPgrassWet <- stack(list.files(path = paste0(FeedQuantityOutdir), pattern="grassWet",full.names = T))
  DMPgrassmeanWet <- mean(iDMPgrassWet, na.rm = T)
  writeRaster(DMPgrassmeanWet, paste0(FeedQuantityOutdir, "/DMPgrassWetmean_", year, ".tif"), overwrite = TRUE)
  
  iDMPgrassDry <- stack(list.files(path = paste0(FeedQuantityOutdir), pattern="grassDry",full.names = T))
  DMPgrassmeanDry <- mean(iDMPgrassDry, na.rm = T)
  writeRaster(DMPgrassmeanDry, paste0(FeedQuantityOutdir, "/DMPgrassDrymean_", year, ".tif"), overwrite = TRUE)
  
  iDMPbrowse <- stack(list.files(path = paste0(FeedQuantityOutdir), pattern="browse",full.names = T))
  DMPbrowsemean <- mean(iDMPbrowse, na.rm = T)
  writeRaster(DMPbrowsemean, paste0(FeedQuantityOutdir, "/DMPbrowsemean_", year, ".tif"), overwrite = TRUE)
  
  iDMPCropGrowing <- stack(list.files(path = paste0(FeedQuantityOutdir), pattern="crop",full.names = T))
  DMPcropmean <- mean(iDMPCropGrowing, na.rm = T)
  writeRaster(DMPcropmean, paste0(FeedQuantityOutdir, "/DMPcropmean_", year, ".tif"), overwrite = TRUE)
  
  iDMPAftermath <- stack(list.files(path = paste0(FeedQuantityOutdir), pattern="aftermath",full.names = T))
  DMPAftermean <- mean(iDMPAftermath, na.rm = T)
  writeRaster(DMPAftermean, paste0(FeedQuantityOutdir, "/DMPaftermean_", year, ".tif"), overwrite = TRUE)
  gc()
  
})
```

Finally, we copy relevant data to specific folders for use in the next section using this R script: <https://github.com/ilri/ruminant-feed-balance/blob/main/src/2Feed-geoprocessing/8Copydata.R>

# Estimating feed balances

## Mapping aggregation regions and zones

The ecological and feed distribution map provided by the Federal Ministry of Agriculture and Food Security (FMAFS) contain eight zones. We map them into ecological regions and zones as shown in Table \ref{tbl:ecological-zones}.

```{=tex}
\begin{table}[h]
\caption{Ecological zones and regions}
\label{tbl:ecological-zones}
\centering
\begin{tabular}{p{5cm}p{4cm}p{4cm}}
\toprule
\textbf{Ecological and Feed Distribution Zones} & \textbf{Ecological Regions} & \textbf{Ecological Zones} \\ 
\midrule
Sahel Savannah & Dry Savannah & (Agro)pastoral sahel \\ 
Sudan Savannah & Dry Savannah & Northern mixed \\ 
Southern Guinea Savannah & Wet Savannah & Southern mixed \\ 
Northern Guinea Savannah & Wet Savannah & Central mixed \\ 
Mountain Vegetations & Wet Savannah & Central mixed \\ 
Lowland Rainfall & Forest & Forest mixed \\ 
Fresh Water Swamp Forest & Forest & Forest mixed \\ 
Mangrove & Forest & Forest mixed \\ 
\bottomrule
\end{tabular}
\end{table}
```
The R script for creating the new categories of ecological regions and zones is available at: <https://github.com/ilri/ruminant-feed-balance/blob/main/src/3Balance-estimates/Nigeria/0_PrepareRegions.R>

As shown in @fig-ecoregions, the ecological regions are classified into `Dry Savannah`, `Forest`, and `Wet Savannah`.

```{r echo=FALSE, warning=FALSE, message=FALSE, fig.cap="Ecological Regions of Nigeria", fig.pos="H"}
#| label: fig-ecoregions
library(sf)
library(ggplot2)

aoi_ECORegion <- read_sf("./src/3Balance-estimates/Nigeria/SpatialData/intermediate/regions.gpkg")
region_colors <- c("Dry Savannah" = "yellow", "Forest" = "darkgreen", "Wet Savannah" = "lightgreen")
ggplot() +
  geom_sf(data = aoi_ECORegion, aes(fill = ECORegion)) +
  scale_fill_manual(values = region_colors) +
  coord_sf(xlim = c(2.1, 15.1), ylim = c(3.8, 14.3), expand = FALSE) +
  labs(fill = "Ecological zones", title = "")
```

Similarly, as shown in @fig-ecozones, the ecological regions are further classified into `(Agro)pastoral sahel`, `Northern mixed`, `Central mixed`, `Southern mixed` and `Forest mixed`.

```{r echo=FALSE, warning=FALSE, message=FALSE, fig.cap="Ecological Zones of Nigeria", fig.pos="H"}
#| label: fig-ecozones
library(sf)
library(ggplot2)

aoi_ECOZones <- read_sf("./src/3Balance-estimates/Nigeria/SpatialData/intermediate/zones.gpkg")
zones_colors <- c("(Agro)pastoral sahel" = "yellow", "Northern mixed" = "yellow3", "Central mixed" = "green", "Southern mixed" = "lightgreen", "Forest mixed" = "darkgreen")
ggplot() +
  geom_sf(data = aoi_ECOZones, aes(fill = ECOZone)) +
  scale_fill_manual(values = zones_colors) +
  coord_sf(xlim = c(2.1, 15.1), ylim = c(3.8, 14.3), expand = FALSE) +
  labs(fill = "Ecological zones", title = "")
```

## Feed energy concentration

We merge `feedQuality` dataset with the `crop_classification` dataset using a common id, and map some specific feed items to new codes. The resulting clean `feedQuality` dataset is saved as a CSV file. The R script is available at: <https://github.com/ilri/ruminant-feed-balance/blob/main/src/3Balance-estimates/Nigeria/0_ProcessSSAfeedDB.R>

We then calculate feed metabolizable energy (ME) by combining various datasets, including crop parameters, feed quality data, and spatial data related to crop type and distribution. We process feed quality data i.e., calculate dry matter intake (DMI), adjust for other variables like CP, NDF, and IVDMD. We generate summary statistics and group the data by crop type. NB: Missing or inconsistent data is addressed by imputing values or setting fixed values for certain crops (See Table \ref{tbl:feed-quality-data}).

```{=tex}
\begin{table}[H]
\caption{Feed quality for different feed items}
\label{tbl:feed-quality-data}
\centering
\begin{tabular}{p{3cm}p{1.5cm}p{1.5cm}p{1.5cm}p{1.5cm}p{1.5cm}p{1.5cm}p{1.5cm}}
\toprule
\textbf{Feed item} & \textbf{ME (SD)} & \textbf{ME (Min)} & \textbf{ME (Max)} & \textbf{ME (Mean)} & \textbf{CP} & \textbf{NDF} & \textbf{IVDMD} \\ 
\midrule
Banana* & 0.5 & 7.7 & 9.0 & 8.3 & 5.8 & 42.1 & 51.8 \\ 
Barley* & 0.8 & 5.8 & 10.3 & 7.3 & 8.0 & 68.6 & 47.7 \\ 
Bean* & 0.6 & 7.0 & 9.8 & 8.7 & 6.6 & 68.1 & 56.1 \\ 
Sugarbeet* & 1.4 & 7.0 & 13.4 & 9.9 & 19.8 & 39.1 & 67.9 \\ 
Chickpea* & 0.6 & 6.8 & 10.8 & 8.2 & 4.7 & 57.5 & 51.6 \\ 
Lentils* & 0.5 & 7.0 & 10.3 & 8.8 & 9.2 & 51.7 & 56.2 \\ 
Pigeon pea* & 0.4 & 7.6 & 9.1 & 8.2 & 16.2 & 53.0 & 53.5 \\ 
Plantain* & 0.5 & 7.7 & 9.0 & 8.3 & 5.8 & 42.1 & 51.8 \\ 
Soybean* & 0.6 & 7.0 & 10.9 & 8.6 & 7.7 & 66.6 & 57.0 \\ 
Sugarcane* & 0.3 & 7.2 & 8.0 & 7.5 & 6.9 & 73.8 & 49.4 \\ 
Wheat* & 1.2 & 5.7 & 11.0 & 7.3 & 4.3 & 72.1 & 48.7 \\ 
Cowpea* & NA & 7.7 & 10.5 & 9.3 & 12.2 & 46.8 & 61.0 \\ 
Maize* & NA & 7.9 & 7.9 & 7.9 & 5.2 & 69.1 & 51.8 \\ 
Pearl millet* & NA & 6.7 & 8.4 & 7.6 & 5.2 & 68.6 & 50.0 \\ 
Rice* & NA & 5.9 & 8.4 & 6.8 & 5.5 & 65.6 & 49.0 \\ 
Sesame* & NA & 9.2 & 9.5 & 9.3 & 19.3 & 24.2 & 59.1 \\ 
Sorghum* & NA & 6.6 & 7.9 & 7.5 & 4.5 & 71.9 & 48.5 \\ 
Sweet potato* & NA & 8.9 & 9.6 & 9.3 & 13.3 & 48.4 & 60.6 \\ 
Small millet & NA & 6.7 & 8.4 & 7.6 & 5.2 & 68.6 & NA \\ 
Other cereals* & NA & 6.7 & 8.4 & 7.6 & 5.2 & 68.6 & NA \\ 
Other legumes* & NA & 7.7 & 10.5 & 9.3 & 5.2 & 68.6 & NA \\ 
Cassava* & NA & 4.2 & 4.2 & 4.2 & NA & NA & NA \\ 
Yam* & NA & 4.2 & 4.2 & 4.2 & NA & NA & NA \\ 
Other roots and tubers* & NA & 4.2 & 4.2 & 4.2 & NA & NA & NA \\ 
Groundnut* & 1.6 & 7.2 & 9.3 & 8.5 & NA & NA & NA \\ 
Natural pasture & NA & 6.2 & 6.8 & 6.5 & NA & NA & NA \\ 
\bottomrule
\end{tabular}
\begin{tablenotes}
\item[*] \textit{* Denotes crop residue.}
\end{tablenotes}
\end{table}
```

We calculate crop-specific harvest index (See Table \ref{tbl:crop-harvest-index}), utilization and ME values. We then weight the values by the proportion of crop per pixel and aggregate the results to administrative level 2.

```{=tex}
\begin{table}[H]
\caption{Harvest index for feedable crops}
\label{tbl:crop-harvest-index}
\centering
\begin{tabular}{p{3cm}p{1.5cm}}
\toprule
\textbf{Crop name} & \textbf{Harvest index} \\ 
\midrule
Banana  &   0.44 \\
Barley  &   0.534 \\
Bean    &   0.555 \\
Sugarbeet   &   1 \\
Chickpea    &   0.31 \\
Lentils &   0.246 \\
Pigeon Pea  &   0.294 \\
Plantain    &   0.44 \\
Soybean &   0.422 \\
Sugarcane   &   0.19 \\
Wheat   &   0.36575 \\
Cowpea  &   0.335 \\
Maize   &   0.3785 \\
Pearl millet    &   0.188 \\
Rice    &   0.4613 \\
Sesame  &   0.01 \\
Sorghum &   0.458 \\
Sweet potato    &   0.27 \\
Small Millet    &   0.391 \\
Other cereals   &   0.188 \\
Other legumes   &   0.3 \\
Cassava &   0.8696 \\
Yam &   0.24102 \\
Other roots and tubers  &   0.27 \\
Groundnut   &   0.45 \\
\bottomrule
\end{tabular}
\end{table}
```
The R script for completing the above steps is available at: <https://github.com/ilri/ruminant-feed-balance/blob/main/src/3Balance-estimates/Nigeria/0b_Prepare_cropME.R>

We analyse feed availability across multiple years (2020–2023) by integrating spatial data related to crop type and distribution, grass and browse productivity and feed utilization. Cropping days are calculated, with thresholds applied to exclude areas with insufficient cropping duration. Dry days are derived as the inverse of cropping days (365 - cropping days). Metabolizable energy (ME) is estimated for various feed sources including crop residues, grass, and browse, incorporating region-specific productivity and feed quality parameters. Feed statistics, such as mean ME values for crops, grass, and browse, are extracted and aggregated within defined regions and zones.

The R script for completing the above steps is available at: <https://github.com/ilri/ruminant-feed-balance/blob/main/src/3Balance-estimates/Nigeria/1_ProcessFeedMEregional.R>, and <https://github.com/ilri/ruminant-feed-balance/blob/main/src/3Balance-estimates/Nigeria/1b_ProcessFeedMEregional_min_max.R>

## Livestock energy requirements

### Requirements for base year

Using Gridded Livestock of the World (GLW) \citep{gilbertGlobalDistributionData2018}, livestock energy requirements are calculated by disaggregating herds and flocks into species and categories including `sheep`, `goats` (adult and young), `cattle` (bulls, steers, cows, heifers, and calves), `horses`, and `donkeys` and calculating seasonal metabolizable energy (ME) requirements for livestock, including energy requirements for `maintenance` , `locomotion`, `lactation` and `growth`, equations ([-@eq-maintenance]), ([-@eq-locomotion]), ([-@eq-lactation]), ([-@eq-growth-positive]) and ([-@eq-growth-negative]).

$$
\text{MERM}_{i,j,k} = K \times S \times M \times \frac{(0.26 \times \text{MLW}_k^{0.75}) \times \exp^{(-0.008 \times A_k)}}{(0.02 \times \text{ME}_{i,j}) + 0.5}
$$ {#eq-maintenance}

Where, \textit{MERM} is the maintenance energy requirement in megajoules per day (MJ/day) of animal \textit{k}, region \textit{i} and season \textit{j}, \textit{K} is a metabolic weight coefficient (1.3, an intermediate value for Bos taurus and Bos indicus cattle breeds), \textit{S} a constant for females and castrated cattle (1), \textit{M} a constant for presence or absence of milk in the diet, \textit{MLW} is the mean live weight of animal \textit{k} in region \textit{i} in season \textit{j}, and \textit{A} is the age of animal \textit{k} in region \textit{i} in season \textit{j} and \textit{ME} is metabolizable energy content of the diet in megajoules of metabolizable energy per kilogram of dry matter (MJ ME/kg DM) in region \textit{i} and season \textit{j}.

$$
\text{MERT}_{i,j,k} = \text{WD}_{i,j,k} \times \text{MLW}_{i,j,k} \times 0.0026
$$ {#eq-locomotion}

Where, \textit{MERT} is the locomotion energy requirement in megajoules per day (MJ/day) of animal \textit{k}, region \textit{i} and season \textit{j}, \textit{WD} is distance travelled by animal \textit{k} per day in kilometres in region \textit{i} and season \textit{j}, \textit{MLW} is the mean live weight of animal \textit{k} in region \textit{i} in season \textit{j}.

For adult cows, goats and sheep: 

$$
\text{MERL}_{i,j,k} = \frac{\text{DMY}_{i,j,k} \times \text{ECM}_{i,j,k}}{\left(0.02 \times \text{ME}_{i,j}\right) + 0.04}
$$ {#eq-lactation}

Where, \textit{MERL} is the energy requirement for lactation for animal \textit{k}, in MJ/day in region \textit{i} and season \textit{j}, \textit{DMY} is the daily milk production for animal \textit{k} in kilogram in region \textit{i} and season \textit{j}, \textit{ECM} is the energy content of milk from animal \textit{k} in region \textit{i} and season \textit{j} and \textit{ME} is metabolizable energy content of the diet in megajoules of metabolizable energy per kilogram of dry matter (MJ ME/kg DM) in region \textit{i} and season \textit{j}.

Energy expended for weight gain or loss as per equations ([-@eq-growth-positive]) and ([-@eq-growth-negative]).

$$
\text{MERG+}_{i,j,k} = \frac{\text{DWG}_{i,j,k} \times 0.92 \times \text{EC}_{i,j,k}}{0.043 \times \text{ME}_{i,j}}
$$ {#eq-growth-positive}

If energy is lost:-

$$
\text{MERG-}_{i,j,k} = \frac{\text{DWG}_{i,j,k} \times 0.92}{0.8}
$$ {#eq-growth-negative}

Where, \textit{MERG+/-} is the energy expended for live weight gain/loss for animal \textit{k}, in MJ/day in region \textit{i} and season \textit{j}, \textit{DWG} is the average daily weight gain or loss for animal \textit{k} in kilograms in region \textit{i} and season \textit{j}, EC is an assumed energy content of tissue taken of 18 MJ/kg \citep{csiroNutrientRequirementsDomesticated2007} and \textit{ME} is metabolizable energy content of the diet in megajoules of metabolizable energy per kilogram of dry matter (MJ ME/kg DM) in region \textit{i} and season \textit{j}.

Finally, total energy requirements for the base year for each livestock category, season and region are calculated. The R script for completing the above steps is available at: <https://github.com/ilri/ruminant-feed-balance/blob/main/src/3Balance-estimates/Nigeria/2_Lvst_requirements.R>, and <https://github.com/ilri/ruminant-feed-balance/blob/main/src/3Balance-estimates/Nigeria/2b_Lvst_requirements_min_max.R>

### Estimating feed balances

We sum livestock energy requirements for cattle, shoats (sheep/goats), and horses/donkeys and adjust the requirements by incorporating population changes from FAOSTAT data, equation ([-@eq-energy-requirements-species]).

$$
\text{LvstME}_y = \sum \left( H_y + \left( C_y \times (1 + \Delta P_C) \right) + \left( S_y \times (1 + \Delta P_S) \right) \right)
$$ {#eq-energy-requirements-species}

Where, \textit{LvstME} is the total energy requirement for all livestock in year \textit{y}, \textit{H} is the energy requirement for horses and donkeys in MJ in year \textit{y}, \textit{C} is the energy requirement for cattle in MJ in year \textit{y}, \textit{S} is the energy requirement for shoats in MJ in year \textit{y}, \textit{\(\Delta {P}_C\)} is the annual percentage change in cattle population and \textit{\(\Delta {P}_S\)} is the annual percentage change in shoats population.

We combine feed availability data into a single dataset for annual analysis, equation ([-@eq-energy-feed]).

$$
\text{FeedME}_y = \sum_{y=2020}^{2023} \text{MEF}_{y}
$$ {#eq-energy-feed}

Where, \textit{FeedME} is the total energy from feeds for year \textit{y}, \textit{MEF} is the energy from feeds for year \textit{y}.

We aggregate and extract feed energy concentration and livestock energy requirements data across three spatial zones: `country`, `region`, and `state`across the time series. The zones were defined by spatial data representing Nigeria’s administrative boundaries and ecological zones. Finally we calculate feed balance for each spatial zone and year as the ratio of feed energy concentration to livestock energy requirements, equation ([-@eq-feed-balance]).

$$
\text{FeedBal}_{y,i} = \frac{\text{FeedME}_{y,i}}{\text{LvstME}_{y,i}}
$$ {#eq-feed-balance}

Where, \textit{FeedBal} is the feed balance for year \textit{y} in region \textit{i}, \textit{FeedME} is the feed energy concentration for year \textit{y} in region \textit{i} and \textit{LvstME} is the livestock energy requirement for year \textit{y} in region \textit{i}.

The R script for completing the above steps is available at: <https://github.com/ilri/ruminant-feed-balance/blob/main/src/3Balance-estimates/Nigeria/3_ProcessTimeseriesRegionalAdeqTotals.R>, and <https://github.com/ilri/ruminant-feed-balance/blob/main/src/3Balance-estimates/Nigeria/3b_ProcessTimeseriesRegionalAdeqTotals_min_max.R>

### Visualising results

The total biomass available as animal feed in Nigeria was `294.08` million tonnes (Mt) of DM, averaging `2.21` tonnes DM per hectare in the dry savannah, `3.73` tonnes DM per hectare in the forest region and `3.83` tonnes DM per hectare in the wet savannah.

In the dry savannah, grass contributed `28.17%` of the total metabolizable energy (ME), with crop residue, browse, and biomass from postharvest growth contributing `46.67%`, `0.40%`, and `24.76%`, respectively. In the forest zone, grass accounted for `45.22%`, while with crop residue, browse, and biomass from postharvest growth contributing `7.95%`, `42.41%` and `4.42%`. In the wet savannah, grass provided `40.94%` of the total ME, crop residue contributed `32.80%` and `9.05%`, and biomass from postharvest growth `17.21%` (@fig-Total_ME).

![Metabolisable energy (ME) availability for ruminant livestock from feed sources in 2023 (TJ ME/year). Uncertainty due to ME concentration represented by error bands.](img/Total_ME.png){#fig-Total_ME fig-pos="H"}

The adequacy of locally available metabolisable energy for ruminant livestock was most constrained in the (Agro)pastoral sahel (68% of requirements locally available) (@fig-Requirements_me).

![Energy balance of ruminant livestock from feeds over time by zone (J ME available / J ME required/year). Uncertainty due to ME concentration and livestock parameters represented by grey bands. Dashed line at 1 where feed is sufficient for all estimated requirements (maintenance, growth, locomotion, draught power, gestation and lactation) for the zone.](img/Requirements_me.png){#fig-Requirements_me fig-pos="H"}

The R script for producing the above visualization and more is available at: <https://github.com/ilri/ruminant-feed-balance/blob/main/src/3Balance-estimates/Nigeria/4_MapAdequacy.R>

\newpage

# References

::: {#refs}
:::