Introducción

En este cuaderno se ilustrara como hacer un mapa tematico que nos muestren la participacion municipal de la producción de cultivos de cereales para el departamento del Tolima. 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.

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)

Lectura de Archivos relacionados con cultivos, municipios y ciudades.

Teniendo en cuenta que el departamento trabajado es Tolima, se deben revisar las rutas de los archivos,los nombres de los archivos y los nombres de las variables. Nos aseguramos que los siguientes archivos esten guardados en el computador: - Un shapefile con los municipios del departamento (Tolima) - Un archivo csv con las estadísticas de EVA 2020 para el grupo de cultivos (Cereales)

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] "Tolima_cereales_2022.csv" "worldcities.csv"

Cuales son los shapefiles guardados en ese directorio.

list.files('BM')
## [1] "Tol_MPIO.shp.CPG" "Tol_MPIO.shp.dbf" "Tol_MPIO.shp.prj" "Tol_MPIO.shp.qix"
## [5] "Tol_MPIO.shp.shp" "Tol_MPIO.shp.shx" "Tol_MPIO.shp.xml"

Ahora, procedamos a leer los archivos basados en EVA obtenidos del primer cuaderno:

(Cereales = read_csv("Tolima_cereales_2022.csv",show_col_types = FALSE))

Este, lee los municipios de el departamento (Tolima).

(mun.tmp =  st_read('BM/Tol_MPIO.shp.shp'))
## Reading layer `Tol_MPIO.shp' from data source 
##   `C:\Users\SEBAS\Downloads\GB2\CUADERNO_EVA\BM\Tol_MPIO.shp.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 46 features and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.313035
## Geodetic CRS:  WGS 84

Ahora, se seleccionan algunos atributos para limpiar el objeto:

mun.tmp %>% select(MPIO_CCNCT, MPIO_CNMBR, MPIO_NAREA) -> municipios
municipios

Ahora, 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:

(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.cities
col.cities

Se filtran para Tolima:

col.cities %>%
  filter(admin_name== "Tolima")->Tol.cities
Tol.cities

Los valores de las coordenadas se almacenan en los atributos “lng” (longitud) y “lat” (latitud).

Para convertir ciudades en un objeto espacial (un objeto de característica simple, en este caso):

sf.cities <-  st_as_sf(x = Tol.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) <- 4326

4. Datos de subconjunto relevantes para nuestro departamento.

Como estamos interesados en un solo departamento (Tolima), necesitamos crear una unión espacial:

sf.cities.joined <- st_join(sf.cities, municipios, join = st_within)

Para mostrar el objeto unido:

sf.cities.joined

Tenga 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 Tolima):

Tol.ct = dplyr::filter(sf.cities.joined, admin_name=='Tolima')
Tol.ct

5. Hacer un Mapa para el grupo De cultivo Hortalizas como el que mas importancia le hemos dado.

Para hacer el mapa es necesario instalar los paquetes necesarios:

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 Cereales, 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.

class(Cereales$COD_MUN)
## [1] "numeric"
class(municipios$MPIO_CCNCT)
## [1] "character"

Por lo tanto, necesitamos arreglar esto:

Cereales$COD_MUN = as.character(Cereales$COD_MUN)
head(municipios)

Ahora si se tiene todo listo para hacer la union:

munic_Cereales = left_join(municipios, Cereales, by = c("MPIO_CCNCT" = "COD_MUN"))

Ahora necesitamos reemplazar 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_Cereales2 <- munic_Cereales %>% replace(is.na(.),0)

Y el resultado es:

munic_Cereales2

Tenga en cuenta que munic_Cereales2 incluye ahora el atributo max_prod que representa la cantidad total de toneladas producidas por cultivos de Cereales 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_Cereales2$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_Cereales2 <-  munic_Cereales2 %>%
  mutate(faktor_class = factor(cut(max_prod, breaks$brks, include.lowest = T), labels = lab_vec))

munic_Cereales2$Produccion = munic_Cereales2$faktor_clas
munic_Cereales2$mid <- sf::st_centroid(munic_Cereales2$geometry)

Obtenemos longitud:

LONG = st_coordinates(munic_Cereales2$mid)[,1]

Obtenemos latitud:

LAT = st_coordinates(munic_Cereales2$mid)[,2]
ggplot(data = munic_Cereales2) +
   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: 9 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_Cereales2) + 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 Cereales 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 viewing
Ht_map
## Credits not supported in view mode.
## Symbol shapes other than circles or icons are not supported in view mode.

6. Haz un nuevo mapa para el segundo grupo más grande de cultivos.

Por ultimo podemos expresar la produccion maxima de nuestros cultivos de cereales, por medio de un mapa mas completo, interactivo y atractivo para su uso.

# see example at https://r-tmap.github.io/tmap-book/layout.html
facet = "max_prod"
  Cereales_map =  
  tm_shape(munic_Cereales2) + tm_polygons(facet) + tm_text(text = "MPIO_CNMBR", size = 0.7, fontfamily = "sans") +
  tm_shape(Tol.ct) + tm_symbols(shape = 2, col = "red", size = 0.20) +
  tm_credits("Data source: UPRA (2020)", fontface = "bold") +
  tm_layout(main.title = "Produccion maxima de Cereales 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 viewing
Cereales_map
## Credits not supported in view mode.
## Symbol shapes other than circles or icons are not supported in view mode.

7. Bibliografia

Lizarazo, I., 2022. Disponible en https://rpubs.com/ials2un/thematic_maps_v2

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.4   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              highr_0.10               pkgconfig_2.0.3         
## [19] KernSmooth_2.23-21       RColorBrewer_1.1-3       leaflet_2.2.0           
## [22] lifecycle_1.0.3          farver_2.1.1             compiler_4.3.1          
## [25] munsell_0.5.0            terra_1.7-55             codetools_0.2-19        
## [28] leafsync_0.1.0           stars_0.6-4              htmltools_0.5.6         
## [31] class_7.3-22             sass_0.4.7               yaml_2.3.7              
## [34] pillar_1.9.0             crayon_1.5.2             jquerylib_0.1.4         
## [37] ellipsis_0.3.2           cachem_1.0.8             lwgeom_0.2-13           
## [40] wk_0.8.0                 abind_1.4-5              mime_0.12               
## [43] tidyselect_1.2.0         digest_0.6.33            fastmap_1.1.1           
## [46] grid_4.3.1               colorspace_2.1-0         cli_3.6.1               
## [49] magrittr_2.0.3           base64enc_0.1-3          dichromat_2.0-0.1       
## [52] XML_3.99-0.14            utf8_1.2.3               leafem_0.2.3            
## [55] e1071_1.7-13             withr_2.5.0              scales_1.2.1            
## [58] sp_2.1-0                 bit64_4.0.5              rmarkdown_2.24          
## [61] bit_4.0.5                png_0.1-8                hms_1.1.3               
## [64] evaluate_0.21            knitr_1.43               tmaptools_3.1-1         
## [67] viridisLite_0.4.2        s2_1.1.4                 rlang_1.1.1             
## [70] Rcpp_1.0.11              glue_1.6.2               DBI_1.1.3               
## [73] rstudioapi_0.15.0        vroom_1.6.3              jsonlite_1.8.7          
## [76] R6_2.5.1                 units_0.8-3