FGV 
# Global chunk options
knitr::opts_chunk$set(warning = FALSE, message = FALSE, fig.width = 12, fig.height = 6, out.width = "100%")

# Packages
library(sf); library(sp); library(dplyr); library(raster)
library(madrat); library(magclass); library(magpie4); library(luplot); library(gdx2)
library(terra); library(ggplot2); library(scales); library(geobr); library(rmdformats)

# Default palette for maps
pal_magpie <- hcl.colors(20, "YlOrBr", rev = TRUE)
# Cluster map and BRA filter with numeric lon/lat
map_file <- readRDS("clustermap_rev4.121_c200_67420_h12.rds")

BRA <- subset(map_file, country == "BRA")
coords <- strsplit(BRA$cell, "\\.")
mat <- do.call(rbind, coords)
BRA$lon <- as.numeric(gsub("p", ".", mat[,1]))
BRA$lat <- as.numeric(gsub("p", ".", mat[,2]))
BRA$iso <- mat[,3]

# Shapes
bra    <- geobr::read_country(showProgress = FALSE) |> st_transform(4326)
biomes <- geobr::read_biomes(showProgress = FALSE)  |> st_transform(4326)
state <- geobr::read_state(showProgress = FALSE) |> st_transform(4326)
# Formats "2010" -> "y2010.y2010"
fmt_year <- function(y) paste0("y", y, ".y", y)

# Converts magclass -> data.frame with consistent names
mag_to_df <- function(obj) {
  df <- as.data.frame(obj, rev = 2)
  names(df) <- c("lon", "lat", "iso", "year.start", "year.end", "emis_source", "pollutants", "value")
  df
}

# Standardized GWP (SAR 100y)
gwp_sar <- c(co2 = 1, ch4 = 21, n2o = 310)

mt_to_Gg <- function(x) x * 1000
mt_to_kt <- function(x) x * 1000

# Applies GWP and converts to Gg CO2e/year
apply_gwp <- function(df, gas = "ch4", gwp_table = gwp_sar) {
  g <- tolower(gas)
  if (!g %in% names(gwp_table)) stop("Gas not found in GWP table.")
  df |>
    mutate(value_co2e = value * gwp_table[[g]])
}

# Rasterizes a data.frame (lon,lat,value) and returns a grid per year (stacked list)
gridify_by_year <- function(d, bra_sf) {
  split(d, d$year.end) |>
    lapply(function(x) {
      r  <- terra::rast(x[, c("lon","lat","value")], type = "xyz", crs = "EPSG:4326")
      r  <- terra::crop(r, terra::vect(bra_sf))
      r  <- terra::mask(r, terra::vect(bra_sf))
      gd <- as.data.frame(r, xy = TRUE, na.rm = TRUE)
      names(gd) <- c("lon","lat","value")
      gd$year.end <- unique(x$year.end)
      gd
    }) |>
    dplyr::bind_rows()
}

# Sums the total by source within a subset (country, yearS) and returns top-k names (dim k of magclass)
top_k_sources <- function(dis_obj, country = "BRA", yearS, k = 2, pattern = NULL) {
  ks <- magclass::getItems(dis_obj, dim = 3)
  if (!is.null(pattern)) ks <- grep(pattern, ks, value = TRUE)
  totals <- sapply(ks, function(p) sum(as.numeric(dis_obj[country, yearS, p]), na.rm = TRUE))
  totals <- totals[totals > 0]
  names(sort(totals, decreasing = TRUE))[seq_len(min(k, length(totals)))]
}

# Extracts a specific source and returns a df with lon/lat/iso/years/value (Mt gas), plus 'source'
extract_one_source <- function(dis_obj, country, yearS, src) {
  obj <- dis_obj[country, yearS, src]
  d   <- mag_to_df(obj)
  d$lon   <- as.numeric(gsub("p", ".", d$lon))
  d$lat   <- as.numeric(gsub("p", ".", d$lat))
  d$source <- src
  d
}

# Map figure faceted by year
plot_faceted_map <- function(gdf, bra_sf, biomes_sf, title, fill_lab = "kt CO2e", pal = pal_magpie) {
  ggplot() +
    geom_raster(data = gdf, aes(lon, lat, fill = value)) +
    geom_sf(data = biomes_sf, fill = NA, color = "grey35", linewidth = 0.25) +
    geom_sf(data = bra_sf,    fill = NA, color = "grey15", linewidth = 0.4) +
    coord_sf(xlim = c(-74, -34), ylim = c(-34, 6), expand = FALSE) +
    facet_wrap(~ year.end) +
    scale_fill_gradientn(colors = pal, name = fill_lab) +
    labs(title = title, caption = "Source: MAgPIE / Default version (SSP2)") +
    theme_minimal(base_size = 12) +
    theme(
      legend.position = "right",
      panel.grid      = element_blank(),
      plot.title      = element_text(face = "bold", size = 26, hjust = 0)
    )
}

# Map faceted by source (single year)
plot_by_source <- function(gdf, bra_sf, biomes_sf, title, fill_lab = "kt CO2e", pal = pal_magpie, ncol = 2) {
  ggplot() +
    geom_raster(data = gdf, aes(lon, lat, fill = value)) +
    geom_sf(data = biomes_sf, fill = NA, color = "grey35", linewidth = 0.25) +
    geom_sf(data = bra_sf,    fill = NA, color = "grey15", linewidth = 0.4) +
    coord_sf(xlim = c(-74, -34), ylim = c(-34, 6), expand = FALSE) +
    facet_wrap(~ source, ncol = ncol) +
    scale_fill_gradientn(colors = pal, name = fill_lab) +
    labs(title = title, caption = "Source: MAgPIE / Default version (SSP2)") +
    theme_minimal(base_size = 12) +
    theme(
      legend.position = "right",
      panel.grid      = element_blank(),
      plot.title      = element_text(face = "bold", size = 22)
    )
}

Introduction

Applying global land-use and emissions models in national contexts requires caution, since each country presents its own dynamics in terms of productive activities, biomes, and emission patterns. In the case of Brazil, marked by a strong presence of the agricultural sector and by the diversity of its ecosystems, incorporating these specificities is essential to ensure the consistency of the analyses.

MAgPIE (Model of Agricultural Production and its Impact on the Environment) is a global partial equilibrium model widely used to simulate interactions among land use, agricultural production, international trade, and greenhouse gas (GHG) emissions. However, its default configuration does not include specific adjustments for the Brazilian context, which raises the need to assess how closely the resulting outputs align with officially reported values.

This report seeks to evaluate the adherence of MAgPIE, in its default configuration, to the emissions data made available by Brazil in its national GHG inventories. The analysis focuses on the main gases associated with the agricultural sector and land use (methane (\(CH_4\)), nitrous oxide (\(N_2 O\)), and carbon dioxide (\(CO_2\)) for the years 2010, 2015, and 2020, as reported in the Fourth National Inventory of Anthropogenic Emissions and Removals, submitted to the United Nations Framework Convention on Climate Change (UNFCCC).

The comparison aims to assess how closely the trajectories simulated by the model align with observed data, highlighting potential discrepancies and relevant spatial patterns. Based on this, the goal is to identify pathways for methodological adjustments that could improve the use of MAgPIE in the Brazilian context, ensuring greater alignment with national specificities and greater robustness in future analyses.

Methods

To obtain the results presented below, MAgPIE was run under the default version, using the generated GDX file (fulldata.gdx) as the basis for the analysis. Extraction in R was performed using functions provided by the magpie4 package, particularly the Emissions() function, which allows greenhouse gas emissions to be retrieved directly from the GDX by gas type, aggregation level, and unit. This function returns an object in magclass format, in which the dimensions correspond to regions, time periods, and emission sources. In this report, the type parameter was set according to the gas of interest ("ch4", "n2o", "co2"), the level parameter was set to "reg" (MAgPIE’s regional level), and the unit chosen was "GWP100AR5", ensuring comparability with the official values of the National Inventory, since in that case emissions are reported in \(CO_2\) equivalent (\(CO_2e\)) according to the AR5 100-year equivalence factors.

From this regional-level extraction, it was necessary to redistribute the results to the cell level, comprising 2,901 cells within Brazilian territory. This conversion was performed using the toolAggregate() function from the madrat package, which uses a mapping between regions and cells combined with a spatial weight. In this report, the variable chosen as weight was pasture area, in order to reflect the spatial distribution of agricultural emissions.

Regarding emission categories, MAgPIE’s structure differentiates sources by source.gas pairs (for example, ent_ferm.ch4 for enteric fermentation, awms.ch4 for animal waste management, or rice.ch4 for irrigated rice). This organization allowed both the summation of emissions to obtain the total for each gas and the decomposition by main sources. In addition, by applying a conversion from megatonnes (Mt) to kilotonnes (kt), the results were adjusted so they could be directly compared to the observed values from the Fourth National Inventory, as reported by SIRENE.

Results

The results for methane and nitrous oxide are presented below.

Methane

The figure below shows the spatial distribution of total \(CH_4\) emissions in Brazil (in kt \(CO_2e\), GWP100-AR5) for 2010, 2015, and 2020, obtained from the MAgPIE model estimation, default version. Each panel corresponds to one year; each grid cell results from the redistribution of regional emissions to the cell level, weighted by pasture area. The color scale (darker shades = higher values) is common to all three panels, allowing visual comparison across years. The map shows the total for the gas (sum of sources), without sectoral decomposition, overlaid on biome boundaries.

# Land use (pasture as weight)
x <- read.magpie("cell.land_0.5.mz")[,, "past"]
x[is.na(x)] <- 0

gas    <- "ch4"
years  <- c(2010, 2015, 2020)
yearS  <- fmt_year(years)                 # "y2010.y2010", ...
 
# 1) Emissions by gas at regional level (reg), in Mt gas
emis_gas      <- Emissions(gdx = "fulldata.gdx", type = gas, level = "reg", unit = "GWP100AR5", cumulative = FALSE)
emis_gas <- emis_gas[,c("y2010", "y2015", "y2020"),]

# 2) Aggregate from region -> cell using the map/weight
dis_emis_gas  <- madrat::toolAggregate(emis_gas, map_file, x, from = "region", to = "cell")

# 3) For all 'pollutants' of the gas, extract BRA x years and stack
ks <- magclass::getItems(dis_emis_gas, dim = 3)
ks <- grep(paste0("\\.", gas, "$"), ks, value = TRUE) 

dfs <- lapply(ks, function(k) mag_to_df(dis_emis_gas["BRA", yearS, k]))
df_gas_total <- dplyr::bind_rows(dfs)

# 4) Sum across sources to get the total per cell/year
df_gas_sum <- aggregate(value ~ lon + lat + iso + year.start + year.end,
                        data = df_gas_total, FUN = sum)

# 5) Join with BRA grid, apply GWP (SAR) and prepare for gridding
df <- BRA |>
  left_join(df_gas_sum, by = c("lon","lat","iso")) |>
  mutate(kt_co2eq = mt_to_kt(value))

df_map <- df |>
  dplyr::select(lon, lat, year.end, value = kt_co2eq) |>
  filter(is.finite(lon), is.finite(lat), is.finite(value))

# 6) Rasterize by year and plot
g_year <- gridify_by_year(df_map, bra_sf = bra)
p_ch4  <- plot_faceted_map(g_year, bra, biomes, title = "CH4 - Brazil (kt CO2e, GWP100-AR5)")
p_ch4
Spatial distribution of $CH_4$ (kt $CO_2e$, GWP100AR5) — 2010–2020.

Spatial distribution of \(CH_4\) (kt \(CO_2e\), GWP100AR5) — 2010–2020.

The figure below shows the spatial distribution of \(CH_4\) by source for 2020 (in kt CO₂e, GWP100-AR5), highlighting the two main contributions according to MAgPIE for Brazil. Each panel corresponds to a specific source among the top emitters in the selected year (for example, enteric fermentation and animal waste management/“AWMS”, as returned by the model), with the same color scale across panels (darker shades = higher emissions).

# Find top-2 (BRA, 2020)
year_one <- fmt_year(2020)
src_top2 <- top_k_sources(dis_emis_gas, country = "BRA", yearS = year_one, k = 2,
                          pattern = paste0("\\.", gas, "$"))

# Extract, apply GWP and rasterize by source
df_src <- lapply(src_top2, function(s) extract_one_source(dis_emis_gas, "BRA", year_one, s)) |>
  dplyr::bind_rows() |>
  mutate(kt_co2eq = mt_to_kt(value)) |>
  transmute(lon, lat, source, value = kt_co2eq) |>
  filter(is.finite(lon), is.finite(lat), is.finite(value))

g_src <- split(df_src, df_src$source) |>
  lapply(function(d) {
    r  <- terra::rast(d[, c("lon","lat","value")], type = "xyz", crs = "EPSG:4326")
    r  <- terra::crop(r, terra::vect(bra))
    r  <- terra::mask(r, terra::vect(bra))
    out <- as.data.frame(r, xy = TRUE, na.rm = TRUE)
    names(out) <- c("lon","lat","value")
    out$source <- unique(d$source)
    out
  }) |>
  dplyr::bind_rows()

labels_dict <- c(
  "awms.ch4"     = "Animal waste management",
  "ent_ferm.ch4" = "Enteric fermentation",
  "rice.ch4"     = "Irrigated rice"
)
lab <- labels_dict[src_top2]; lab[is.na(lab)] <- src_top2[is.na(lab)]
g_src$source <- factor(g_src$source, levels = src_top2, labels = lab)

p_ch4_src <- plot_by_source(g_src, bra, biomes, title = "CH4 - Brazil (kt CO2e, GWP100-AR5) — main sources (2020)")
p_ch4_src
Spatial distribution of $CH_4$ considering the main emission sources (kt $CO_2e$, GWP100AR5) — 2020.

Spatial distribution of \(CH_4\) considering the main emission sources (kt \(CO_2e\), GWP100AR5) — 2020.

The figure presents a comparison between total \(CH_4\) emissions in Brazil, as estimated by the MAgPIE model in its default version (SSP2), and the official values reported in the Fourth National Emissions Inventory, for the years 2010, 2015, and 2020.

# CH4 totals for Brazil from MAgPIE (kt CO2e, AR5)
magpie_br <- df_map |>
  dplyr::mutate(year = as.integer(gsub("[^0-9]", "", year.end))) |>
  dplyr::filter(year %in% c(2010, 2015, 2020)) |>
  dplyr::group_by(year) |>
  dplyr::summarise(value_kt = sum(value, na.rm = TRUE), .groups = "drop") |>
  dplyr::mutate(source = "MAgPIE")

# CH4 totals for Brazil from SIRENE
sirene_br <- tibble::tibble(
  year = c(2010, 2015, 2020),
  value_kt = c(432565, 434079, 449829),
  source = "SIRENE"
)

# Merge and plot
df_comp <- dplyr::bind_rows(magpie_br, sirene_br)

ggplot(df_comp, aes(x = factor(year), y = value_kt, fill = source)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.75) +
  geom_text(aes(label = scales::comma(value_kt)),
            position = position_dodge(width = 0.8),
            vjust = -0.2, size = 3.5) +
  scale_y_continuous(labels = scales::label_comma()) +
  scale_fill_manual(values = c("MAgPIE" = pal_magpie[5], "SIRENE" = pal_magpie[15])) +
  labs(
    title = "CH4 — Brazil (kt CO2e, GWP100-AR5): MAgPIE vs. SIRENE)",
    x = NULL, y = "kt CO2e",
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")
Total $CH_4$ emissions in Brazil (kt $CO_2e$, GWP100-AR5) in 2010, 2015, and 2020, comparing the MAgPIE model results (default scenario, SSP2) with the values reported by the Fourth National Emissions Inventory (SIRENE).

Total \(CH_4\) emissions in Brazil (kt \(CO_2e\), GWP100-AR5) in 2010, 2015, and 2020, comparing the MAgPIE model results (default scenario, SSP2) with the values reported by the Fourth National Emissions Inventory (SIRENE).

Nitrous Oxide

The figure below shows the spatial distribution of total \(N_2O\) emissions in Brazil (in kt CO₂e, GWP100–AR5), as estimated by the MAgPIE model in its default configuration, for the years 2010, 2015, and 2020. Unlike the methane emissions case, each grid cell results from the redistribution of regional emissions to the cell level, weighted by cropland area.

According to MAgPIE’s estimates, the highest concentration of \(N_2O\) emissions is located mainly in the Center-South region, covering areas of the states of Mato Grosso do Sul, Paraná, São Paulo, Minas Gerais, and part of Goiás.

# Land use (pasture as weight)
x <- read.magpie("cell.land_0.5.mz")[,, "crop"]
x[is.na(x)] <- 0

pal_magpie <- hcl.colors(20, "Reds", rev = TRUE)

gas    <- "n2o"
years  <- c(2010, 2015, 2020)
yearS  <- fmt_year(years)                 # "y2010.y2010", ...

# 1) Emissions by gas at regional level (reg), in Mt gas
emis_gas      <- Emissions(gdx = "fulldata.gdx", type = gas, level = "reg", unit = "GWP100AR5")
emis_gas <- emis_gas[,c("y2010", "y2015", "y2020"),]

# 2) Aggregate from region -> cell using the map/weight
dis_emis_gas  <- madrat::toolAggregate(emis_gas, map_file, x, from = "region", to = "cell")

# 3) For all 'pollutants' of the gas, extract BRA x years and stack
ks <- magclass::getItems(dis_emis_gas, dim = 3)
ks <- grep(paste0("\\.", gas, "$"), ks, value = TRUE) 

dfs_n <- lapply(ks, function(k) mag_to_df(dis_emis_gas["BRA", yearS, k]))
df_gas_total <- dplyr::bind_rows(dfs_n)

# 4) Sum across sources to get the total per cell/year
df_gas_sum <- aggregate(value ~ lon + lat + iso + year.start + year.end,
                        data = df_gas_total, FUN = sum)

# 5) Join with BRA grid, apply GWP (SAR) and prepare for gridding
df_n2o <- BRA |>
  left_join(df_gas_sum, by = c("lon","lat","iso")) |>
  mutate(kt_co2eq = mt_to_kt(value))

df_map_n2o <- df_n2o |>
  dplyr::select(lon, lat, year.end, value = kt_co2eq) |>
  filter(is.finite(lon), is.finite(lat), is.finite(value))

# 6) Rasterize by year and plot
g_year_n2o <- gridify_by_year(df_map_n2o, bra_sf = bra)
p_n2o  <- plot_faceted_map(g_year_n2o, bra, biomes, title = "N2O - Brazil (kt CO2e, GWP100-AR5")
p_n2o
Spatial distribution of $N_2O$ (kt $CO_2e$, GWP100AR5) — 2010–2020.

Spatial distribution of \(N_2O\) (kt \(CO_2e\), GWP100AR5) — 2010–2020.

The model incorporates multiple \(N_2O\) sources (33 sources), including: soil organic matter, inorganic fertilizers and manure applied to crops and pastures, animal waste management, decomposition and burning of agricultural residues, rice cultivation, as well as nitrous oxide from plant carbon, litter and other vegetation, etc.

The figure below shows the spatial distribution of \(N_2O\) in Brazil in 2020 (kt CO2e, GWP100-AR5), considering the two main emission sources according to the MAgPIE model: pasture manure and animal waste management. According to the estimates, in 2020 Brazil would have emitted 113,298.6 ktCO_2e of nitrous oxide. Of this total, 26% is attributed to pasture manure and 20% to animal waste management.

The databases used include FAO, the International Fertilizer Association (IFA), and specialized scientific literature, with methodological details described in Bodirsky et al. (2012).

# Find top-2 n2o (BRA, 2020)
year_one_n <- fmt_year(2020)
src_top2_n <- top_k_sources(dis_emis_gas, country = "BRA", yearS = year_one_n, k = 2,
                          pattern = paste0("\\.", gas, "$"))

# Extract, apply GWP and rasterize by source
df_src_n <- lapply(src_top2_n, function(s) extract_one_source(dis_emis_gas, "BRA", 
                                                              year_one_n, s)) |>
  dplyr::bind_rows() |>
  mutate(kt_co2eq = mt_to_kt(value)) |>
  transmute(lon, lat, source, value = kt_co2eq) |>
  filter(is.finite(lon), is.finite(lat), is.finite(value))

g_src_n <- split(df_src_n, df_src_n$source) |>
  lapply(function(d) {
    r  <- terra::rast(d[, c("lon","lat","value")], type = "xyz", crs = "EPSG:4326")
    r  <- terra::crop(r, terra::vect(bra))
    r  <- terra::mask(r, terra::vect(bra))
    out <- as.data.frame(r, xy = TRUE, na.rm = TRUE)
    names(out) <- c("lon","lat","value")
    out$source <- unique(d$source)
    out
  }) |>
  dplyr::bind_rows()

labels_dict <- c(
  "man_past.n2o" = "Pasture manure",
  "awms.n2o"  = "Animal waste management"
  )

lab <- labels_dict[src_top2_n]; lab[is.na(lab)] <- src_top2_n[is.na(lab)]
g_src_n$source <- factor(g_src_n$source, levels = src_top2_n, labels = lab)

p_n2o_src <- plot_by_source(g_src_n, bra, biomes, title = "N2O - Brazil (kt CO2e, GWP100-AR5) - main sources (2020)")
p_n2o_src
Spatial distribution of $N_2O$ considering the main emission sources (kt $CO_2e$, GWP100AR5) — 2020.

Spatial distribution of \(N_2O\) considering the main emission sources (kt \(CO_2e\), GWP100AR5) — 2020.

The figure compares the total \(N_2O\) emissions in Brazil estimated by the MAgPIE model (default version) with the official values from the Fourth National Emissions Inventory for the years 2010, 2015, and 2020. Overall, the estimates were relatively close. However, it is important to note that the disaggregation of emission-related activities is done differently in each case, which can lead to different conclusions regarding the main emission sources.

In SIRENE, within the agricultural sector, the main emission source is agricultural soil management, which includes synthetic fertilizers, organic fertilizers, direct manure deposition, agricultural residues, nitrogen mineralization associated with soil carbon loss, and management of organic soils. Pasture manure, identified by MAgPIE as one of the main sources, is included in this grouping.

The databases and references used by MAgPIE and SIRENE differ significantly, although both adopt emission factors in accordance with IPCC (2006) guidelines.

# N2O totals for Brazil from MAgPIE (kt CO2e, AR5)
magpie_br_n <- df_map_n2o |>
  dplyr::mutate(year = as.integer(gsub("[^0-9]", "", year.end))) |>
  dplyr::filter(year %in% c(2010, 2015, 2020)) |>
  dplyr::group_by(year) |>
  dplyr::summarise(value_kt = sum(value, na.rm = TRUE), .groups = "drop") |>
  dplyr::mutate(source = "MAgPIE")

# N2O totals for Brazil from SIRENE
sirene_br_n <- tibble::tibble(
  year = c(2010, 2015, 2020),
  value_kt = c(130140,  140231,  156306), # excludes energy, waste, industrial processes
  source = "SIRENE"
  )

# Merge and plot
df_comp_n <- dplyr::bind_rows(magpie_br_n, sirene_br_n)

ggplot(df_comp_n, aes(x = factor(year), y = value_kt, fill = source)) +
  geom_col(position = position_dodge(width = 0.8), width = 0.75) +
  geom_text(aes(label = scales::comma(value_kt)),
            position = position_dodge(width = 0.8),
            vjust = -0.2, size = 3.5) +
  scale_y_continuous(labels = scales::label_comma()) +
  scale_fill_manual(values = c("MAgPIE" = pal_magpie[5], "SIRENE" = pal_magpie[15])) +
  labs(
    title = "N20 - Brazil (kt CO2e, GWP100-AR5): MAgPIE vs. SIRENE",
    x = NULL, y = "kt CO2e",
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top")
Total $N_2O$ emissions in Brazil (kt $CO_2e$, GWP100-AR5) in 2010, 2015, and 2020, comparing the MAgPIE model results (default scenario, SSP2) with the values reported by the Fourth National Emissions Inventory (SIRENE).

Total \(N_2O\) emissions in Brazil (kt \(CO_2e\), GWP100-AR5) in 2010, 2015, and 2020, comparing the MAgPIE model results (default scenario, SSP2) with the values reported by the Fourth National Emissions Inventory (SIRENE).