La expansión forestal reciente de Concepción deja una señal territorial medible: más verdor y humedad del dosel, menor temperatura superficial y mayor evapotranspiración. Este artículo integra cartografía, teledetección, un estudio de eventos ponderado y el contexto exportador nacional, con resultados recalculados durante la generación del sitio.

1 Resumen

La expansión de plantaciones forestales sobre pasturas modifica la cobertura del suelo, el intercambio de energía y agua y la organización productiva del territorio. Este trabajo cuantifica esas dimensiones en el departamento de Concepción, Paraguay, durante 2015–2024. Se identificaron transiciones anuales mediante MapBiomas Paraguay y se validaron con la cartografía oficial del Instituto Forestal Nacional (INFONA). Se construyó un panel balanceado de 533 unidades tratadas y 533 controles de pastura persistente, con 10.660 observaciones anuales y 127.920 mensuales. Las trayectorias ambientales integraron Sentinel-2, Landsat, MODIS, CHIRPS y JRC Global Surface Water. Para reducir el desequilibrio pretratamiento se aplicaron pesos de solapamiento estimados mediante regresión logística y diferencias ponderadas de cambios por cohorte, con el año anterior a la conversión como referencia. El componente económico procesó 120 archivos mensuales y 42.663 registros físicos de exportación del capítulo 44 de la Nomenclatura Común del Mercosur.

Se validaron 533 transiciones que abarcaron 3.425,94 ha; el 87,81 % de la superficie correspondió a 2022–2024 y el 94,42 % se concentró en Sargento José Félix López y Paso Barreto. En el año de conversión, respecto de pasturas persistentes ponderadas, el NDVI aumentó 0,168, el EVI 0,119 y el NDMI 0,191; la temperatura superficial disminuyó 2,21 °C y la evapotranspiración aumentó 243,18 mm/año. Entre 2015 y 2024, el valor FOB nacional del capítulo 44 creció 36,7 % y el volumen 31,3 %, mientras el valor unitario aumentó 4,1 %. Los resultados muestran una firma biofísica consistente; la atribución causal permanece preliminar por la incertidumbre espacial, las diferencias climáticas tardías y la ausencia de procedencia subnacional en aduanas.

Palabras clave: cambio de uso del suelo; plantaciones forestales; teledetección; evapotranspiración; temperatura superficial; estudio de eventos; Paraguay.

Mensaje central. La conversión observada está asociada con un dosel más verde y húmedo, una superficie aproximadamente 2–3 °C más fría y una evapotranspiración anual 200–340 mm mayor. La señal es coherente y cuantitativamente importante; su interpretación causal requiere una etapa inferencial adicional.

2 Introducción

La forestación productiva altera simultáneamente la cobertura del suelo, los flujos de energía y agua y las cadenas de valor. En paisajes dominados por pasturas, el establecimiento de árboles aumenta la altura, rugosidad y área foliar, puede reducir la temperatura superficial diurna y modificar la partición de la precipitación. Los estudios globales también advierten que una mayor evapotranspiración puede acompañarse de menor rendimiento hídrico de cuencas (Farley et al., 2005; Jackson et al., 2005). La dirección y magnitud de la respuesta dependen del clima, suelo, edad del rodal, manejo y escala de observación.

Concepción constituye un caso relevante por la rapidez de la expansión reciente y por su estructura productiva. En 2024, el producto interno bruto departamental creció 7,3 %. El agregado de ganadería, actividad forestal, pesca y minería creció 18,3 %, aportó 2,0 puntos porcentuales al crecimiento y representó 12,1 % de la estructura regional; la manufactura representó 29,2 % (Banco Central del Paraguay, 2026). Estas categorías no aíslan la silvicultura, pero sitúan el estudio en un territorio con dinamismo primario e industrial.

El problema empírico exige distinguir plantaciones nuevas de coberturas arbóreas anteriores, construir un contrafactual razonable y evitar atribuir exportaciones nacionales a un departamento sin trazabilidad de origen. Este trabajo combina validación cartográfica, controles de pastura persistente y ponderación por soporte común. El componente aduanero se utiliza como contexto nacional y no como evidencia del origen de los productos.

2.1 Pregunta de investigación e hipótesis

La pregunta principal es: ¿qué cambios térmicos, hídricos y espectrales se observan después de la conversión de pasturas a plantaciones forestales en Concepción, frente a pasturas persistentes comparables, y qué contexto productivo ofrece la dinámica nacional de exportaciones forestales?

  1. H1, biofísica. Después de la conversión, las unidades forestadas presentan mayor NDVI, EVI y NDMI, menor temperatura superficial y mayor evapotranspiración que pasturas persistentes de soporte común.
  2. H2, productiva. Entre 2015 y 2024, las exportaciones forestales paraguayas se expanden principalmente por volumen y mantienen una composición concentrada en bioenergía y productos de transformación limitada. Esta hipótesis es descriptiva y no atribuye embarques a Concepción.

3 Objetivos

3.1 Objetivo general

Cuantificar la expansión reciente de plantaciones forestales sobre pasturas en Concepción y estimar su asociación con indicadores de vegetación, humedad, temperatura superficial y evapotranspiración, integrando un contexto económico nacional reproducible para 2015–2024.

3.2 Objetivos específicos

  1. Delimitar y validar espacialmente las transiciones anuales de pastura a plantación forestal.
  2. Construir controles de pastura persistente y un panel ambiental sin duplicaciones.
  3. Evaluar el balance pretratamiento y estimar trayectorias relativas al año de conversión.
  4. Describir el valor, volumen, valor unitario, composición y destinos de las exportaciones paraguayas del capítulo 44.
  5. Diferenciar los hallazgos respaldados por los datos de las inferencias que requieren validación adicional.

4 Materiales y métodos

4.1 Diseño, área y unidad de análisis

El área de estudio fue el departamento de Concepción, Paraguay. El periodo principal abarcó 2015–2024; la cartografía de uso del suelo incluyó 2025 para comprobar la persistencia de los controles. La unidad tratada fue un polígono con transición anual de pastura a plantación forestal, identificado con MapBiomas Paraguay y validado por intersección con el inventario del INFONA (Instituto Forestal Nacional, 2026; MapBiomas Paraguay, 2026). Por cada tratamiento se generó un control circular en el mismo distrito, sobre píxeles clasificados como pastura durante 2015–2025, fuera de una franja de 1.000 m de los tratamientos, sin solapamiento y con superficie máxima de 5 ha.

4.2 Fuentes e indicadores ambientales

El procesamiento geoespacial utilizó UTM 21S (EPSG:32721). Sentinel-2 aportó NDVI, EVI y NDMI; Landsat 8/9, temperatura superficial; MOD16A2GF, ET y PET; CHIRPS, precipitación; y JRC Global Surface Water, contexto hidrológico. La documentación metodológica de los índices y productos se apoya en Tucker (1979), Huete et al. (2002), Gao (1996), Mu et al. (2011), Funk et al. (2015) y Pekel et al. (2016). Las extracciones se ejecutaron en Google Earth Engine (Gorelick et al., 2017).

Sea ρNIR\rho_{NIR} la reflectancia del infrarrojo cercano, ρR\rho_R la del rojo, ρB\rho_B la del azul y ρSWIR\rho_{SWIR} la del infrarrojo de onda corta. Los índices fueron

NDVI=ρNIR−ρRρNIR+ρR,EVI=2.5ρNIR−ρRρNIR+6ρR−7.5ρB+1, NDVI=\frac{\rho_{NIR}-\rho_R}{\rho_{NIR}+\rho_R},\qquad EVI=2.5\frac{\rho_{NIR}-\rho_R}{\rho_{NIR}+6\rho_R-7.5\rho_B+1},

NDMI=ρNIR−ρSWIRρNIR+ρSWIR. NDMI=\frac{\rho_{NIR}-\rho_{SWIR}}{\rho_{NIR}+\rho_{SWIR}}.

Para la unidad ii y el año tt, la media zonal fue Yit=∑paipYpt/∑paipY_{it}=\sum_p a_{ip}Y_{pt}/\sum_p a_{ip}, donde aipa_{ip} es el área de intersección entre el píxel pp y la unidad.

4.3 Ponderación y estudio de eventos

El análisis incluyó las cohortes 2020–2024, con tres años pretratamiento completos para Sentinel-2. Si DiD_i indica tratamiento y 𝐗i\mathbf X_i reúne niveles, cambios y tendencias previas, superficie, distancias y efectos de distrito y cohorte, la probabilidad de tratamiento se estimó mediante

logit⁡(pi)=α+𝐗i𝖳𝛃,pi=Pr⁡(Di=1∣𝐗i). \operatorname{logit}(p_i)=\alpha+\mathbf X_i^{\mathsf T}\boldsymbol\beta, \qquad p_i=\Pr(D_i=1\mid\mathbf X_i).

Los pesos de solapamiento (Li et al., 2018) fueron

wi=Di(1−pi)+(1−Di)pi. w_i=D_i(1-p_i)+(1-D_i)p_i.

El estimando corresponde a la diferencia media en cambios de la población con soporte común. Para cohorte gg, tiempo relativo kk y referencia g−1g-1:

ΔYigk=Yi,g+k−Yi,g−1, \Delta Y_{igk}=Y_{i,g+k}-Y_{i,g-1},

μ̂dgk=∑i:Gi=g,Di=dwiΔYigk∑i:Gi=g,Di=dwi,τ̂gk=μ̂1gk−μ̂0gk. \widehat\mu_{dgk}=\frac{\sum_{i:G_i=g,D_i=d}w_i\Delta Y_{igk}} {\sum_{i:G_i=g,D_i=d}w_i},\qquad \widehat\tau_{gk}=\widehat\mu_{1gk}-\widehat\mu_{0gk}.

Las estimaciones se agregaron como τ̂k=∑gωgkτ̂gk\widehat\tau_k=\sum_g\omega_{gk}\widehat\tau_{gk}, con ωgk=n1gk/∑hn1hk\omega_{gk}=n_{1gk}/\sum_hn_{1hk}. Esta construcción sigue la lógica de comparaciones grupo-tiempo para adopción escalonada (Callaway & Sant’Anna, 2021), aunque la implementación es diagnóstica. La varianza condicional de una media fue

V̂(μ̂)=∑iwi2(ΔYi−μ̂)2(∑iwi)2, \widehat V(\widehat\mu)=\frac{\sum_iw_i^2(\Delta Y_i-\widehat\mu)^2}{(\sum_iw_i)^2},

y los intervalos normales se calcularon como τ̂±1.96SÊ\widehat\tau\pm1.96\widehat{SE}. No incorporan la estimación de los pesos, dependencia espacial ni corrección por pruebas múltiples.

4.4 Balance y soporte efectivo

Para la covariable jj, el balance se midió mediante

SMDj=X‾1jw−X‾0jw(s1j2,w+s0j2,w)/2, SMD_j=\frac{\bar X_{1j}^{w}-\bar X_{0j}^{w}} {\sqrt{(s_{1j}^{2,w}+s_{0j}^{2,w})/2}},

con criterio |SMDj|≤0.10|SMD_j|\leq0.10. El tamaño efectivo fue ESSd=(∑Di=dwi)2/∑Di=dwi2ESS_d=(\sum_{D_i=d}w_i)^2/\sum_{D_i=d}w_i^2.

4.5 Análisis económico

Se procesaron 120 archivos mensuales del portal de la DNIT (Dirección Nacional de Ingresos Tributarios, 2026). De 13.360.960 registros fuente se seleccionaron 42.663 ítems físicos del capítulo 44. Para el año tt, UVt=FOBt/QtUV_t=FOB_t/Q_t, donde QtQ_t son toneladas netas. La participación del grupo cc fue sct=FOBct/FOBts_{ct}=FOB_{ct}/FOB_t y la concentración de destinos se resumió mediante HHI=∑m(100sm)2HHI=\sum_m(100s_m)^2. El cociente FOB/tonelada es un valor unitario agregado, no el precio de un producto homogéneo.

4.6 Reproducibilidad y alcance

El documento comienza en archivos analíticos compactos incluidos junto al .Rmd y vuelve a calcular tablas, estimaciones y figuras. La adquisición de rásteres, las consultas autenticadas a Earth Engine y la descarga mensual de aduanas se conservan en el proyecto maestro. El HTML es estático y autocontenido.

Alcance inferencial. Las estimaciones representan evidencia asociativa fuerte. El diseño reduce desequilibrios observados, pero todavía no justifica una afirmación causal definitiva.

5 Resultados

5.1 Integridad de los datos analíticos

transiciones <- read.csv(data_path("transiciones_uso_suelo.csv"), check.names = FALSE)
panel <- read.csv(gzfile(data_path("panel_anual.csv.gz")), check.names = FALSE)
pesos <- read.csv(data_path("pesos_solapamiento.csv"), check.names = FALSE)
export_m <- read.csv(gzfile(data_path("exportaciones_mensuales.csv.gz")), check.names = FALSE)
bcp <- read.csv(data_path("contexto_bcp_2024.csv"), check.names = FALSE)
balance_publicado <- read.csv(data_path("balance_publicado.csv"), check.names = FALSE)
eventos_publicados <- read.csv(data_path("eventos_publicados.csv"), check.names = FALSE)
mapa_distritos <- read.csv(data_path("mapa_distritos.csv"), check.names = FALSE)
mapa_unidades <- read.csv(gzfile(data_path("mapa_unidades.csv.gz")), check.names = FALSE)
manifest <- read.csv(data_path("MANIFEST_SHA256.csv"), check.names = FALSE)

transiciones$validado_infona <- tolower(as.character(transiciones$validado_infona)) == "true"

stopifnot(
  nrow(panel) == 10660,
  length(unique(panel$id_unidad)) == 1066,
  sum(duplicated(panel[c("id_unidad", "anio")])) == 0,
  sum(transiciones$validado_infona) == 533,
  sum(mapa_unidades$tratado == 1) == 533,
  sum(mapa_unidades$tratado == 0) == 533,
  all(2015:2024 %in% unique(panel$anio)),
  all(2015:2024 %in% unique(export_m$anio))
)

auditoria <- data.frame(
  Control = c(
    "Archivos incluidos en el manifiesto", "Unidades del panel", "Filas anuales",
    "Duplicados unidad-año", "Transiciones validadas", "Controles persistentes",
    "Meses aduaneros representados"
  ),
  Resultado = c(
    nrow(manifest), length(unique(panel$id_unidad)), nrow(panel),
    sum(duplicated(panel[c("id_unidad", "anio")])),
    sum(transiciones$validado_infona), sum(mapa_unidades$tratado == 0),
    length(unique(sprintf("%04d-%02d", export_m$anio, export_m$mes)))
  ),
  Esperado = c(nrow(manifest), 1066, 10660, 0, 533, 533, 120)
)
auditoria$Estado <- ifelse(auditoria$Resultado == auditoria$Esperado, "OK", "REVISAR")
knitr::kable(auditoria, caption = "Controles de integridad ejecutados durante el render.")
Controles de integridad ejecutados durante el render.
Control Resultado Esperado Estado
Archivos incluidos en el manifiesto 9 9 OK
Unidades del panel 1066 1066 OK
Filas anuales 10660 10660 OK
Duplicados unidad-año 0 0 OK
Transiciones validadas 533 533 OK
Controles persistentes 533 533 OK
Meses aduaneros representados 120 120 OK

El paquete contiene 9 archivos de datos con tamaño y SHA-256 registrados. La publicación no incorpora los rásteres originales ni los 120 archivos fuente de aduanas; conserva las unidades analíticas necesarias para recalcular los resultados mostrados.

5.2 Expansión forestal validada

tratamientos <- transiciones[transiciones$validado_infona, ]

cohortes <- aggregate(
  superficie_ha ~ anio_conversion,
  data = tratamientos,
  FUN = sum
)
poligonos <- aggregate(id ~ anio_conversion, data = tratamientos, FUN = length)
cohortes <- merge(cohortes, poligonos, by = "anio_conversion", all = TRUE)
names(cohortes) <- c("Año", "Superficie_ha", "Polígonos")
cohortes$Participación_pct <- 100 * cohortes$Superficie_ha / sum(cohortes$Superficie_ha)
cohortes <- cohortes[order(cohortes$Año), ]

distritos <- aggregate(
  cbind(superficie_ha, conteo = rep(1, nrow(tratamientos))) ~ distrito,
  data = tratamientos,
  FUN = sum
)
names(distritos) <- c("Distrito", "Superficie_ha", "Polígonos")
distritos$Participación_pct <- 100 * distritos$Superficie_ha / sum(distritos$Superficie_ha)
distritos <- distritos[order(-distritos$Superficie_ha), ]

superficie_total <- sum(cohortes$Superficie_ha)
reciente_pct <- 100 * sum(cohortes$Superficie_ha[cohortes$Año >= 2022]) / superficie_total
dos_distritos_pct <- sum(head(distritos$Participación_pct, 2))

tabla_cohortes <- transform(
  cohortes,
  Superficie_ha = fmt_num(Superficie_ha),
  Participación_pct = fmt_pct(Participación_pct, 2)
)
names(tabla_cohortes) <- c("Año", "Superficie (ha)", "Polígonos", "Participación")
knitr::kable(tabla_cohortes, align = c("r", "r", "r", "r"),
             caption = "Transiciones validadas de pastura a plantación forestal por cohorte.")
Transiciones validadas de pastura a plantación forestal por cohorte.
Año Superficie (ha) Polígonos Participación
2017 6,21 1 0,18 %
2018 1,26 1 0,04 %
2019 9,27 4 0,27 %
2020 48,60 10 1,42 %
2021 352,26 35 10,28 %
2022 1.315,71 115 38,40 %
2023 674,10 158 19,68 %
2024 1.018,53 209 29,73 %

Se validaron 533 transiciones, equivalentes a 3.425,94 ha. Las cohortes 2022–2024 concentran 87,81 % del área; Sargento José Félix López y Paso Barreto reúnen 94,42 %.

El máximo anual de superficie ocurrió en 2022, mientras que 2024 registró el mayor número de polígonos. La combinación de concentración espacial y crecimiento del número de unidades muestra que la expansión reciente involucró tanto áreas extensas como la multiplicación de polígonos menores.

par(mar = c(1.5, 1.5, 3.2, 1), family = "sans")
plot(
  NA,
  xlim = range(mapa_distritos$x), ylim = range(mapa_distritos$y), asp = 1,
  axes = FALSE, xlab = "", ylab = "",
  main = "Unidades tratadas y controles de pastura persistente"
)
for (part in unique(mapa_distritos$parte)) {
  z <- mapa_distritos[mapa_distritos$parte == part, ]
  z <- z[order(z$orden), ]
  polygon(z$x, z$y, col = "#F4F6F2", border = "#AAB2B8", lwd = 0.7)
}
controls_map <- mapa_unidades[mapa_unidades$tratado == 0, ]
treated_map <- mapa_unidades[mapa_unidades$tratado == 1, ]
points(controls_map$x, controls_map$y, pch = 16, cex = 0.45, col = alpha_color(color$agua, 0.42))
points(treated_map$x, treated_map$y, pch = 16, cex = 0.55, col = alpha_color(color$tierra, 0.68))
legend(
  "bottomleft", bty = "n", pch = 16,
  col = c(color$tierra, color$agua),
  legend = c("Transición validada", "Pastura persistente"),
  pt.cex = 0.9, cex = 0.88
)
mtext("Fuente: elaboración reproducible con MapBiomas e INFONA", side = 1, line = 0.25, cex = 0.75, col = color$gris)
Figura 1. Unidades tratadas y controles de pastura persistente en Concepción. Las coordenadas del paquete público fueron redondeadas a 100 m para la visualización departamental.

Figura 1. Unidades tratadas y controles de pastura persistente en Concepción. Las coordenadas del paquete público fueron redondeadas a 100 m para la visualización departamental.

5.3 Balance pretratamiento

excluded <- c(
  "id_unidad", "tratado", "distrito", "cohorte_balance",
  "propensity_score", "peso_solapamiento"
)
covariables <- setdiff(names(pesos), excluded)

smd_one <- function(variable) {
  t <- pesos$tratado == 1
  x_t <- pesos[[variable]][t]
  x_c <- pesos[[variable]][!t]
  w_t <- pesos$peso_solapamiento[t]
  w_c <- pesos$peso_solapamiento[!t]
  var_t <- mean((x_t - mean(x_t))^2)
  var_c <- mean((x_c - mean(x_c))^2)
  before <- (mean(x_t) - mean(x_c)) / sqrt((var_t + var_c) / 2)
  after <- (wmean(x_t, w_t) - wmean(x_c, w_c)) /
    sqrt((wvar(x_t, w_t) + wvar(x_c, w_c)) / 2)
  data.frame(variable = variable, smd_sin_ponderar = before, smd_ponderado = after)
}

balance <- do.call(rbind, lapply(covariables, smd_one))
balance$estado <- ifelse(abs(balance$smd_ponderado) <= 0.10, "OK", "REVISAR")

ess <- function(w) sum(w)^2 / sum(w^2)
ess_t <- ess(pesos$peso_solapamiento[pesos$tratado == 1])
ess_c <- ess(pesos$peso_solapamiento[pesos$tratado == 0])

balance_check <- merge(
  balance,
  balance_publicado[c("variable", "smd_sin_ponderar", "smd_ponderado")],
  by = "variable", suffixes = c("_r", "_publicado")
)
max_balance_diff <- max(
  abs(balance_check$smd_sin_ponderar_r - balance_check$smd_sin_ponderar_publicado),
  abs(balance_check$smd_ponderado_r - balance_check$smd_ponderado_publicado)
)
stopifnot(max_balance_diff < 1e-10)

balance_show <- balance[order(-abs(balance$smd_sin_ponderar)), ]
balance_show <- head(balance_show, 10)
balance_show$smd_sin_ponderar <- fmt_num(balance_show$smd_sin_ponderar, 3)
balance_show$smd_ponderado <- fmt_num(balance_show$smd_ponderado, 3)
names(balance_show) <- c("Variable", "SMD sin ponderar", "SMD ponderado", "Estado")
knitr::kable(balance_show, caption = "Diez mayores desequilibrios iniciales y su resultado ponderado.")
Diez mayores desequilibrios iniciales y su resultado ponderado.
Variable SMD sin ponderar SMD ponderado Estado
8 distancia_industria_km 0,864 0,008 OK
7 distancia_agua_m 0,657 0,003 OK
2 ndvi -0,627 -0,006 OK
10 ndvi_tendencia_m3_m1 -0,526 -0,015 OK
3 lst_c -0,516 -0,009 OK
1 superficie_ha 0,374 0,072 OK
11 lst_c_cambio_m2_m1 -0,334 -0,004 OK
5 pet 0,267 0,002 OK
12 lst_c_tendencia_m3_m1 -0,220 -0,003 OK
4 et -0,193 -0,001 OK

La ponderación redujo el máximo |SMD| de 0,864 a 0,072. Las 18 covariables quedaron por debajo de 0,10. El tamaño efectivo fue 229,4 tratamientos y 270,7 controles.

5.4 Respuesta biofísica después de la conversión

analysis <- merge(
  panel,
  pesos[c("id_unidad", "peso_solapamiento")],
  by = "id_unidad", all = FALSE
)
analysis$cohorte_balance <- ifelse(
  analysis$tratado == 1,
  analysis$anio_conversion,
  analysis$cohorte_referencia
)

estimate_group_time <- function(data, outcome, cohort, event_time) {
  target_year <- cohort + event_time
  base_year <- cohort - 1
  if (event_time == -1 || target_year < 2015 || target_year > 2024) return(NULL)
  z <- data[data$cohorte_balance == cohort & data$anio %in% c(base_year, target_year),
            c("id_unidad", "tratado", "anio", "peso_solapamiento", outcome)]
  base <- z[z$anio == base_year, c("id_unidad", "tratado", "peso_solapamiento", outcome)]
  target <- z[z$anio == target_year, c("id_unidad", "tratado", "peso_solapamiento", outcome)]
  names(base)[4] <- "base"
  names(target)[4] <- "target"
  wide <- merge(base, target, by = c("id_unidad", "tratado", "peso_solapamiento"))
  wide <- wide[complete.cases(wide), ]
  if (nrow(wide) == 0) return(NULL)

  group_stats <- lapply(c(1, 0), function(group) {
    g <- wide[wide$tratado == group, ]
    if (nrow(g) == 0) return(NULL)
    change <- g$target - g$base
    w <- g$peso_solapamiento
    m <- wmean(change, w)
    variance <- sum(w^2 * (change - m)^2) / sum(w)^2
    list(mean = m, variance = variance, n = nrow(g), sum_w = sum(w))
  })
  if (is.null(group_stats[[1]]) || is.null(group_stats[[2]])) return(NULL)
  estimate <- group_stats[[1]]$mean - group_stats[[2]]$mean
  se <- sqrt(group_stats[[1]]$variance + group_stats[[2]]$variance)
  data.frame(
    variable = outcome, cohorte = cohort, event_time = event_time,
    anio_resultado = target_year, estimacion = estimate, error_estandar = se,
    ic95_inferior = estimate - 1.96 * se, ic95_superior = estimate + 1.96 * se,
    n_tratados = group_stats[[1]]$n, n_controles = group_stats[[2]]$n,
    suma_pesos_tratados = group_stats[[1]]$sum_w,
    suma_pesos_controles = group_stats[[2]]$sum_w
  )
}

outcomes <- c("ndvi", "evi", "ndmi", "lst_c", "et", "pet", "precipitacion_mm")
event_grid <- expand.grid(
  variable = outcomes,
  cohorte = 2020:2024,
  event_time = -4:4,
  stringsAsFactors = FALSE
)
event_grid <- event_grid[event_grid$event_time != -1, ]
group_time_list <- lapply(seq_len(nrow(event_grid)), function(i) {
  estimate_group_time(
    analysis,
    event_grid$variable[i], event_grid$cohorte[i], event_grid$event_time[i]
  )
})
group_time <- do.call(rbind, group_time_list[!vapply(group_time_list, is.null, logical(1))])

aggregate_events <- function(z) {
  total <- sum(z$n_tratados)
  if (total < 30) return(NULL)
  w <- z$n_tratados / total
  estimate <- sum(w * z$estimacion)
  se <- sqrt(sum((w * z$error_estandar)^2))
  data.frame(
    variable = z$variable[1], event_time = z$event_time[1],
    estimacion = estimate, error_estandar = se,
    ic95_inferior = estimate - 1.96 * se,
    ic95_superior = estimate + 1.96 * se,
    cohortes_con_soporte = length(unique(z$cohorte)),
    n_tratados_suma_cohortes = total
  )
}
event_split <- split(group_time, interaction(group_time$variable, group_time$event_time, drop = TRUE))
eventos <- do.call(rbind, lapply(event_split, aggregate_events))
eventos <- eventos[order(eventos$variable, eventos$event_time), ]
row.names(eventos) <- NULL

event_check <- merge(
  eventos,
  eventos_publicados,
  by = c("variable", "event_time"), suffixes = c("_r", "_publicado")
)
max_event_diff <- max(abs(event_check$estimacion_r - event_check$estimacion_publicado))
stopifnot(max_event_diff < 1e-10)

La reestimación R reproduce los coeficientes publicados con una diferencia máxima de 7.8e-14. Esto verifica el cálculo desde el panel y los pesos, no una copia de la tabla final.

event_table <- eventos[eventos$event_time %in% 0:3 & eventos$variable %in% c("ndvi", "evi", "ndmi", "lst_c", "et"), ]
event_wide <- reshape(
  event_table[c("variable", "event_time", "estimacion")],
  idvar = "event_time", timevar = "variable", direction = "wide"
)
support <- eventos[eventos$variable == "ndvi" & eventos$event_time %in% 0:3,
                   c("event_time", "n_tratados_suma_cohortes")]
event_wide <- merge(event_wide, support, by = "event_time")
event_wide <- event_wide[order(event_wide$event_time), ]
event_wide[-c(1, ncol(event_wide))] <- lapply(event_wide[-c(1, ncol(event_wide))], fmt_num, digits = 3)
names(event_wide) <- c("Tiempo", "ET", "EVI", "LST", "NDMI", "NDVI", "Tratados con soporte")
knitr::kable(event_wide, caption = "Diferencias ponderadas de cambios respecto de t = -1.")
Diferencias ponderadas de cambios respecto de t = -1.
Tiempo ET EVI LST NDMI NDVI Tratados con soporte
0 243,183 0,119 -2,212 0,191 0,168 527
1 341,743 0,157 -2,723 0,273 0,239 318
2 257,571 0,116 -2,405 0,244 0,210 160
3 198,845 0,084 -2,394 0,176 0,163 45

5.4.1 Explorador dinámico de variables

El gráfico siguiente permite recorrer la trayectoria estimada desde cuatro años antes hasta tres años después de la conversión. Seleccione una variable, mueva el control temporal o utilice la reproducción automática. La línea se revela progresivamente para mostrar cómo cambia la diferencia entre plantaciones y pasturas persistentes.

FIGURA INTERACTIVA

Trayectoria biofísica alrededor de la conversión

Estimaciones ponderadas e intervalos de confianza del 95 %. t = -1 es el periodo de referencia.

Datos recalculados en R
Evolución de la variable seleccionada Gráfico interactivo del estudio de eventos.
MomentoAño de conversión
Estimación—
Intervalo 95 %—
Soporte—

Diferencia estimada Intervalo 95 % Sin diferencia 
plot_event <- function(variable, title, ylab, col) {
  z <- eventos[eventos$variable == variable, ]
  z <- rbind(
    z[, c("event_time", "estimacion", "ic95_inferior", "ic95_superior")],
    data.frame(event_time = -1, estimacion = 0, ic95_inferior = 0, ic95_superior = 0)
  )
  z <- z[order(z$event_time), ]
  ylim <- range(z$ic95_inferior, z$ic95_superior)
  plot(
    z$event_time, z$estimacion, type = "o", pch = 16, lwd = 2,
    col = col, xlab = "Años desde la conversión", ylab = ylab,
    main = title, xaxt = "n", ylim = ylim
  )
  axis(1, at = -4:3)
  abline(h = 0, col = color$gris, lwd = 1)
  abline(v = -0.5, col = color$gris, lty = 2)
  has_interval <- z$ic95_superior > z$ic95_inferior
  arrows(z$event_time[has_interval], z$ic95_inferior[has_interval],
         z$event_time[has_interval], z$ic95_superior[has_interval],
         angle = 90, code = 3, length = 0.035, col = col)
  grid(nx = NA, ny = NULL, col = color$claro)
  box(bty = "l")
}
par(mfrow = c(2, 2), mar = c(4.2, 4.6, 3, 1), oma = c(0, 0, 2, 0), family = "sans")
plot_event("ndvi", "A  NDVI", "Diferencia en el cambio", color$verde)
plot_event("lst_c", "B  Temperatura superficial", "Diferencia (°C)", color$rojo)
plot_event("et", "C  Evapotranspiración", "Diferencia (mm/año)", color$agua)
plot_event("pet", "D  Evapotranspiración potencial", "Diferencia (mm/año)", color$naranja)
mtext("Plantaciones frente a pasturas persistentes", outer = TRUE, cex = 1.25, font = 2)
Figura 2. Estudio de eventos ponderado. Los puntos son diferencias de cambios respecto de t = -1; las barras representan intervalos normales de 95 %.

Figura 2. Estudio de eventos ponderado. Los puntos son diferencias de cambios respecto de t = -1; las barras representan intervalos normales de 95 %.

Resultado central. En el año de conversión, el NDVI aumentó 0,168, el EVI 0,119 y el NDMI 0,191. La temperatura superficial disminuyó 2,21 °C y la evapotranspiración aumentó 243,18 mm/año. Un año después, las diferencias alcanzaron -2,72 °C y 341,74 mm/año.

En t = 2 y t = 3 la dirección se mantuvo, aunque el soporte disminuyó de 527 tratamientos en el año de conversión a 45 en el tercer año. Solo una de 21 comparaciones agregadas previas fue significativa al 5 %: NDVI en t = -4. Las diferencias previas restantes fueron mayormente compatibles con cero. La precipitación difirió en t = 2 y t = 3, por lo que puede confundir la respuesta tardía de ET y de los índices espectrales.

5.5 Exportaciones forestales y contexto productivo

annual <- aggregate(
  cbind(registros, kg_neto, fob_usd) ~ anio,
  data = export_m,
  FUN = sum
)
annual$toneladas <- annual$kg_neto / 1000
annual$usd_tonelada <- annual$fob_usd / annual$toneladas

processed <- aggregate(
  fob_usd ~ anio,
  data = export_m[export_m$grado_procesamiento %in% c("Semiprocesado", "Elaborado"), ],
  FUN = sum
)
names(processed)[2] <- "fob_procesado"
annual <- merge(annual, processed, by = "anio", all.x = TRUE)
annual$participacion_procesados_pct <- 100 * annual$fob_procesado / annual$fob_usd
annual <- annual[order(annual$anio), ]

composition <- aggregate(fob_usd ~ anio + grado_procesamiento, data = export_m, FUN = sum)
year_totals <- aggregate(fob_usd ~ anio, data = composition, FUN = sum)
names(year_totals)[2] <- "total"
composition <- merge(composition, year_totals, by = "anio")
composition$participacion <- 100 * composition$fob_usd / composition$total

destinations <- aggregate(fob_usd ~ pais_destino, data = export_m, FUN = sum)
destinations <- destinations[order(-destinations$fob_usd), ]
destinations$share <- 100 * destinations$fob_usd / sum(destinations$fob_usd)
hhi <- sum(destinations$share^2)

fob_growth <- 100 * (annual$fob_usd[annual$anio == 2024] / annual$fob_usd[annual$anio == 2015] - 1)
volume_growth <- 100 * (annual$toneladas[annual$anio == 2024] / annual$toneladas[annual$anio == 2015] - 1)
unit_growth <- 100 * (annual$usd_tonelada[annual$anio == 2024] / annual$usd_tonelada[annual$anio == 2015] - 1)
peak_year <- annual$anio[which.max(annual$fob_usd)]

selected_years <- annual[annual$anio %in% c(2015, 2018, 2020, 2022, 2024), ]
export_table <- data.frame(
  Año = selected_years$anio,
  `FOB (millones USD)` = fmt_num(selected_years$fob_usd / 1e6),
  `Volumen (miles t)` = fmt_num(selected_years$toneladas / 1e3),
  `USD/t` = fmt_num(selected_years$usd_tonelada),
  `Procesados (% FOB)` = fmt_num(selected_years$participacion_procesados_pct)
)
knitr::kable(export_table, caption = "Indicadores seleccionados de exportaciones paraguayas del capítulo 44.")
Indicadores seleccionados de exportaciones paraguayas del capítulo 44.
Año FOB..millones.USD. Volumen..miles.t. USD.t Procesados….FOB.
2015 70,46 148,88 473,27 36,03
2018 71,96 149,49 481,39 25,23
2020 57,31 126,13 454,38 27,73
2022 105,71 233,52 452,67 31,60
2024 96,31 195,45 492,76 31,51

Entre 2015 y 2024, el FOB aumentó 36,7 %, el volumen 31,3 % y el valor unitario 4,1 %. El máximo ocurrió en 2022. La expansión fue predominantemente volumétrica.

En el periodo completo se exportaron 1,618 millones de toneladas por USD 772,42 millones. El índice HHI de destinos fue 844, compatible con baja concentración agregada. Estos datos son nacionales y no contienen departamento ni especie de origen.

par(mfrow = c(2, 2), mar = c(4.2, 4.5, 3, 1), family = "sans")
barplot(annual$fob_usd / 1e6, names.arg = annual$anio, col = color$bosque,
        ylab = "Millones de USD", main = "A  Valor FOB", las = 2, cex.names = 0.75)
barplot(annual$toneladas / 1e3, names.arg = annual$anio, col = color$agua,
        ylab = "Miles de toneladas", main = "B  Volumen", las = 2, cex.names = 0.75)
plot(annual$anio, annual$usd_tonelada, type = "o", pch = 16, lwd = 2,
     col = color$tierra, xlab = "Año", ylab = "USD por tonelada", main = "C  Valor unitario")
grid(nx = NA, ny = NULL, col = color$claro)

grades <- c("Bioenergia", "Primario", "Semiprocesado", "Elaborado")
comp_matrix <- sapply(annual$anio, function(y) {
  z <- composition[composition$anio == y, ]
  setNames(z$participacion, z$grado_procesamiento)[grades]
})
comp_matrix[is.na(comp_matrix)] <- 0
barplot(comp_matrix, names.arg = annual$anio, col = c("#5B8E7D", "#D6A85F", "#7EA6C4", "#A46B9B"),
        ylab = "% del FOB", main = "D  Composición", las = 2, cex.names = 0.75)
legend("top", legend = grades, fill = c("#5B8E7D", "#D6A85F", "#7EA6C4", "#A46B9B"),
       bty = "n", cex = 0.68, ncol = 2)
Figura 3. Exportaciones paraguayas del capítulo 44, 2015–2024. Los datos son nacionales y no identifican el departamento ni la especie de origen.

Figura 3. Exportaciones paraguayas del capítulo 44, 2015–2024. Los datos son nacionales y no identifican el departamento ni la especie de origen.

5.6 Síntesis integrada

par(mfrow = c(2, 2), mar = c(4.2, 4.5, 3.1, 4.2), family = "sans")

bp <- barplot(cohortes$Superficie_ha, names.arg = cohortes$Año, col = color$bosque,
              ylab = "Superficie (ha)", main = "A  Expansión validada", las = 2, cex.names = 0.78)
par(new = TRUE)
plot(bp, cohortes$Polígonos, type = "o", pch = 16, lwd = 2, col = color$tierra,
     axes = FALSE, xlab = "", ylab = "", ylim = c(0, max(cohortes$Polígonos) * 1.15))
axis(4, col.axis = color$tierra)
mtext("Polígonos", side = 4, line = 2.5, col = color$tierra)

ndvi_event <- eventos[eventos$variable == "ndvi", ]
ndvi_event <- rbind(ndvi_event, data.frame(
  variable = "ndvi", event_time = -1, estimacion = 0, error_estandar = 0,
  ic95_inferior = 0, ic95_superior = 0, cohortes_con_soporte = NA,
  n_tratados_suma_cohortes = NA
))
ndvi_event <- ndvi_event[order(ndvi_event$event_time), ]
plot(ndvi_event$event_time, ndvi_event$estimacion, type = "o", pch = 16, lwd = 2,
     col = color$verde, xlab = "Años desde la conversión", ylab = "Diferencia en cambio",
     main = "B  Respuesta NDVI", xaxt = "n")
axis(1, at = -4:3)
abline(h = 0, col = color$gris)
ndvi_interval <- ndvi_event$ic95_superior > ndvi_event$ic95_inferior
arrows(ndvi_event$event_time[ndvi_interval], ndvi_event$ic95_inferior[ndvi_interval],
       ndvi_event$event_time[ndvi_interval], ndvi_event$ic95_superior[ndvi_interval],
       angle = 90, code = 3, length = 0.03, col = color$verde)

bp2 <- barplot(annual$fob_usd / 1e6, names.arg = annual$anio, col = alpha_color(color$bosque, 0.75),
               ylab = "FOB (millones USD)", main = "C  Exportaciones nacionales", las = 2, cex.names = 0.75)
par(new = TRUE)
plot(bp2, annual$toneladas / 1e3, type = "o", pch = 16, lwd = 2, col = color$agua,
     axes = FALSE, xlab = "", ylab = "", ylim = c(0, max(annual$toneladas / 1e3) * 1.15))
axis(4, col.axis = color$agua)
mtext("Miles de toneladas", side = 4, line = 2.5, col = color$agua)

barplot(comp_matrix, names.arg = annual$anio,
        col = c("#5B8E7D", "#D6A85F", "#7EA6C4", "#A46B9B"),
        ylab = "% del FOB", main = "D  Composición del FOB", las = 2, cex.names = 0.75)
legend("top", legend = grades, fill = c("#5B8E7D", "#D6A85F", "#7EA6C4", "#A46B9B"),
       bty = "n", cex = 0.65, ncol = 2)
Figura 4. Síntesis reproducida en R: expansión validada, respuesta NDVI, exportaciones nacionales y composición del FOB.

Figura 4. Síntesis reproducida en R: expansión validada, respuesta NDVI, exportaciones nacionales y composición del FOB.

5.7 Contexto regional del BCP

bcp_selected <- bcp[
  bcp$actividad %in% c("PIB regional", "Ganaderia, Forestal, Pesca y Mineria", "Manufactura") &
    bcp$tipo_indicador %in% c("crecimiento_pib", "crecimiento_vab", "estructura_pib", "contribucion_crecimiento"),
  c("tipo_indicador", "actividad", "valor", "unidad")
]
names(bcp_selected) <- c("Indicador", "Actividad", "Valor", "Unidad")
bcp_selected$Valor <- fmt_num(bcp_selected$Valor, 1)
knitr::kable(bcp_selected, caption = "Contexto económico de Concepción en 2024 según el BCP.")
Contexto económico de Concepción en 2024 según el BCP.
Indicador Actividad Valor Unidad
1 crecimiento_pib PIB regional 7,3 porcentaje_interanual
3 crecimiento_vab Ganaderia, Forestal, Pesca y Mineria 18,3 porcentaje_interanual
4 crecimiento_vab Manufactura 6,9 porcentaje_interanual
10 contribucion_crecimiento Ganaderia, Forestal, Pesca y Mineria 2,0 puntos_porcentuales
11 contribucion_crecimiento Manufactura 2,0 puntos_porcentuales
14 estructura_pib Ganaderia, Forestal, Pesca y Mineria 12,1 porcentaje
15 estructura_pib Manufactura 29,2 porcentaje

El agregado “Ganadería, Forestal, Pesca y Minería” no permite aislar el aporte forestal. De forma análoga, las exportaciones del capítulo 44 no pueden atribuirse a Concepción ni a los polígonos estudiados.

6 Discusión

6.1 Una firma biofísica convergente

La coincidencia de indicadores independientes constituye el hallazgo principal. NDVI y EVI aumentan con el vigor y densidad de la cubierta; NDMI indica una señal más húmeda del dosel; la reducción de la temperatura superficial muestra mayor enfriamiento de la superficie; y la ET confirma un aumento del flujo de agua hacia la atmósfera. La convergencia reduce la posibilidad de que el patrón dependa de una sola métrica o colección satelital.

La mayor respuesta en t = 1 es compatible con el establecimiento inicial del dosel, aunque la etiqueta anual puede integrar fechas de implantación diferentes. La resolución de MODIS también produce mezcla entre el polígono y su entorno. Los resultados describen una modificación del balance superficial, no una mejora hídrica integral: mayor ET puede coexistir con menor disponibilidad de agua aguas abajo (Farley et al., 2005; Jackson et al., 2005).

6.2 Alcance de la comparación contrafactual

La reducción del máximo |SMD| de 0,864 a 0,072 muestra que la ponderación hizo comparables los grupos en las covariables observadas. La persistencia de la pastura y la separación espacial disminuyen la contaminación directa, y los placebos previos son mayormente compatibles con trayectorias paralelas. Persisten factores no observados, dependencia espacial e incertidumbre del primer paso. Por ello, la evidencia sustenta una asociación fuerte y temporalmente alineada, no una prueba causal definitiva.

6.3 Territorio y cadena de valor

La expansión validada en Concepción y el máximo nacional de exportaciones se concentraron alrededor de 2022. La coincidencia es compatible con un entorno productivo favorable, pero no identifica un vínculo de origen y destino. La alta participación de bioenergía y el crecimiento principalmente volumétrico sugieren que el aumento exportador no estuvo acompañado por una transformación equivalente hacia productos de mayor elaboración. La interpretación debe complementarse con precios específicos, costos logísticos y trazabilidad territorial.

6.4 Implicancias para el monitoreo

Sargento José Félix López y Paso Barreto son áreas prioritarias para mediciones de campo. Una red de seguimiento debería combinar edad y especie del rodal, precipitación local, humedad del suelo, nivel freático, caudal, propiedades edáficas y manejo. La teledetección aporta cobertura espacial; las mediciones terrestres son necesarias para interpretar mecanismos y validar escalas.

7 Limitaciones

  1. Las fuentes verifican plantaciones forestales, pero no demuestran que todos los polígonos sean eucalipto.
  2. Las resoluciones espaciales varían entre 10 m y aproximadamente 5 km.
  3. Los intervalos no incorporan la estimación de los pesos, agrupamiento espacial ni multiplicidad.
  4. La precipitación difiere entre grupos en t = 2 y t = 3.
  5. El soporte disminuye hasta 45 tratamientos en t = 3.
  6. ET y NDMI no sustituyen mediciones de caudal, acuífero o humedad profunda.
  7. Las exportaciones nacionales no pueden atribuirse a Concepción ni a las unidades estudiadas.
  8. La categoría regional del BCP combina varias actividades y no identifica la silvicultura por separado.

8 Conclusiones

La expansión de plantaciones forestales sobre pasturas en Concepción fue reciente, espacialmente concentrada y mensurable. Se validaron 533 transiciones por 3.425,94 ha; casi nueve décimas partes del área corresponden a 2022–2024 y dos distritos reúnen 94,42 %.

La hipótesis biofísica recibió apoyo consistente. Frente a pasturas persistentes ponderadas, la transición estuvo asociada con aumentos de NDVI, EVI y NDMI, una reducción aproximada de 2–3 °C en temperatura superficial y un aumento de 200–340 mm/año en ET. La hipótesis productiva también fue respaldada en su dimensión descriptiva: el crecimiento exportador nacional de 2015–2024 estuvo impulsado principalmente por volumen. No existe evidencia para atribuir ese crecimiento a Concepción ni a las plantaciones detectadas.

El estudio demuestra una asociación espacio-temporal fuerte y establece una línea de base reproducible. Una afirmación causal requiere métodos formales para adopción escalonada, errores agrupados o bootstrap espacial, ajuste climático explícito y validación de campo.

Conclusión en una frase. Después de la pastura, la nueva cobertura forestal deja una huella observable: más verdor y humedad del dosel, una superficie más fría y un mayor flujo de agua hacia la atmósfera; el desafío es medir cómo esa huella se traduce en disponibilidad hídrica y desarrollo productivo territorial.

9 Declaraciones

9.1 Contribución del autor

Diego Bernardo Meza: conceptualización, metodología, curación de datos, análisis formal, programación, visualización y redacción.

9.2 Consideraciones éticas

El estudio utiliza cartografía, imágenes satelitales y registros económicos secundarios. No incluye participantes humanos, animales de experimentación ni datos personales.

9.3 Financiamiento y declaración de intereses

No se consigna una subvención específica en el proyecto reproducible. El autor declara su afiliación con Paracel S.A.; esta relación institucional se informa para que sea considerada al interpretar el trabajo. Los resultados se derivan de fuentes públicas y procedimientos documentados. Cualquier interés adicional deberá declararse antes de un envío editorial.

9.4 Agradecimientos

Se reconoce al Banco Central del Paraguay, la Dirección Nacional de Ingresos Tributarios, el Instituto Forestal Nacional, MapBiomas Paraguay y las agencias responsables de las colecciones satelitales por mantener fuentes que permiten la investigación reproducible.

10 Reproducibilidad y disponibilidad

Los datos analíticos compactos, el código R y los manifiestos SHA-256 acompañan al documento. Las tablas y figuras se reconstruyen durante el render. Los rásteres originales, las consultas autenticadas y los archivos mensuales de aduanas permanecen en el proyecto maestro. El HTML final es autocontenido y no requiere rutas locales para visualizarse.

Ver información de la sesión R 
sessionInfo()
## R version 4.4.3 (2025-02-28 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19044)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8    LC_MONETARY=Spanish_Spain.utf8
## [4] LC_NUMERIC=C                   LC_TIME=Spanish_Spain.utf8    
## 
## time zone: Etc/GMT+4
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.35     R6_2.6.1          fastmap_1.2.0     xfun_0.52         cachem_1.1.0     
##  [6] knitr_1.50        htmltools_0.5.8.1 rmarkdown_2.29    lifecycle_1.0.4   cli_3.6.2        
## [11] sass_0.4.10       jquerylib_0.1.4   compiler_4.4.3    rstudioapi_0.16.0 tools_4.4.3      
## [16] evaluate_1.0.4    bslib_0.9.0       yaml_2.3.8        rlang_1.1.6       jsonlite_2.0.0

Estado: versión científica completa preparada y verificada localmente. No publicada. Antes de una eventual presentación deben confirmarse la categoría temática, las reglas editoriales, la declaración de financiamiento y cualquier dato de correspondencia requerido.

11 Referencias

Banco Central del Paraguay. (2026). Cuentas regionales anuales 2024. Banco Central del Paraguay. https://www.bcp.gov.py/documents/20117/1935346/Informe%2BCRA_2024.pdf
Callaway, B., & Sant’Anna, P. H. C. (2021). Difference-in-differences with multiple time periods. Journal of Econometrics, 225(2), 200-230. https://doi.org/10.1016/j.jeconom.2020.12.001
Dirección Nacional de Ingresos Tributarios. (2026). Portal de datos abiertos de comercio exterior. https://datosabiertos.aduana.gov.py/ddaa/app/index.html
Farley, K. A., Jobbágy, E. G., & Jackson, R. B. (2005). Effects of afforestation on water yield: A global synthesis with implications for policy. Global Change Biology, 11(10), 1565-1576. https://doi.org/10.1111/j.1365-2486.2005.01011.x
Funk, C., Peterson, P., Landsfeld, M., Pedreros, D., Verdin, J., Shukla, S., Husak, G., Rowland, J., Harrison, L., Hoell, A., & Michaelsen, J. (2015). The climate hazards infrared precipitation with stations—A new environmental record for monitoring extremes. Scientific Data, 2, 150066. https://doi.org/10.1038/sdata.2015.66
Gao, B. (1996). NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sensing of Environment, 58(3), 257-266. https://doi.org/10.1016/S0034-4257(96)00067-3
Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., & Moore, R. (2017). Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sensing of Environment, 202, 18-27. https://doi.org/10.1016/j.rse.2017.06.031
Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., & Ferreira, L. G. (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83(1–2), 195-213. https://doi.org/10.1016/S0034-4257(02)00096-2
Instituto Forestal Nacional. (2026). Visor geoespacial institucional. https://visor.infona.gov.py
Jackson, R. B., Jobbágy, E. G., Avissar, R., Roy, S. B., Barrett, D. J., Cook, C. W., Farley, K. A., Maitre, D. C. le, McCarl, B. A., & Murray, B. C. (2005). Trading water for carbon with biological carbon sequestration. Science, 310(5756), 1944-1947. https://doi.org/10.1126/science.1119282
Li, F., Morgan, K. L., & Zaslavsky, A. M. (2018). Balancing covariates via propensity score weighting. Journal of the American Statistical Association, 113(521), 390-400. https://doi.org/10.1080/01621459.2016.1260466
MapBiomas Paraguay. (2026). Colección 3 de mapas anuales de cobertura y uso del suelo. https://paraguay.mapbiomas.org/herramientas/
Mu, Q., Zhao, M., & Running, S. W. (2011). Improvements to a MODIS global terrestrial evapotranspiration algorithm. Remote Sensing of Environment, 115(8), 1781-1800. https://doi.org/10.1016/j.rse.2011.02.019
Pekel, J.-F., Cottam, A., Gorelick, N., & Belward, A. S. (2016). High-resolution mapping of global surface water and its long-term changes. Nature, 540(7633), 418-422. https://doi.org/10.1038/nature20584
Tucker, C. J. (1979). Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing of Environment, 8(2), 127-150. https://doi.org/10.1016/0034-4257(79)90013-0