Mapas de producción maxima de cereales y frutales en el ejecafetero
Nicolas Cely
2023-05-06
Introducción
Este es el cuaderno R Markdown numero 2 del curso de geomatica basica 2023. este ilustra como hacer mapas teamticos que muestren la participación municipal de los dos grupos mas importantes para los departamentos de Caldas, Quindio y Risaralda. Se usara como fuente principal de información los archivos CSV guardados en el cuaderno EVA, así como un shapefile de municipios obtenidos en clase.
Configuración
Aqui se instalan y cargan todas las librerias necesarias para realizar este trabajo
library(dplyr)
library(readr)
library(sf)## Warning: package 'sf' was built under R version 4.2.3Lectura de archivos
Cargar archivos relacionados con muncipios, cultivos y ciudades de nuestros de departamentos. En este se trabajaran los departamentos de Caldas, Risaralda y Quindio.
Antes de empezar se deben tener los siguientes archivos en el computador
- Un shapefile con los municipios de los departamentos de interes.
- Un archivo CSV con las estadisticas EVA 2020
ademas se necesita un archivo con las ciudades de Colombia. Puede descargarse en https://simplemaps.com/data/co-cities.
Primero enumeraremos los archivos disponibles
Archivos CSV guardados en el directorio de trabajo:
list.files("./Datos", pattern=c('csv'))## [1] "co.csv" "eje_cereales.csv" "eje_frutales.csv" "EVA.csv"Archivos shapefile guardados en el directorio de trabajo
list.files("./Datos", pattern=c('shp'))## [1] "muni_eje.shp"Ahora, se leeran los archivos obtenidos a partir del EVA en cuaderno 1.
(Cereales = read_csv("C:/Users/LENOVO/Desktop/Geomatica/R/Datos/eje_cereales.csv",show_col_types = FALSE))## New names:
## • `` -> `...1`## # A tibble: 53 × 5
## ...1 grupo cod_mun mun prodmax
## <dbl> <chr> <dbl> <chr> <dbl>
## 1 1 CEREALES 17380 LA DORADA 27840
## 2 2 CEREALES 66001 PEREIRA 5000
## 3 3 CEREALES 17616 RISARALDA 3400
## 4 4 CEREALES 17042 ANSERMA 1778
## 5 5 CEREALES 63470 MONTENEGRO 1566
## 6 6 CEREALES 17174 CHINCHINA 1315
## 7 7 CEREALES 66687 SANTUARIO 1150
## 8 8 CEREALES 17877 VITERBO 1080
## 9 9 CEREALES 63401 LA TEBAIDA 1071
## 10 10 CEREALES 17001 MANIZALES 720
## # … with 43 more rowsLeer los municipios de los departamentos
(mun.eje = st_read('C:/Users/LENOVO/Desktop/Geomatica/R/Datos/muni_eje.shp'))## Reading layer `muni_eje' from data source
## `C:\Users\LENOVO\Desktop\Geomatica\R\Datos\muni_eje.shp' using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 14 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -76.3595 ymin: 4.135801 xmax: -74.6489 ymax: 5.7626
## Geodetic CRS: WGS 84## Simple feature collection with 51 features and 14 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -76.3595 ymin: 4.135801 xmax: -74.6489 ymax: 5.7626
## Geodetic CRS: WGS 84
## First 10 features:
## fid ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 TYPE_2 ENGTYPE_2
## 1 1 53 COL Colombia 8 Caldas 348 Aguadas Municipio Municipality
## 2 2 53 COL Colombia 8 Caldas 349 Anserma Municipio Municipality
## 3 3 53 COL Colombia 8 Caldas 350 Aranzazú Municipio Municipality
## 4 4 53 COL Colombia 8 Caldas 351 Belalcázar Municipio Municipality
## 5 5 53 COL Colombia 8 Caldas 352 Chinchiná Municipio Municipality
## 6 6 53 COL Colombia 8 Caldas 353 Filadelfia Municipio Municipality
## 7 7 53 COL Colombia 8 Caldas 354 La Dorada Municipio Municipality
## 8 8 53 COL Colombia 8 Caldas 355 La Merced Municipio Municipality
## 9 9 53 COL Colombia 8 Caldas 356 Manizales Municipio Municipality
## 10 10 53 COL Colombia 8 Caldas 357 Manzanares Municipio Municipality
## NL_NAME_2 VARNAME_2 Area codigo geometry
## 1 <NA> <NA> 501.64063 17013 POLYGON ((-75.3399 5.458, -...
## 2 <NA> <NA> 196.73857 17616 POLYGON ((-75.6908 5.152099...
## 3 <NA> Aranzazé 140.63766 17050 POLYGON ((-75.3997 5.2202, ...
## 4 <NA> <NA> 114.27085 17088 POLYGON ((-75.7841 5.032401...
## 5 <NA> <NA> 125.39270 17174 POLYGON ((-75.6204 4.979101...
## 6 <NA> <NA> 202.22182 17272 POLYGON ((-75.5888 5.207699...
## 7 <NA> <NA> 583.68634 17380 POLYGON ((-74.8504 5.290101...
## 8 <NA> <NA> 90.16428 17388 POLYGON ((-75.6437 5.362199...
## 9 <NA> <NA> 443.87241 17001 POLYGON ((-75.636 5.0011, -...
## 10 <NA> <NA> 171.87523 17433 POLYGON ((-75.1662 5.1453, ...Seleccionaremos algunos atributos para limpiar el objeto
mun.eje %>% select(codigo, NAME_2, Area) -> municipioshead (municipios)## Simple feature collection with 6 features and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -75.9068 ymin: 4.9149 xmax: -75.3399 ymax: 5.730701
## Geodetic CRS: WGS 84
## codigo NAME_2 Area geometry
## 1 17013 Aguadas 501.6406 POLYGON ((-75.3399 5.458, -...
## 2 17616 Anserma 196.7386 POLYGON ((-75.6908 5.152099...
## 3 17050 Aranzazú 140.6377 POLYGON ((-75.3997 5.2202, ...
## 4 17088 Belalcázar 114.2709 POLYGON ((-75.7841 5.032401...
## 5 17174 Chinchiná 125.3927 POLYGON ((-75.6204 4.979101...
## 6 17272 Filadelfia 202.2218 POLYGON ((-75.5888 5.207699...cargar el archivo CSV de las ciudades
(cities = read_csv("./Datos/co.csv"))## Rows: 1102 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): city, country, iso2, admin_name, capital
## dbl (4): lat, lng, population, population_proper
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.## # A tibble: 1,102 × 9
## city lat lng country iso2 admin_name capital popul…¹ popul…²
## <chr> <dbl> <dbl> <chr> <chr> <chr> <chr> <dbl> <dbl>
## 1 Bogotá 4.61 -74.1 Colombia CO Bogotá primary 9464000 7963000
## 2 Medellín 6.24 -75.6 Colombia CO Antioquia admin 2529403 2529403
## 3 Cali 3.44 -76.5 Colombia CO Valle del Ca… admin 2471474 2471474
## 4 Barranquilla 11.0 -74.8 Colombia CO Atlántico admin 1274250 1274250
## 5 Cartagena 10.4 -75.5 Colombia CO Bolívar admin 1036412 1036412
## 6 Soacha 4.58 -74.2 Colombia CO Cundinamarca minor 995268 995268
## 7 Palermo 2.89 -75.4 Colombia CO Huila minor 800000 800000
## 8 Cúcuta 7.91 -72.5 Colombia CO Norte de San… admin 750000 750000
## 9 Soledad 10.9 -74.8 Colombia CO Atlántico minor 698852 342556
## 10 Pereira 4.81 -75.7 Colombia CO Risaralda admin 590554 590554
## # … with 1,092 more rows, and abbreviated variable names ¹population,
## # ²population_properConvertir ciudades en un objeto espacial
sf.cities <- st_as_sf(x = cities,
coords = c("lng", "lat"))head (sf.cities)## Simple feature collection with 6 features and 7 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -76.5197 ymin: 3.44 xmax: -74.0705 ymax: 10.9639
## CRS: NA
## # A tibble: 6 × 8
## city country iso2 admin…¹ capital popul…² popul…³ geometry
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <POINT>
## 1 Bogotá Colomb… CO Bogotá primary 9464000 7963000 (-74.0705 4.6126)
## 2 Medel… Colomb… CO Antioq… admin 2529403 2529403 (-75.5748 6.2447)
## 3 Cali Colomb… CO Valle … admin 2471474 2471474 (-76.5197 3.44)
## 4 Barra… Colomb… CO Atlánt… admin 1274250 1274250 (-74.7964 10.9639)
## 5 Carta… Colomb… CO Bolívar admin 1036412 1036412 (-75.5253 10.4236)
## 6 Soacha Colomb… CO Cundin… minor 995268 995268 (-74.2144 4.5781)
## # … with abbreviated variable names ¹admin_name, ²population,
## # ³population_properAñadir CRS
st_crs(sf.cities) <- 4326Agrupar datos relevantes
Crear union espacial entre el sf.cities y municipios
sf.cities.joined <- st_join(sf.cities, municipios, join = st_within)head(sf.cities.joined)## Simple feature collection with 6 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -76.5197 ymin: 3.44 xmax: -74.0705 ymax: 10.9639
## Geodetic CRS: WGS 84
## # A tibble: 6 × 11
## city country iso2 admin…¹ capital popul…² popul…³ geometry
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <POINT [°]>
## 1 Bogotá Colomb… CO Bogotá primary 9464000 7963000 (-74.0705 4.6126)
## 2 Medel… Colomb… CO Antioq… admin 2529403 2529403 (-75.5748 6.2447)
## 3 Cali Colomb… CO Valle … admin 2471474 2471474 (-76.5197 3.44)
## 4 Barra… Colomb… CO Atlánt… admin 1274250 1274250 (-74.7964 10.9639)
## 5 Carta… Colomb… CO Bolívar admin 1036412 1036412 (-75.5253 10.4236)
## 6 Soacha Colomb… CO Cundin… minor 995268 995268 (-74.2144 4.5781)
## # … with 3 more variables: codigo <dbl>, NAME_2 <chr>, Area <dbl>, and
## # abbreviated variable names ¹admin_name, ²population, ³population_properFiltrar filas corespondientes a los departamentos de Caldas, Quindio y risaralda.
municipios.eje = dplyr::filter(sf.cities.joined, admin_name=='Caldas'|
admin_name=='Risaralda'|
admin_name=='Quindio')head(municipios.eje)## Simple feature collection with 6 features and 10 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -75.7028 ymin: 4.8143 xmax: -74.6647 ymax: 5.4538
## Geodetic CRS: WGS 84
## # A tibble: 6 × 11
## city country iso2 admin…¹ capital popul…² popul…³ geometry
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <POINT [°]>
## 1 Perei… Colomb… CO Risara… admin 590554 590554 (-75.6946 4.8143)
## 2 Maniz… Colomb… CO Caldas admin 434403 434403 (-75.4847 5.0661)
## 3 Dosqu… Colomb… CO Risara… minor 179301 179301 (-75.6761 4.8361)
## 4 La Do… Colomb… CO Caldas minor 78540 78540 (-74.6647 5.4538)
## 5 Santa… Colomb… CO Risara… minor 73028 73028 (-75.6211 4.8672)
## 6 Riosu… Colomb… CO Caldas minor 62296 62296 (-75.7028 5.4208)
## # … with 3 more variables: codigo <dbl>, NAME_2 <chr>, Area <dbl>, and
## # abbreviated variable names ¹admin_name, ²population, ³population_properHacer un mapa con los grupos de cultivo mas importantes
instalar las librerias ““tmap”“,”“ggplot2”““classint”
library(htmltools)## Warning: package 'htmltools' was built under R version 4.2.3library(tmap)## Warning: package 'tmap' was built under R version 4.2.3library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.2.3library(ggrepel)## Warning: package 'ggrepel' was built under R version 4.2.3library(classInt)## Warning: package 'classInt' was built under R version 4.2.3Es importante hacer una union entre “municipios” y “cereales”, ya que este primero no tiene atributos de EVA. Lo anterior permite obtener datos importantes en un unico objeto.
El atributo en comun de “municipios” y “cereales” es el codigo municipal, sin embargo, se debe revisar que este sea del mismo tipo para ambos objetos.
Clase de codigo municipal en “municipios”
class(municipios$codigo)## [1] "numeric"Clase codigo municipal en “cereales”
class(Cereales$cod_mun)## [1] "numeric"“municipios y”cereales” presentan el mismo tipo de objeto, por ende se pueden unir
union:
munycereales = left_join(municipios, Cereales, by = c("codigo" = "cod_mun"))head (munycereales)## Simple feature collection with 6 features and 7 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -75.9068 ymin: 4.9149 xmax: -75.3399 ymax: 5.730701
## Geodetic CRS: WGS 84
## codigo NAME_2 Area ...1 grupo mun prodmax
## 1 17013 Aguadas 501.6406 43 CEREALES AGUADAS 93
## 2 17616 Anserma 196.7386 3 CEREALES RISARALDA 3400
## 3 17050 Aranzazú 140.6377 29 CEREALES ARANZAZU 256
## 4 17088 Belalcázar 114.2709 11 CEREALES BELALCAZAR 675
## 5 17174 Chinchiná 125.3927 6 CEREALES CHINCHINA 1315
## 6 17272 Filadelfia 202.2218 23 CEREALES FILADELFIA 385
## geometry
## 1 POLYGON ((-75.3399 5.458, -...
## 2 POLYGON ((-75.6908 5.152099...
## 3 POLYGON ((-75.3997 5.2202, ...
## 4 POLYGON ((-75.7841 5.032401...
## 5 POLYGON ((-75.6204 4.979101...
## 6 POLYGON ((-75.5888 5.207699...Ahora munycereales , tiene el atributo de prodmax, el cual representa la cantidad total de toneladas producidas por cultivos de cereales en 2020.
Con la union de objetos obtenida, otorgaremos colores a los municipios con el siguinete codigo.
breaks <- classIntervals(munycereales$prodmax, n = 6, style = 'fisher')
lab_vec <- vector(length = length(breaks$brks)-1)
rounded_breaks <- round(breaks$brks,2)
lab_vec[1] <- paste0('[', rounded_breaks[1],' - ', rounded_breaks[2],']')
for(i in 2:(length(breaks$brks) - 1)){
lab_vec[i] <- paste0('(',rounded_breaks[i], ' - ', rounded_breaks[i+1], ']')
}Ahora, haremos el mapa estatico
munycereales <- munycereales %>%
mutate(faktor_class = factor(cut(prodmax, breaks$brks, include.lowest = T), labels = lab_vec))
# cambiar el nombre del atributo
munycereales$prodmax = munycereales$faktor_classObtención de coordenadas del centroide (mid)
munycereales$mid <- sf::st_centroid(munycereales$geometry)obtención de los valores de longitud
LONG = st_coordinates(munycereales$mid)[,1]Obtención de los valores de latitud
LAT = st_coordinates(munycereales$mid)[,2]Visualización mapa
ggplot(data = munycereales) +
geom_sf(aes(fill = prodmax)) +
geom_label_repel(aes(x = LONG, y = LAT, label = munycereales$NAME_2),
label.padding = unit(0.05,"lines"),
label.r = unit(0.025, "lines"),
label.size = 0.05)## Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
Segunda mapa para grupo de cultivo importante en la región del eje cafetero
En este apartado se utilizaran los datos del EVA para el grupo cultivable de los frutales.
Cargar archivo:
(Frutales = read_csv("C:/Users/LENOVO/Desktop/Geomatica/R/Datos/eje_frutales.csv",show_col_types = FALSE))## New names:
## • `` -> `...1`## # A tibble: 53 × 5
## ...1 grupo cod_mun mun prodmax
## <dbl> <chr> <dbl> <chr> <dbl>
## 1 1 FRUTALES 63470 MONTENEGRO 49410
## 2 2 FRUTALES 63401 LA TEBAIDA 29750
## 3 3 FRUTALES 17042 ANSERMA 29120
## 4 4 FRUTALES 17513 PACORA 26460
## 5 5 FRUTALES 66001 PEREIRA 24560
## 6 6 FRUTALES 63594 QUIMBAYA 20343
## 7 7 FRUTALES 17174 CHINCHINA 20316
## 8 8 FRUTALES 17380 LA DORADA 19734
## 9 9 FRUTALES 63548 PIJAO 19691
## 10 10 FRUTALES 17013 AGUADAS 16845
## # … with 43 more rowsse debe conocer la clase del atributo con el que desea hacer el join, en este caso es codigo municipal
class(Frutales$cod_mun)## [1] "numeric"class(municipios$codigo)## [1] "numeric"En esta ocasión los dos atributos de interes son numericos, por ende, se puede realizar la union de estos.
Union:
muniyfrutales = left_join(municipios, Frutales, by = c("codigo" = "cod_mun"))head (muniyfrutales)## Simple feature collection with 6 features and 7 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -75.9068 ymin: 4.9149 xmax: -75.3399 ymax: 5.730701
## Geodetic CRS: WGS 84
## codigo NAME_2 Area ...1 grupo mun prodmax
## 1 17013 Aguadas 501.6406 10 FRUTALES AGUADAS 16845
## 2 17616 Anserma 196.7386 16 FRUTALES RISARALDA 11100
## 3 17050 Aranzazú 140.6377 32 FRUTALES ARANZAZU 3000
## 4 17088 Belalcázar 114.2709 33 FRUTALES BELALCAZAR 2940
## 5 17174 Chinchiná 125.3927 7 FRUTALES CHINCHINA 20316
## 6 17272 Filadelfia 202.2218 43 FRUTALES FILADELFIA 1200
## geometry
## 1 POLYGON ((-75.3399 5.458, -...
## 2 POLYGON ((-75.6908 5.152099...
## 3 POLYGON ((-75.3997 5.2202, ...
## 4 POLYGON ((-75.7841 5.032401...
## 5 POLYGON ((-75.6204 4.979101...
## 6 POLYGON ((-75.5888 5.207699...construcción del mapa
facet = "prodmax"
frut_map =
tm_shape(muniyfrutales) + tm_polygons(facet) + tm_text(text = "NAME_2", size = 0.7, fontfamily = "sans") +
tm_shape(municipios.eje) + tm_symbols(shape = 2, col = "blue", size = 0.20) +
tm_credits("Data source: UPRA (2020)", fontface = "bold") +
tm_layout(main.title = "Produccion de frutales en 2020",
main.title.fontface = "bold.italic",
legend.title.fontfamily = "monospace") +
tm_scale_bar(position = c("left", "bottom"))Visualización
tmap_mode("view")## tmap mode set to interactive viewingfrut_map## Credits not supported in view mode.## Symbol shapes other than circles or icons are not supported in view mode.Reproducibilidad
Este trabajo se realizó tomando como guía el cuaderno: Lizarazo, I. 2022, My second R Notebook: How to make thematic maps.
Cite este trabajo así: Cely, N., 2023, Mapas de producción maxima de cereales y frutales en el ejecafetero. Disponible en https://rpubs.com/ncely/1038734
sessionInfo()## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
##
## Matrix products: default
##
## 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
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] classInt_0.4-9 ggrepel_0.9.3 ggplot2_3.4.2 tmap_3.3-3
## [5] htmltools_0.5.5 sf_1.0-12 readr_2.1.4 dplyr_1.1.0
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.2 bit64_4.0.5 vroom_1.6.1
## [4] jsonlite_1.8.4 viridisLite_0.4.0 bslib_0.4.0
## [7] highr_0.9 sp_1.6-0 yaml_2.3.5
## [10] pillar_1.8.1 lattice_0.20-45 glue_1.6.2
## [13] digest_0.6.29 RColorBrewer_1.1-3 leaflet.providers_1.9.0
## [16] colorspace_2.0-3 XML_3.99-0.14 pkgconfig_2.0.3
## [19] raster_3.6-20 stars_0.6-1 s2_1.1.3
## [22] scales_1.2.0 terra_1.7-29 tzdb_0.3.0
## [25] tibble_3.1.8 proxy_0.4-27 generics_0.1.3
## [28] farver_2.1.1 ellipsis_0.3.2 cachem_1.0.6
## [31] withr_2.5.0 leafsync_0.1.0 cli_3.3.0
## [34] mime_0.12 magrittr_2.0.3 crayon_1.5.2
## [37] evaluate_0.15 fansi_1.0.3 lwgeom_0.2-11
## [40] class_7.3-20 tools_4.2.2 hms_1.1.2
## [43] lifecycle_1.0.3 stringr_1.5.0 munsell_0.5.0
## [46] compiler_4.2.2 jquerylib_0.1.4 e1071_1.7-13
## [49] rlang_1.1.0 units_0.8-2 grid_4.2.2
## [52] tmaptools_3.1-1 dichromat_2.0-0.1 rstudioapi_0.14
## [55] htmlwidgets_1.6.2 crosstalk_1.2.0 leafem_0.2.0
## [58] base64enc_0.1-3 rmarkdown_2.14 wk_0.7.2
## [61] gtable_0.3.0 codetools_0.2-18 abind_1.4-5
## [64] DBI_1.1.3 markdown_1.1 R6_2.5.1
## [67] knitr_1.39 fastmap_1.1.0 bit_4.0.5
## [70] utf8_1.2.2 KernSmooth_2.23-20 stringi_1.7.8
## [73] parallel_4.2.2 Rcpp_1.0.10 vctrs_0.5.2
## [76] png_0.1-8 leaflet_2.1.2 tidyselect_1.2.0
## [79] xfun_0.31