Mapa de producción de Hortalizas en el departamento de Boyacá
Gabriel Andres Gutierrez Gomez
Introducción
En este cuaderno se ilustra como hacer mapas tematicos que nos muestren la participacion municipal de la producción de cultivos de hortalizas para el departamento de Boyacá. Se usa como fuente principal de datos, el archivo csv, guardado en el cuaderno EVA, así como el shapefile de los municipios obtenidos en clase.
Configuración
Se instalaran y cargaran las bibliotecas de R necesarias.
#run the following lines from the command line:
#install.packages('dplyr')
#install.packages('readxl')
#install.packages('sf)library(sf)
library(dplyr)
library(readr)Lectura de Archivos relacionados con cultivos, municipios y ciudades.
Teniendo en cuenta que el mucipio trabajado es bayacá, ses deben revisar las rutas de los archivos,los nombres de los archivos y los nombres de las variables. Se asegura que los siguientes archivos esten guardados en el computador: - Un shapefile con los municipios del departamento (Boyacá) - Un archivo csv con las estadísticas de EVA 2020 para el grupo de cultivos (Hortalizas)
Además de esto, se necesitaran un archivo con ciudades de Colombia, este se puede descargar desde https://simplemaps.com/data/co-cities. Este archivo descargado queda de forma co.csv y se lleva al el directorio de trabajo. Se enumeran los archivos disponibles, y se verifican cuales son los de tipo csv guardados en el directorio de trabajo.
list.files( pattern=c('csv'))## [1] "Boy_Hortalizas_2022.csv" "co.csv"
## [3] "Eva_Boyaca.csv" "worldcities.csv"Cuales son los shapefiles guardados en ese directorio.
#Se cambia la siguiente linea en el directorio de datos#
list.files('BM')## [1] "Boy_Mun.cpg" "Boy_Mun.dbf" "Boy_Mun.prj" "Boy_Mun.qmd" "Boy_Mun.shp"
## [6] "Boy_Mun.shx"Ahora, procedamos a leer los archivos basados en EVA obtenidos del primer cuaderno:
(Hortalizas = read_csv("Boy_Hortalizas_2022.csv",show_col_types = FALSE))Este, lee los municipios de el departamento (Boyacá).
# este archivo shapefile debe tener un código EPSG 4326
# este archivo shapefile se obtuvo en clase usando QGIS
(mun.tmp = st_read('BM/Boy_Mun.shp'))## Reading layer `Boy_Mun' from data source
## `C:\Users\gagug\OneDrive\Escritorio\GEOMATICA-20230824T155303Z-001\GEOMATICA\GB2\CuaerdoEVA\BM\Boy_Mun.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 246 features and 8 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS: WGS 84Ahora, se seleccionan algunos atributos para limpiar el objeto:
# comprobar la salida del objeto en el último fragmento y
# cambiar los nombres de los atributos según sus propios datos
mun.tmp %>% select(MPIO_CCDGO, MPIO_CNMBR, AREA) -> municipiosmunicipiosAhora, se lee el archivo cvs de las ciudades, para este es necesario ir filtrando desde las ciudades del mundo hasta llegar a las de nuestro departamento:
# Este archivo se descargó de simplemaps como se describe en la parte de arriba
(cities = read_csv("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.Se filtran para Colombia:
cities %>%
filter(iso3== "COL")->col.citiescol.citiesSe filtran para Boyacá:
col.cities %>%
filter(admin_name== "Boyacá")->Boy.citiesBoy.citiesLos valores de las coordenadas se almacenan en los atributos “lng” (latitud) y “lat” (longitud).
Para convertir ciudades en un objeto espacial (un objeto de característica simple, en este caso):
#Convertir data frame a un objeto de caracteristica simple
sf.cities <- st_as_sf(x = Boy.cities,
coords = c("lng", "lat"))sf.cities¡Tenga en cuenta el CRS que falta!
Agreguemos el CRS usando el código EPSG:
st_crs(sf.cities) <- 43264. Datos de subconjunto relevantes para nuestro departamento.
Como estamos interesados en un solo departamento (Boyacá), necesitamos crear una unión espacial:
# buscar puntos (ciudades) dentro de polígonos (nuestros municipios)
sf.cities.joined <- st_join(sf.cities, municipios, join = st_within)Para mostrar el objeto unido:
sf.cities.joinedTenga en cuenta que obtuvimos un marco de datos sf con cada fila de ciudades adjuntadas con las columnas de nuestros municipios. Las ciudades ubicadas en un departamento diferente tienen muchas NA.
Ahora filtramos las filas que corresponden a nuestro departamento (en este caso Boyacá):
Boy.ct = dplyr::filter(sf.cities.joined, admin_name=='Boyacá')Boy.ct5. Hacer un Mapa para el grupo De cultivo Hortalizas como el que mas importancia le hemos dado.
Para hacer rl mapa es necesario installar los paquetes necesarios:
# Se ejecutan las siguientes líneas desde la ventana de comandos:
#install.packages("tmap")
#install.packages("ggplot2")
#install.packages("classInt")Ahora las librerias:
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 Hortalizas, 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.
### Revisamos de acuerdo con nuestros datos propios
class(Hortalizas$COD_MUN)## [1] "numeric"class(municipios$MPIO_CCDGO)## [1] "character"Por lo tanto, necesitamos arreglar esto:
Hortalizas$COD_MUN = as.character(Hortalizas$COD_MUN)Pero antes de unirlos debemos darnos cuenta que el atributo de la union sea igual y al revisar en nuestra tabla municipios el codigo de la ciudad se presenta solo con las 3 ultimas cifras, esto se arregla compensando las cifras faltantes por medio de la siguiente función:
municipios$str = "15"municipios$codigo=paste(municipios$str, municipios$MPIO_CCDGO, sep = "")head(municipios)Ahora si se tiene todo listo para hacer la union:
munic_Hortalizas = left_join(municipios, Hortalizas, by = c("codigo" = "COD_MUN"))Ahora necesitamos remplazar todos aquellos datos que aparacen como no existentes NA por valores de 0 pra que no intervengan en la creación de nuestro mapa, esto se hace por medio de:
munic_Hortalizas2 <- munic_Hortalizas %>% replace(is.na(.),0)Y el resultado es:
munic_Hortalizas2Tenga en cuenta que munic_Hortalizas2 incluye ahora el atributo max_prod que representa la cantidad total de toneladas producidas por cultivos de Hortalizas en 2022.
El siguiente código personaliza los colores de los municipios. Más información [aquí] (https://stackoverflow.com/questions/57177312/trying-to-plot-in-tmap-shapefile-with-attribute).
breaks <- classIntervals(munic_Hortalizas2$max_prod, n = 6, style = 'fisher')
#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 para el mapa sintactico:
munic_Hortalizas2 <- munic_Hortalizas2 %>%
mutate(faktor_class = factor(cut(max_prod, breaks$brks, include.lowest = T), labels = lab_vec))
# Cambiar el nombre del atributo
munic_Hortalizas2$Produccion = munic_Hortalizas2$faktor_clas# Este código crea un nuevo campo ("mind") con coordenadas centroides
munic_Hortalizas2$mid <- sf::st_centroid(munic_Hortalizas2$geometry)# Obtener los valores de longitud
LONG = st_coordinates(munic_Hortalizas2$mid)[,1]# Obtener los valores de latitud
LAT = st_coordinates(munic_Hortalizas2$mid)[,2]ggplot(data = munic_Hortalizas2) +
geom_sf(aes(fill = Produccion)) +
geom_label_repel(aes(x = LONG, y = LAT, label = MPIO_CNMBR),
label.padding = unit(0.05,"lines"),
label.r = unit(0.025, "lines"),
label.size = 0.05)## Warning: ggrepel: 238 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
El mapa no muetra la información de muchos municipios mas alla de denotar su producción por medio de una barra de colores, es por eso que se ha creado el siguiente mapa:
# see example at https://r-tmap.github.io/tmap-book/layout.html
facet = "Produccion"
Ht_map =
tm_shape(munic_Hortalizas2) + tm_polygons(facet) + tm_text(text = "MUNICIPIO", size = 0.7, fontfamily = "sans") +
tm_shape(municipios) + tm_symbols(shape = 2, col = "red", size = 0.20) +
tm_credits("Data source: UPRA (2020)", fontface = "bold") +
tm_layout(main.title = "Produccion de Hortalizas 2022",
main.title.fontface = "bold.italic",
legend.title.fontfamily = "monospace") +
tm_scale_bar(position = c("left", "bottom"))
tmap_mode("view")## tmap mode set to interactive viewingHt_map## Credits not supported in view mode.## Symbol shapes other than circles or icons are not supported in view mode.