Este cuaderno va a Ilustrar acerca de cómo hacer mapas temáticos que muestren la participación municipal de dos grupos de cultivos del departamento de Santander (frutales y leguminosas). Usaremos como fuente principal los archivos csv guardados en el cuaderno EVA, así como un shapefile de municipios obtenidos en clase.
Primero, es necesario instalar las librerias que usaremos.
#run the following lines from the command line:
#install.packages('dplyr')
#install.packages('readxl')
#install.packages('sf)
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(readr)
Antes de comenzar es necesario asegurarse de que se tiene guardados estos archivos en la computadora. Un archivo shapefiles con los municipios de su departamento. Un archivo csv con estadísticas de EVA 2020 para el grupo de cultivos seleccionado. Además, necesitarás un archivo con ciudades de Colombia. Puede descargarlo desde https://simplemaps.com/data/co-cities. Luego, mueva el archivo co.csv descargado a su directorio de trabajo.
¿Cuáles son los archivos csv guardados en el directorio de trabajo?
list.files("./datos", pattern=c('csv'))
## [1] "co.csv" "EVA_STDR.csv"
## [3] "stder_frutales_2022.csv" "stder_leguminosas_2022.csv"
## [5] "worldcities.csv"
¿Cuáles son los archivos de shapefiles guardados en el directorio de trabajo?
list.files("./datos")
## [1] "BaseEVA_Agricola20192022.xlsx"
## [2] "co.csv"
## [3] "Cuaderno.otroGrupo.Rmd"
## [4] "EVA_STDR.csv"
## [5] "license.txt"
## [6] "Mun-colombia.dbf"
## [7] "Mun-colombia.qmd"
## [8] "santander-mun.cpg"
## [9] "santander-mun.dbf"
## [10] "santander-mun.prj"
## [11] "santander-mun.shp"
## [12] "santander-mun.shx"
## [13] "SANTANDER_EVA.qgz"
## [14] "simplemaps_worldcities_basicv1.76 (1).zip"
## [15] "stder_frutales_2022.csv"
## [16] "stder_leguminosas_2022.csv"
## [17] "worldcities.csv"
## [18] "worldcities.xlsx"
Ahora, vamos a leer los archivos EVA obtenidos en el primer cuaderno:
# adjust the filepath acording to your dato
(frutales = read_csv("./datos/stder_frutales_2022.csv",show_col_types = FALSE))
Ahora lee los municipios del departamento
# this shapefile must have a 4326 EPSG code
# this shapefile was obtained in class using QGIS
(mun.tmp = st_read('./datos/santander-mun.shp'))
## Reading layer `santander-mun' from data source
## `D:\CuadernoAguacate\datos\santander-mun.shp' using driver `ESRI Shapefile'
## Simple feature collection with 87 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.52895 ymin: 5.719626 xmax: -72.5161 ymax: 8.119248
## Geodetic CRS: MAGNA-SIRGAS
Seleccionemos algunos atributos para limpiar el objeto:
# check the objet output in the last chunk and
# change attribute names according to your own data
mun.tmp %>% select(MPIO_CDPMP, MPIO_CNMBR, AREA) -> municipios
municipios
Ahora, vamos a leer el archivo csv de las ciudades del mundo:
# This file was downloaded from simplemaps as described above
(worldcities = read_csv("./datos/worldcities.csv"))
## Rows: 44691 Columns: 11
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (7): city, city_ascii, country, iso2, iso3, admin_name, capital
## dbl (4): lat, lng, population, id
##
## ℹ 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.
Para convertir ciudades en un objeto espacial (un objeto de característica simple, en este caso
# Convert data frame to simple feature object
worldcities <- st_as_sf(x = worldcities,
coords = c("lng", "lat"))
worldcities
Vamos a filtrar las ciudades del mundo para colombia.
colombia.cities = dplyr::filter(worldcities, iso3=='COL')
colombia.cities
Agreguemos el CRS usando el código EPSG:
st_crs(colombia.cities) <- 4326
colombia.cities <- st_transform(colombia.cities, st_crs(municipios))
Como estamos interesados en un solo departamento, necesitamos crear una unión espacial:
# find points (cities) within polygons (our municipalities)
sf.cities.joined <- st_join(colombia.cities, municipios, join = st_within)
Para mostrarlos juntos
sf.cities.joined
Ahora filtramos las filas que corresponden a nuestro departamento: Santander.
stder.cities = dplyr::filter(colombia.cities, admin_name=='Santander')
stder.cities
Obtuvimos 28 ciudades ubicadas en Santander.
Ahora, realizaremos un mapa coroplético de estos datos. Usaremos las bibliotecas tmap, ggplot2 y classInt.
# run the following lines from the command window:
#install.packages("tmap")
#install.packages("ggplot2")
#install.packages("classInt")
library(tmap)
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, were retired in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(ggplot2)
library(ggrepel)
library(classInt)
Recordemos que el objeto municipios no tiene atributos EVA. En cambio, el objeto frutales, que contiene estadísticas de cultivos, es un objeto no espacial. Nuestra siguiente tarea es unir el objeto de estadísticas a los objetos espaciales para tener los datos relevantes en un solo objeto.
Para poder realizar la unión, necesitamos una clave compartida, es decir, un atributo común. En este caso lo tenemos en ambos objetos, ese es el código municipal. Sin embargo, tenga en cuenta que tanto sus nombres como sus tipos de datos son diferentes.
### check according to your own data
class(frutales$COD_MUN)
## [1] "numeric"
class(municipios$MPIO_CDPMP)
## [1] "character"
Vamo a arreglarlo
frutales$COD_MUN = as.character(frutales$COD_MUN)
Ahora vamos a unirlos
munic_frutales = left_join(municipios, frutales, by = c("MPIO_CDPMP" = "COD_MUN"))
Resultado
munic_frutales
El siguiente código personaliza los colores de los municipios. Más información
breaks <- classIntervals(munic_frutales$max_prod, n = 6, style = 'fisher')
## Warning in classIntervals(munic_frutales$max_prod, n = 6, style = "fisher"):
## var has missing values, omitted in finding classes
#label breaks
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, el mapa estatico
munic_frutales <- munic_frutales %>%
mutate(faktor_class = factor(cut(max_prod, breaks$brks, include.lowest = T), labels = lab_vec))
# change attribute name
munic_frutales$Produccion = munic_frutales$faktor_class
# This code creates a new field ("mid") with centroid coordinates
munic_frutales$mid <- sf::st_centroid(munic_frutales$geometry)
# Get the longitude values
LONG = st_coordinates(munic_frutales$mid)[,1]
# Get the latitude values
LAT = st_coordinates(munic_frutales$mid)[,2]
ggplot(data = munic_frutales) +
geom_sf(aes(fill = Produccion)) +
geom_label_repel(aes(x = LONG, y = LAT, label = MUNICIPIO),
label.padding = unit(0.05,"lines"),
label.r = unit(0.025, "lines"),
label.size = 0.05)
## Warning: Removed 1 rows containing missing values (`geom_label_repel()`).
## Warning: ggrepel: 49 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
# see example at https://r-tmap.github.io/tmap-book/layout.html
facet = "max_prod"
frutales_map =
tm_shape(munic_frutales) + tm_polygons(facet) + tm_text(text = "MPIO_CNMBR", size = 0.7, fontfamily = "sans") +
tm_shape(stder.cities) + tm_symbols(shape = 2, col = "red", 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"))
tmap_mode("view")
## tmap mode set to interactive viewing
frutales_map
## Credits not supported in view mode.
## Symbol shapes other than circles or icons are not supported in view mode.
frutales_map <- iconv(frutales_map, from = "ISO-8859-1", to = "UTF-8")
Utilizaremos los datos del EVA del grupo de “Leguminosas”. Para este paso necesitamos realizar pasos similares a los anteriores:
(leguminosas = read_csv("./datos/stder_leguminosas_2022.csv",show_col_types = FALSE))
leguminosas$COD_MUN = as.character(leguminosas$COD_MUN)
#municipios$MPIO_CDPMP = as.character(municipios$MPIO_CDPMP)
munic_leguminosas = left_join(municipios, leguminosas, by = c("MPIO_CDPMP" = "COD_MUN"))
munic_leguminosas
# see example at https://r-tmap.github.io/tmap-book/layout.html
facet = "max_prod"
leguminosas_map =
tm_shape(munic_leguminosas) + tm_polygons(facet) + tm_text(text = "MPIO_CNMBR", size = 0.7, fontfamily = "sans") +
tm_shape(stder.cities) + tm_symbols(shape = 2, col = "red", size = 0.20) +
tm_credits("Data source: UPRA (2020)", fontface = "bold") +
tm_layout(main.title = "Produccion de leguminosas en 2020",
main.title.fontface = "bold.italic",
legend.title.fontfamily = "monospace") +
tm_scale_bar(position = c("left", "bottom"))
tmap_mode("view")
## tmap mode set to interactive viewing
leguminosas_map
## Credits not supported in view mode.
## Symbol shapes other than circles or icons are not supported in view mode.
sessionInfo()
## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
##
## 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
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] classInt_0.4-10 ggrepel_0.9.3 ggplot2_3.4.3 tmap_3.3-4
## [5] readr_2.1.4 dplyr_1.1.3 sf_1.0-14
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.4 xfun_0.40 bslib_0.5.1
## [4] raster_3.6-23 htmlwidgets_1.6.2 lattice_0.21-8
## [7] tzdb_0.4.0 leaflet.providers_1.13.0 vctrs_0.6.3
## [10] tools_4.3.1 crosstalk_1.2.0 generics_0.1.3
## [13] parallel_4.3.1 tibble_3.2.1 proxy_0.4-27
## [16] fansi_1.0.4 pkgconfig_2.0.3 KernSmooth_2.23-21
## [19] RColorBrewer_1.1-3 leaflet_2.2.0 lifecycle_1.0.3
## [22] farver_2.1.1 compiler_4.3.1 munsell_0.5.0
## [25] terra_1.7-46 codetools_0.2-19 leafsync_0.1.0
## [28] stars_0.6-4 htmltools_0.5.6 class_7.3-22
## [31] sass_0.4.7 yaml_2.3.7 pillar_1.9.0
## [34] crayon_1.5.2 jquerylib_0.1.4 ellipsis_0.3.2
## [37] cachem_1.0.8 lwgeom_0.2-13 wk_0.8.0
## [40] abind_1.4-5 mime_0.12 tidyselect_1.2.0
## [43] digest_0.6.33 fastmap_1.1.1 grid_4.3.1
## [46] colorspace_2.1-0 cli_3.6.1 magrittr_2.0.3
## [49] base64enc_0.1-3 dichromat_2.0-0.1 XML_3.99-0.14
## [52] utf8_1.2.3 leafem_0.2.3 e1071_1.7-13
## [55] withr_2.5.0 scales_1.2.1 sp_2.1-0
## [58] bit64_4.0.5 rmarkdown_2.24 bit_4.0.5
## [61] png_0.1-8 hms_1.1.3 evaluate_0.21
## [64] knitr_1.44 tmaptools_3.1-1 viridisLite_0.4.2
## [67] s2_1.1.4 rlang_1.1.1 Rcpp_1.0.11
## [70] glue_1.6.2 DBI_1.1.3 rstudioapi_0.15.0
## [73] vroom_1.6.3 jsonlite_1.8.7 R6_2.5.1
## [76] units_0.8-3