Mapa de producción de Frutales en el departamento de Huila

1.Introducción

Este es el segundo cuaderno R Notebook en el curso Geomatica Básica 2023 dictada por el ing. Ivan Lizarazo.Aqui ilustraremos cómo hacer mapas temáticos que muestren la participación municipal de los grupos de cultivos más importantes para un el departamento de Huila. Usaremos como fuente principal los archivos csv guardados en el cuaderno EVA, así como un shapefile de municipios obtenidos en clase.

2.Configuración

En este paso, instalaremos y cargaremos las bibliotecas R necesarias.

#ejecute las siguientes líneas desde la línea de comando:
#install.packages('dplyr')
#install.packages('readxl')
#install.packages('sf')

library(sf)
library(dplyr)
library(readr)

3.Leer archivos relacionados con municipios, cultivos y ciudades.

Antes de comenzar, asegúrese de haber guardado los siguientes archivos en su computadora: 1.Un archivo shapefile con los municipios de su departamento (como se produjo en clase) 2.Un archivo csv con estadísticas EVA 2020 para el grupo de cultivos seleccionado (como se produjo en el primer cuaderno )

Además, es necesario tener un expediente con las principales ciudades de Colombia. Puede descargarlo desde https://simplemaps.com/data/co-cities. Luego, debe descomprimirlo y moverlo a su directorio de trabajo.

Comencemos enumerando los archivos disponibles.

¿Cuáles son los archivos csv guardados en el directorio de trabajo?

list.files("./", pattern=c('csv'))

## [1] "co.csv"                  "huila_frutales_2022.csv"
## [3] "worldcities.csv"

Cuales son los shapefiles guardados en ese directorio.

#Se cambia la siguiente linea en el directorio de datos#
list.files()

##  [1] "01.qgz"                                                  
##  [2] "Clase27"                                                 
##  [3] "co.csv"                                                  
##  [4] "Diccionario_Datos_Niveles_Variables_MGN_CNPV2018Int.xlsx"
##  [5] "Geo_Huila.gpkg"                                          
##  [6] "huila_frutales_2022.csv"                                 
##  [7] "HUILA_MUN.gpkg"                                          
##  [8] "MapasTematicos.nb.html"                                  
##  [9] "MapasTematicos.Rmd"                                      
## [10] "MapaTematico.qgz"                                        
## [11] "MGN_ANM_MPIOS.cpg"                                       
## [12] "MGN_ANM_MPIOS.dbf"                                       
## [13] "MGN_ANM_MPIOS.prj"                                       
## [14] "MGN_ANM_MPIOS.qmd"                                       
## [15] "MGN_ANM_MPIOS.sbn"                                       
## [16] "MGN_ANM_MPIOS.sbx"                                       
## [17] "MGN_ANM_MPIOS.shp"                                       
## [18] "MGN_ANM_MPIOS.shp.xml"                                   
## [19] "MGN_ANM_MPIOS.shx"                                       
## [20] "NUEVOHUILA.qgz"                                          
## [21] "NUEVOS_HUILA.cpg"                                        
## [22] "NUEVOS_HUILA.dbf"                                        
## [23] "NUEVOS_HUILA.prj"                                        
## [24] "NUEVOS_HUILA.qmd"                                        
## [25] "NUEVOS_HUILA.shp"                                        
## [26] "NUEVOS_HUILA.shx"                                        
## [27] "SHP_MGN2018_INTGRD_MPIO.zip"                             
## [28] "worldcities.csv"

Ahora procedemos a leer los archivos EVA creados en el anterior cuaderno de R

(frutales = read_csv("huila_frutales_2022.csv",show_col_types = FALSE))

Éste código lee los municipios de el departamento de Huila creado en QGIS en clases anteriores:

# este archivo shapefile debe tener un código EPSG 4326
(mun.tmp =  st_read('NUEVOS_HUILA.shp'))

## Reading layer `NUEVOS_HUILA' from data source 
##   `D:\Sexta Matricula\Geomatica\Mapas tematicos\NUEVOS_HUILA.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 37 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.62466 ymin: 1.552125 xmax: -74.41303 ymax: 3.843208
## Geodetic CRS:  MAGNA-SIRGAS

Ahora, 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(DPTO_CCDGO,MPIO_CCDGO, MPIO_CNMBR, AREA) -> 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 los municipios de nuestro departamento:

#Este archivo se descargó de *simplemaps* como se decribio anteriormente
(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 filtraron los datos para Colombia:

cities %>%
  filter(iso3=="COL")-> col.cites

col.cites

Filtramos por nuetro departamento (Huila).

col.cites %>%
  filter(admin_name == "Huila") -> Huila.cites

Huila.cites

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

#Convertir data frame a un objeto de caracteristica simple
sf.cities <-  st_as_sf(x = Huila.cites, 
                        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 el departamento.

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

sf.cities <- st_transform(sf.cities, st_crs(municipios))

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 NAs.

Ahora filtramos las filas que corresponden a nuestro departamento (en este caso Huila):

huila.cities = dplyr::filter(sf.cities.joined, admin_name=='Huila')

huila.cities

5. Haz un mapa para el grupo de cultivos más importante.

Ahora, haz un mapa coroplético de estos datos. Usaremos las bibliotecas tmap, ggplot2 y classInt.

# 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 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.

### Revisamos de acuerdo con nuestros datos propios
class(frutales$COD_MUN)

## [1] "numeric"

class(municipios$MPIO_CCDGO)

## [1] "character"

frutales$COD_MUN = as.character(frutales$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 = "41"

municipios$codigo=paste(municipios$str, municipios$MPIO_CCDGO, sep = "")

head(municipios)

Ahora si se tiene todo listo para hacer la union:

munic_frutales = left_join(municipios, frutales, 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_frutales2 <- munic_frutales %>% replace(is.na(.),0)

munic_frutales2

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

# Este código crea un nuevo campo ("mind") con coordenadas centroides
munic_frutales2$mid <- sf:: st_centroid(munic_frutales2$geometry)

# Obtener los valores de longitud
LONG = st_coordinates(munic_frutales2$mid) [,1]

# Obtener los valores de latitud
LAT = st_coordinates(munic_frutales2$mid) [,2]

ggplot(data = munic_frutales2) +
  geom_sf(aes(fill = max_prod)) +
  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: 2 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:

 #mira un ejempo en 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(huila.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.

Como vemos en el mapa creado ,Produccion de frutales en 2020, el municipio de Garzón es el que presenta mayor prducción anaual, seguido de Pitalito.

frutales_map <- iconv(frutales_map, from = "ISO-8859-1", to = "UTF-8")

Bibliografia

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

Puede encontrar este cuaderno en el siguiente link < >

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              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.25           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-4