Cómo hacer un mapa de producción agrícola departamento del Cesar
Sthelyn Catalina Medina Trillos
11.11.2026
1. Introducción
Este cuaderno muestra cómo hacer un mapa temático de la producción de un grupo importante de cultivos en un departamento. Se usan dos fuentes principales que se obtuvieron previamente: (i) un archivo con la producción agrícola por municipio (EVA 2025, UPRA); y (ii) un geopackage con los municipios del departamento. Además, se usa un conjunto de datos de vías.
En este caso trabajo con el Departamento del Cesar (código DANE 20), que tiene 25 municipios.
2. Configuración
Primero, cargamos las librerías de R necesarias.
# ejecutar la siguiente línea desde la consola si no tienes los paquetes:
# install.packages(c("tmap", "sf", "dplyr", "readr", "ggplot2", "knitr"))
library(sf)
library(dplyr)
library(readr)
library(tmap)
library(ggplot2)
library(knitr)3. Lectura de los archivos de municipios y cultivos
Voy a trabajar con el Departamento del Cesar. Como mis archivos están en mi portatil y otra carpeta, cambié las rutas, los nombres de archivos y los nombres de las variables según corresponde.
Antes de empezar, guardé los siguientes archivos en mi carpeta
PROYECTO 7:
- un geopackage con los municipios del Cesar
(
Mun_Cesar.gpkg) - un archivo con las estadísticas EVA 2025 del Cesar
(
SOLO_CESAR_2025.csv) - un geopackage con las vías del Cesar
(
Vias_Cesar.gpkg)
Definamos las rutas y empecemos listando los archivos disponibles:
ruta <- "C:/PROYECTO 7"
archivo_eva <- file.path(ruta, "SOLO_CESAR_2025.csv")
archivo_munic <- file.path(ruta, "Mun_Cesar.gpkg")
archivo_roads <- file.path(ruta, "Vias_Cesar.gpkg")
list.files(ruta)## [1] "dinamicas_pro_cesar.Rmd"
## [2] "dinamicas_produccion_de_cultivos_cesar.html"
## [3] "dinamicas_produccion_de_cultivos_cesar.Rmd"
## [4] "EVA_Cesar"
## [5] "EVA_Cesar.zip"
## [6] "Evaluaciones_Agropecuarias_Municipales_EVA_20261110.csv"
## [7] "mapa_produccion_agricolas_cesar.html"
## [8] "mapa_produccion_agricolas_cesar.Rmd"
## [9] "Mun_Cesar.gpkg"
## [10] "palma_aceite.jpg"
## [11] "rsconnect"
## [12] "SOLO_CESAR_2025.csv"
## [13] "Vias_Cesar.gpkg"Ahora leamos el archivo de cultivos y revisemos sus columnas:
# el csv puede venir separado por ";" (Excel en español) o por ","
sep <- ifelse(grepl(";", readLines(archivo_eva, n = 1, warn = FALSE)), ";", ",")
eva <- read_delim(archivo_eva, delim = sep,
locale = locale(encoding = "latin1"),
col_types = cols(.default = "c"), show_col_types = FALSE)
names(eva)## [1] "COD_DPTO" "DPTO" "COD_MUN" "MUN" "CULTIVO"
## [6] "CICLO" "GRUPO" "YEAR" "PERIODO" "SIEMBRA_HA"
## [11] "COSECHAS_HA" "PROD_T" "RDTO_HA"Todas las columnas se leen como texto para evitar problemas con los
decimales (algunos valores vienen con punto y otros con coma). Por eso
convierto PROD_T a numérico. Luego selecciono el grupo de
cultivos con mayor producción, Oleaginosas
(principalmente palma de aceite), y me quedo solo con las columnas que
necesito:
(grupo <- eva %>%
mutate(PROD_T = as.numeric(gsub(",", ".", PROD_T))) %>%
filter(trimws(GRUPO) == "Oleaginosas") %>%
transmute(COD_MUN = trimws(COD_MUN), MUNICIPIO = MUN, GRUPO, CULTIVO,
max_prod = PROD_T))Observa que COD_MUN quedó como texto, igual que el
código de los municipios del mapa.
Ahora leamos los municipios que nos interesan:
municipios <- st_read(archivo_munic)## Multiple layers are present in data source C:\PROYECTO 7\Mun_Cesar.gpkg, reading layer `Municipios'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.## Reading layer `Municipios' from data source `C:\PROYECTO 7\Mun_Cesar.gpkg' using driver `GPKG'
## Simple feature collection with 24 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.15633 ymin: 7.683902 xmax: -72.88663 ymax: 10.8524
## Geodetic CRS: WGS 84# paso los nombres de las columnas a minúscula para que coincidan con el resto del cuaderno
geom_col <- attr(municipios, "sf_column")
municipios <- municipios %>% rename_with(tolower, .cols = -all_of(geom_col))
names(municipios)## [1] "id_0" "iso" "name_0" "id_1" "name_1" "id_2"
## [7] "name_2" "type_2" "engtype_2" "nl_name_2" "varname_2" "geom"La capa de municipios no trae el código DANE, así que lo busco a partir del nombre del municipio, usando la tabla EVA, que sí trae los códigos:
# función para comparar nombres sin tildes ni mayúsculas
quita_tildes <- function(x) toupper(trimws(chartr("áéíóúüñÁÉÍÓÚÜÑ", "aeiouunAEIOUUN", x)))
# esta capa trae el nombre del municipio en la columna name_2
if (!"name_2" %in% names(municipios) && !"mpio_cnmbr" %in% names(municipios)) {
stop("No encuentro la columna con el nombre del municipio, revisa names(municipios)")
}
if (!"mpio_cnmbr" %in% names(municipios)) {
municipios$mpio_cnmbr <- municipios$name_2
}
# tabla con los códigos DANE de los municipios (sale del archivo EVA)
codigos <- eva %>%
distinct(COD_MUN, MUN) %>%
mutate(COD_MUN = trimws(COD_MUN), llave = quita_tildes(MUN))
# la capa no trae el código DANE: lo busco por el nombre del municipio
if (!"mpio_cdpmp" %in% names(municipios)) {
llave_munic <- quita_tildes(municipios$mpio_cnmbr)
municipios$mpio_cdpmp <- codigos$COD_MUN[match(llave_munic, codigos$llave)]
# para los nombres que no coinciden exacto (ej. "Codazzi" vs "Agustín Codazzi")
for (i in which(is.na(municipios$mpio_cdpmp))) {
hit <- which(grepl(llave_munic[i], codigos$llave, fixed = TRUE) |
mapply(grepl, codigos$llave, llave_munic[i], MoreArgs = list(fixed = TRUE)))
if (length(hit) == 1) municipios$mpio_cdpmp[i] <- codigos$COD_MUN[hit]
}
}
municipios$mpio_cdpmp <- as.character(municipios$mpio_cdpmp)
municipiosclass(municipios$mpio_cdpmp)## [1] "character"Observa que mpio_cdpmp es de tipo texto (character).
Revisemos si algún municipio quedó sin código y si falta alguno del
archivo EVA en la capa:
# municipios de la capa a los que no les encontré código (debería salir vacío)
municipios %>% st_drop_geometry() %>% filter(is.na(mpio_cdpmp)) %>% select(mpio_cnmbr)# municipios del archivo EVA que no están en la capa
codigos %>% filter(!COD_MUN %in% municipios$mpio_cdpmp) %>% select(COD_MUN, MUN)(cesar.crs <- sf::st_crs(municipios))## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]Respuesta: la capa de municipios del Cesar tiene 24 municipios. Después de filtrar el grupo de cultivos, quedaron 11 municipios con datos, es decir, 13 municipios del Cesar no tienen registro para este grupo en el archivo.
4. Lectura de las vías del Cesar
Para las vías ya tengo un geopackage recortado al departamento del Cesar que hicimos en las primeras clases con el profe
(sf.roads = st_read(archivo_roads))## Reading layer `vias' from data source `C:\PROYECTO 7\Vias_Cesar.gpkg' using driver `GPKG'
## Simple feature collection with 160 features and 5 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: -74.10105 ymin: 7.684987 xmax: -73.02679 ymax: 10.81996
## Geodetic CRS: WGS 84Observa que sf.roads es una colección de simple
features. Ahora respondo las preguntas:
- ¿Cuál es el tipo de geometría de los elementos? MULTILINESTRING.
- ¿Cuántos elementos hay en la colección? 160 elementos.
Ahora, asegurémonos de que las vías y los municipios tengan el mismo sistema de referencia (CRS):
# tomar nota de la sintaxis de esta función
cesar.roads <- st_transform(sf.roads, crs = cesar.crs)
head(cesar.roads)# el archivo ya viene recortado al Cesar, así que uso todas las vías en el mapa
cesar.roads.main <- cesar.roadsEl Cesar tiene 160 elementos de vías.
5. Unión de los datos de producción con los municipios
Empecemos agregando a la tabla de cultivos un atributo nuevo con el código del municipio como texto:
grupo$mpio_cdpmp = as.character(grupo$COD_MUN)Ahora unamos los datos de producción con nuestros municipios:
(munic_grupo = left_join(municipios, grupo, by = "mpio_cdpmp"))Seleccionemos algunos atributos para mostrar después en el mapa:
(munic_grupo2 <- munic_grupo %>%
dplyr::select(mpio_cdpmp, mpio_cnmbr, max_prod))Antes de hacer el mapa, exploremos los números para saber qué esperar:
nombre_grupo <- unique(grupo$GRUPO)[1]
top5 <- munic_grupo2 %>%
st_drop_geometry() %>%
arrange(desc(max_prod)) %>%
head(5)
kable(top5, caption = paste("Los 5 municipios con mayor producción de",
nombre_grupo, "en el Cesar (2025)"))| mpio_cdpmp | mpio_cnmbr | max_prod |
|---|---|---|
| 20710 | San Alberto | 76326.36 |
| 20770 | San Martín | 36800.00 |
| 20013 | Agustín Codazzi | 30300.00 |
| 20238 | El Copey | 23870.00 |
| 20383 | La Gloria | 20814.00 |
ggplot(top5, aes(x = reorder(mpio_cnmbr, max_prod), y = max_prod)) +
geom_col(fill = "darkgreen") +
coord_flip() +
labs(title = paste("Top 5 municipios en producción de", nombre_grupo),
x = "Municipio", y = "Producción (t)",
caption = "Fuente: UPRA (2025)") +
theme_minimal()
Interpretación: el municipio con mayor producción de oleaginosas en el Cesar es San Alberto, con 76,326 toneladas. Esto representa el 31.1 % de la producción total del departamento para este grupo. Los cinco municipios del top juntos suman el 76.7 %, lo que muestra qué tan concentrada está la producción en pocos municipios.
6. Mapa temático de la producción de cultivos
Ahora hagamos el mapa de producción.
Primero, seleccionamos el modo de visualización para producir un mapa interactivo:
tmap_mode("view")Para las ciudades, uso los centroides de los municipios, un punto por cada municipio:
cesar.cities <- st_centroid(municipios)Luego definimos las capas del mapa y los componentes de la leyenda. Usamos la técnica de cortes naturales de Jenks con 8 clases.
grupo_map =
# CAPA 1: polígonos base con los datos
tm_shape(munic_grupo2) +
# c(longitud, latitud, nivel de zoom) para el Cesar
tm_view(set_view = c(-73.5, 9.5, 8)) +
tm_polygons(
fill = "max_prod",
fill.scale = tm_scale_intervals(
style = "jenks",
n = 8),
fill.legend = tm_legend(title = "Producción (t)")
) +
tm_text(
text = "mpio_cnmbr",
size = 0.55,
fontfamily = "sans",
col = "#222222",
alpha = 0.85
) +
# CAPA 2: vías
tm_shape(cesar.roads.main) +
tm_lines(col = "darkgrey", lwd = 1.5, lty = "solid") +
# CAPA 3: ciudades (centroides de los municipios) en color rojo
tm_shape(cesar.cities) +
tm_symbols(fill = "red", size = 0.4) +
# COMPONENTES DEL DISEÑO
tm_credits("Fuente de los datos: UPRA (2025)", fontface = "bold") +
tm_title(paste("Producción de", tolower(nombre_grupo), "en el Cesar en 2025"),
fontface = "bold.italic",
fontfamily = "monospace") +
tm_scalebar(position = c("left", "bottom"))Ahora veamos el mapa:
grupo_map7. Interpretación del mapa
En el mapa se ve que la producción de oleaginosas principalmente palma de aceite no está repartida de manera uniforme en el Cesar: los municipios con las clases más altas de la leyenda San Alberto, San Martín y Agustín Codazzi concentran la mayor parte de la producción, mientras que los demás quedan en las clases bajas.
Hay 13 municipios sin datos. En este archivo cada municipio aparece con su cultivo principal, así que estos municipios son aquellos donde el cultivo principal pertenece a otro grupo pero no significa que no produzcan nada.
Al superponer las vías y las cabeceras municipales que son los puntos rojos, se puede revisar si los municipios más productivos quedan sobre los corredores viales.
8. Reproducibilidad
sessionInfo()## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22000)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8 LC_CTYPE=Spanish_Colombia.utf8
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] knitr_1.52 ggplot2_4.0.3 tmap_4.4-1 readr_2.2.0 dplyr_1.2.1
## [6] sf_1.1-3
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 xfun_0.62 bslib_0.12.0
## [4] raster_3.6-32 htmlwidgets_1.6.4 lattice_0.22-9
## [7] tzdb_0.5.0 leaflet.providers_3.0.0 vctrs_0.7.3
## [10] tools_4.6.1 crosstalk_1.2.2 generics_0.1.4
## [13] parallel_4.6.1 tibble_3.3.1 proxy_0.4-29
## [16] pkgconfig_2.0.3 KernSmooth_2.23-26 data.table_1.18.6.1
## [19] RColorBrewer_1.1-3 S7_0.2.2 leaflet_2.2.3
## [22] lifecycle_1.0.5 farver_2.1.2 compiler_4.6.1
## [25] terra_1.9-50 codetools_0.2-20 leafsync_0.1.0
## [28] leaflegend_1.3.0 stars_0.7-3 htmltools_0.5.9
## [31] class_7.3-23 sass_0.4.10 yaml_2.3.12
## [34] crayon_1.5.3 pillar_1.11.1 jquerylib_0.1.4
## [37] classInt_0.4-11 cachem_1.1.0 lwgeom_0.2-17
## [40] wk_0.9.5 abind_1.4-8 mime_0.13
## [43] tidyselect_1.2.1 digest_0.6.39 labeling_0.4.3
## [46] fastmap_1.2.0 grid_4.6.1 colorspace_2.1-4
## [49] cli_3.6.6 logger_0.4.3 magrittr_2.0.5
## [52] maptiles_0.13.0 base64enc_0.1-6 XML_3.99-0.25
## [55] cols4all_0.10 leafem_0.2.5 e1071_1.7-17
## [58] withr_3.0.3 scales_1.4.0 bit64_4.8.6
## [61] sp_2.2-4 rmarkdown_2.32 bit_4.6.0
## [64] otel_0.2.0 png_0.1-9 hms_1.1.4
## [67] evaluate_1.0.5 tmaptools_3.3 s2_1.1.13
## [70] rlang_1.3.0 Rcpp_1.1.2 glue_1.8.1
## [73] DBI_1.3.0 vroom_1.7.1 rstudioapi_0.19.0
## [76] jsonlite_2.0.0 R6_2.6.1 spacesXYZ_1.6-0
## [79] units_1.0-19. Referencias
Lizarazo, I. 2026. How to make a map of agricultural production. Disponible en https://rpubs.com/ials2un/MakeEVAmap
UPRA (2025). Evaluaciones Agropecuarias Municipales (EVA).