Mapa de producción de Hortalizas en el departamento de Boyacá

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"                  "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"
[4] "Boy_Mun.qmd" "Boy_Mun.shp" "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 84
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 84
First 10 features:
   fid DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CDPMP      AREA  LATITUD  LONGITUD                       geometry
1    1         15        001      TUNJA      15001 119688918 5.518473 -73.37802 POLYGON ((-73.36346 5.55805...
2    2         15        022    ALMEIDA      15022  57672119 4.954825 -73.38813 POLYGON ((-73.36793 5.01349...
3    3         15        047  AQUITANIA      15047 942146563 5.437416 -72.87149 POLYGON ((-72.76242 5.63856...
4    4         15        051   ARCABUCO      15051 137898588 5.749565 -73.43888 POLYGON ((-73.50487 5.84347...
5    5         15        087      BELEN      15087 163088220 6.005059 -72.89370 POLYGON ((-72.91692 6.08612...
6    6         15        090     BERBEO      15090  58013016 5.232058 -73.10117 POLYGON ((-73.0677 5.27048,...
7    7         15        092  BETEITIVA      15092 101899548 5.920859 -72.84858 POLYGON ((-72.81796 5.97422...
8    8         15        097    BOAVITA      15097 145305291 6.337516 -72.62021 POLYGON ((-72.64907 6.43640...
9    9         15        104     BOYACA      15104  48022868 5.439856 -73.38137 POLYGON ((-73.34806 5.47411...
10  10         15        106    BRICENO      15106  64599703 5.675510 -73.90933 POLYGON ((-73.89118 5.73749...

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(MPIO_CCDGO, MPIO_CNMBR, AREA) -> municipios

municipios

Simple feature collection with 246 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
Geodetic CRS:  WGS 84
First 10 features:
   MPIO_CCDGO MPIO_CNMBR      AREA                       geometry
1         001      TUNJA 119688918 POLYGON ((-73.36346 5.55805...
2         022    ALMEIDA  57672119 POLYGON ((-73.36793 5.01349...
3         047  AQUITANIA 942146563 POLYGON ((-72.76242 5.63856...
4         051   ARCABUCO 137898588 POLYGON ((-73.50487 5.84347...
5         087      BELEN 163088220 POLYGON ((-72.91692 6.08612...
6         090     BERBEO  58013016 POLYGON ((-73.0677 5.27048,...
7         092  BETEITIVA 101899548 POLYGON ((-72.81796 5.97422...
8         097    BOAVITA 145305291 POLYGON ((-72.64907 6.43640...
9         104     BOYACA  48022868 POLYGON ((-73.34806 5.47411...
10        106    BRICENO  64599703 POLYGON ((-73.89118 5.73749...

Ahora, se lee el archivo cvs de las ciudades:

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

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

col.cities

col.cities %>%
  filter(admin_name== "Boyacá")->Boy.cities

Boy.cities

También 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):

#Convertir data frame a un objeto de caracteristica simple
sf.cities <-  st_as_sf(x = Boy.cities, 
                        coords = c("lng", "lat"))

sf.cities

Simple feature collection with 29 features and 9 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -74.4167 ymin: 5.0056 xmax: -72.4167 ymax: 6.1667
CRS:           NA

¡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 (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.joined

Simple feature collection with 58 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -74.4167 ymin: 5.0056 xmax: -72.4167 ymax: 6.1667
Geodetic CRS:  WGS 84

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

Boy.ct = dplyr::filter(sf.cities.joined, admin_name=='Boyacá')

Boy.ct

Simple feature collection with 58 features and 12 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -74.4167 ymin: 5.0056 xmax: -72.4167 ymax: 6.1667
Geodetic CRS:  WGS 84

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

# Se ejecutan las siguientes líneas desde la ventana de comandos:
#install.packages("tmap")
#install.packages("ggplot2")
#install.packages("classInt")

library(tmap)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
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)

Ahora estamos listos para unirlos:

municipios$str = "15"

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

head(municipios)

Simple feature collection with 6 features and 5 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -73.51402 ymin: 4.904038 xmax: -72.65243 ymax: 6.098348
Geodetic CRS:  WGS 84
  MPIO_CCDGO MPIO_CNMBR      AREA                       geometry str codigo
1        001      TUNJA 119688918 POLYGON ((-73.36346 5.55805...  15  15001
2        022    ALMEIDA  57672119 POLYGON ((-73.36793 5.01349...  15  15022
3        047  AQUITANIA 942146563 POLYGON ((-72.76242 5.63856...  15  15047
4        051   ARCABUCO 137898588 POLYGON ((-73.50487 5.84347...  15  15051
5        087      BELEN 163088220 POLYGON ((-72.91692 6.08612...  15  15087
6        090     BERBEO  58013016 POLYGON ((-73.0677 5.27048,...  15  15090

munic_Hortalizas = left_join(municipios, Hortalizas, by = c("codigo" = "COD_MUN"))

Y el resultado es:

munic_Hortalizas2

Simple feature collection with 246 features and 10 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
Geodetic CRS:  WGS 84
First 10 features:
   MPIO_CCDGO MPIO_CNMBR      AREA str codigo MUNICIPIO      GRUPO max_prod                       geometry
1         001      TUNJA 119688918  15  15001     Tunja Hortalizas     5180 POLYGON ((-73.36346 5.55805...
2         022    ALMEIDA  57672119  15  15022   Almeida Hortalizas      180 POLYGON ((-73.36793 5.01349...
3         047  AQUITANIA 942146563  15  15047 Aquitania Hortalizas    65250 POLYGON ((-72.76242 5.63856...
4         051   ARCABUCO 137898588  15  15051  Arcabuco Hortalizas       40 POLYGON ((-73.50487 5.84347...
5         087      BELEN 163088220  15  15087     Belén Hortalizas      765 POLYGON ((-72.91692 6.08612...
6         090     BERBEO  58013016  15  15090    Berbeo Hortalizas      560 POLYGON ((-73.0677 5.27048,...
7         092  BETEITIVA 101899548  15  15092 Betéitiva Hortalizas      256 POLYGON ((-72.81796 5.97422...
8         097    BOAVITA 145305291  15  15097   Boavita Hortalizas      750 POLYGON ((-72.64907 6.43640...
9         104     BOYACA  48022868  15  15104    Boyacá Hortalizas     4500 POLYGON ((-73.34806 5.47411...
10        106    BRICENO  64599703  15  15106         0          0        0 POLYGON ((-73.89118 5.73749...
      faktor_class      Produccion
1  (1860 - 6062.5] (1860 - 6062.5]
2       [0 - 1860]      [0 - 1860]
3  (58305 - 65250] (58305 - 65250]
4       [0 - 1860]      [0 - 1860]
5       [0 - 1860]      [0 - 1860]
6       [0 - 1860]      [0 - 1860]
7       [0 - 1860]      [0 - 1860]
8       [0 - 1860]      [0 - 1860]
9  (1860 - 6062.5] (1860 - 6062.5]
10            <NA>            <NA>

Tenga en cuenta que munic_Hortalizas 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).

munic_Hortalizas2 <- munic_Hortalizas %>% replace(is.na(.),0)

Warning: invalid factor level, NA generatedWarning: invalid factor level, NA generated

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)

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

# see example at https://r-tmap.github.io/tmap-book/layout.html
facet = "max_prod"
cereal_map =  
  tm_shape(munic_Hortalizas2) + tm_polygons(facet) + tm_text(text = "MPIO_CNMBR", size = 0.7, fontfamily = "sans") +
  tm_shape(cities) + tm_symbols(shape = 2, col = "red", size = 0.20) +
  tm_credits("Data source: UPRA (2020)", fontface = "bold") +
  tm_layout(main.title = "Produccion de oleaginosas 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