En este cuaderno vamos a hacer un mapa temático que muestre la participación municipal del grupo de cultivo de Raíces y Tubérculos en el departamento de Antioquia, en el que utilizaremos como fuente principal los archivos csv guardados en el cuaderno EVA, así como un shapefile de municipios obtenido en clase.
rm(list = ls())Primero, como en los anteriores libros, vamos a cargar las librerías necesarias para realizar el trabajo:
library(sf)
library(dplyr)
library(readr)Para comenzar debemos estar seguros de que los archivos que vamos a usar están en el directorio que estamos trabajando:
Un shapefile con los municipios de su departamento, en este caso, los municipios del departamento de Antioquia.
Un archivo csv con las estadísticas EVA 2022 para el grupo de cultivos seleccionado, es decir, Raíces y Tubérculos.
Además, necesitaremos un archivo con las ciudades de Colombia, el cual se puede descargar en el siguiente link: https://simplemaps.com/data/world-cities
Listemos los archivos por orden:
list.files(pattern = c("csv"))[1] "Antiq_Raices_Tuberculos_2022.csv"
[2] "Ciudades_Col.csv"
[3] "EVA_Antioquia.csv"
[4] "worldcities.csv" Ahora veamos los shapefiles guardados en el mismo directorio:
list.files() [1] "Antiq_Raices_Tuberculos_2022.csv"
[2] "Ciudades_Col.csv"
[3] "EVA_Antioquia.csv"
[4] "MA"
[5] "Mapa_Tem.nb.html"
[6] "Mapa_Tem.Rmd"
[7] "Mun_Antiq.cpg"
[8] "Mun_Antiq.dbf"
[9] "Mun_Antiq.prj"
[10] "Mun_Antiq.qmd"
[11] "Mun_Antiq.shp"
[12] "Mun_Antiq.shx"
[13] "simplemaps_worldcities_basicv1.76.zip"
[14] "worldcities.csv"
[15] "worldcities.xlsx" Leeremos los archivos basados en EVA obtenidos del primer cuaderno:
(Raices_Tuberculos = read.csv("Antiq_Raices_Tuberculos_2022.csv"))(Raices_Tuberculos$COD_MUN=as.numeric(Raices_Tuberculos$COD_MUN)) [1] 5480 5264 5686 5664 5674 5893 5895
[8] 5400 5790 5604 5051 5120 5234 5045
[15] 5088 5490 5148 5154 5250 5086 5361
[22] 5736 5837 5670 5756 5495 5697 5440
[29] 5425 5887 5237 5031 5854 5652 5002
[36] 5040 5665 5001 5659 5649 5197 5873
[43] 5206 5890 5607 5034 5101 5172 5266
[50] 5541 5667 5318 5858 5585 5679 5576
[57] 5310 5885 5308 5125 5055 5030 5658
[64] 5647 5579 5036 5147 5079 5282 5591
[71] 5660 5847 5091 5284 5315 5240 5313
[78] 5475 5376 5190 5615 5483 5809 5501
[85] 5021 5360 5142 5145 5792 5093 5364
[92] 5861 5380 5107 5347 5004 5856 5368
[99] 5390class(Raices_Tuberculos$COD_MUN)[1] "numeric"Raices_TuberculosEste archivo lee los municipios del departamento trabajado (Antioquia)
mun.tmp = st_read("./Mun_Antiq.shp")Reading layer `Mun_Antiq' from data source `C:\Users\natal\Desktop\Universidad\Geomatica\RStudio\Mapa_T\Mun_Antiq.shp'
using driver `ESRI Shapefile'
Simple feature collection with 125 features and 9 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -77.12783 ymin: 5.418558 xmax: -73.88128 ymax: 8.873974
Geodetic CRS: WGS 84mun.tmp %>% select(MPIO_CCDGO, MPIO_CNMBR, MPIO_NAREA) -> MunicipiosMunicipiosSimple feature collection with 125 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -77.12783 ymin: 5.418558 xmax: -73.88128 ymax: 8.873974
Geodetic CRS: WGS 84
First 10 features:
MPIO_CCDGO MPIO_CNMBR MPIO_NAREA
1 1 MEDELLÍN 374.83400
2 2 ABEJORRAL 507.14109
3 4 ABRIAQUÍ 296.89405
4 21 ALEJANDRÍA 128.93215
5 30 AMAGÁ 84.13268
6 31 AMALFI 1209.14802
7 34 ANDES 402.46597
8 36 ANGELÓPOLIS 81.87630
9 38 ANGOSTURA 338.50229
10 40 ANORÍ 1413.77572
geometry
1 MULTIPOLYGON (((-75.66974 6...
2 MULTIPOLYGON (((-75.46938 5...
3 MULTIPOLYGON (((-76.08351 6...
4 MULTIPOLYGON (((-75.0332 6....
5 MULTIPOLYGON (((-75.67587 6...
6 MULTIPOLYGON (((-74.92268 6...
7 MULTIPOLYGON (((-75.86822 5...
8 MULTIPOLYGON (((-75.69149 6...
9 MULTIPOLYGON (((-75.27173 6...
10 MULTIPOLYGON (((-74.90935 7...class("Municipios$MPIO_CCDGO")[1] "character"Ahora, lee el archivo csv de las ciudades:
(cities = read_csv("./Ciudades_Col.csv"))Rows: 714 Columns: 11── Column specification ───────────────
Delimiter: ","
chr (7): city, city_ascii, country,...
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.Nótese que los 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):
sf.cities <- st_as_sf(x=cities,
coords = c("lng","lat"))sf.citiesSimple feature collection with 714 features and 9 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -81.7006 ymin: -4.2167 xmax: -67.4667 ymax: 12.5847
CRS: NAObservemos que aún hace falta el CRS.
Vamos a añadirlo
st_crs(sf.cities) <- 4326Aquí, como sólo nos interesa un departamento, tenemos que crear una unión espacial:
sf.cities.joined <- st_join(sf.cities, Municipios, join = st_within)visualizamos el objeto unido:
sf.cities.joinedSimple feature collection with 714 features and 12 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -81.7006 ymin: -4.2167 xmax: -67.4667 ymax: 12.5847
Geodetic CRS: WGS 84Ahora, vamos a filtrar las filas que corresponden a nuestro departamento (en este caso Antioquia):
Antiq.cities = dplyr::filter(sf.cities.joined,admin_name=="Antioquia")Antiq.citiesSimple feature collection with 92 features and 12 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -76.7833 ymin: 5.5833 xmax: -73.9167 ymax: 8.8517
Geodetic CRS: WGS 84Hemos obtenido 92 ciudades que corresponden al departamento de Antioquia.
Ahora, haremos un mapa coropleto de estos datos. Utilizaremos las bibliotecas tmap, ggplot2 y classInt, entonces como hicimos antes, vamos a cargarlas:
# Se ejecutan las siguientes líneas desde la ventana de comandos:
#install.packages("tmap")
#install.packages("ggplot2")
#install.packages("classInt")Ahora las librerías:
library(tmap)library("ggplot2")
library("ggrepel")
library("classInt")Recordemos que el objeto municipios no tiene atributos EVA. En cambio, el objeto de Raíces y Tubérculos, que contiene estadísticas de cultivos, es un objeto no espacial. Lo que sigue es unir el objeto stats a los objetos espaciales para tener los datos relevantes en un único 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, es decir, el codigo municipal. Sin embargo, ten en cuenta que tanto sus nombres como sus tipos de datos son diferentes.
class(Raices_Tuberculos$COD_MUN)[1] "numeric"class(Municipios$MPIO_CCDGO)[1] "character"(Municipios$MPIO_CCDGO=as.character(Municipios$MPIO_CCDGO)) [1] "1" "2" "4" "21" "30"
[6] "31" "34" "36" "38" "40"
[11] "42" "44" "45" "51" "55"
[16] "59" "79" "86" "88" "91"
[21] "93" "101" "107" "113" "120"
[26] "125" "129" "134" "138" "142"
[31] "145" "147" "148" "150" "154"
[36] "172" "190" "197" "206" "209"
[41] "212" "234" "237" "240" "250"
[46] "264" "266" "282" "284" "306"
[51] "308" "310" "313" "315" "318"
[56] "321" "347" "353" "360" "361"
[61] "364" "368" "376" "380" "390"
[66] "400" "411" "425" "440" "467"
[71] "475" "480" "483" "490" "495"
[76] "501" "541" "543" "576" "579"
[81] "585" "591" "604" "607" "615"
[86] "628" "631" "642" "647" "649"
[91] "652" "656" "658" "659" "660"
[96] "664" "665" "667" "670" "674"
[101] "679" "686" "690" "697" "736"
[106] "756" "761" "789" "790" "792"
[111] "809" "819" "837" "842" "847"
[116] "854" "856" "858" "861" "873"
[121] "885" "887" "890" "893" "895"class(Municipios$MPIO_CCDGO)[1] "character"Por lo tanto, tenemos que solucionarlo:
Raices_Tuberculos$COD_MUN = as.character(Raices_Tuberculos$COD_MUN)Ahora si procedemos a unirlos:
Mun_RT = left_join(Municipios, Raices_Tuberculos, by = c("MPIO_CCDGO" = "COD_MUN"))Aquí el resultado:
Mun_RTSimple feature collection with 125 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -77.12783 ymin: 5.418558 xmax: -73.88128 ymax: 8.873974
Geodetic CRS: WGS 84
First 10 features:
MPIO_CCDGO MPIO_CNMBR MPIO_NAREA
1 1 MEDELLÍN 374.83400
2 2 ABEJORRAL 507.14109
3 4 ABRIAQUÍ 296.89405
4 21 ALEJANDRÍA 128.93215
5 30 AMAGÁ 84.13268
6 31 AMALFI 1209.14802
7 34 ANDES 402.46597
8 36 ANGELÓPOLIS 81.87630
9 38 ANGOSTURA 338.50229
10 40 ANORÍ 1413.77572
MUNICIPIO GRUPO max_prod
1 <NA> <NA> NA
2 <NA> <NA> NA
3 <NA> <NA> NA
4 <NA> <NA> NA
5 <NA> <NA> NA
6 <NA> <NA> NA
7 <NA> <NA> NA
8 <NA> <NA> NA
9 <NA> <NA> NA
10 <NA> <NA> NA
geometry
1 MULTIPOLYGON (((-75.66974 6...
2 MULTIPOLYGON (((-75.46938 5...
3 MULTIPOLYGON (((-76.08351 6...
4 MULTIPOLYGON (((-75.0332 6....
5 MULTIPOLYGON (((-75.67587 6...
6 MULTIPOLYGON (((-74.92268 6...
7 MULTIPOLYGON (((-75.86822 5...
8 MULTIPOLYGON (((-75.69149 6...
9 MULTIPOLYGON (((-75.27173 6...
10 MULTIPOLYGON (((-74.90935 7...breaks <- classIntervals(Mun_RT$max_prod, n = 6, style = "fisher")Warning: var has missing values, omitted in finding classesWarning: n greater than number of different finite values\nn reset to number of different finite valuesWarning: n same as number of different finite values\neach different finite value is a separate classError in sVar[1:(length(sVar) - 1)] :
solamente 0's pueden ser mezclados con subscritos negativos