1. Introducción

La agricultura en el departamento de Boyacá, es una de las principales actividades económicas junto con la producción pecuaria y la minería. Boyacá, rico en recursos naturales, con una infraestructura energética, vial y de servicios adecuada y una posición geográfica privilegiada ha centrado su actividad económica en la agricultura tradicional, con un alto porcentaje de la producción dedicada a abastecer el mercado interno, se caracteriza por ser uno de los departamentos con mayor producción papera en el país, aunque también se cultiva maíz, cebada, cebolla, caña, café, frutales, hortalizas y muchos más. Este es un cuaderno de R Markdown que ilustra las estadísticas agricolas para el departamento de Boyacá en Colombia.

2. Funcionalidades SIG

La exploración de estadísticas no espaciales es esencial para comprender lo que está sucediendo en los territorios. Varias librerías R, en particular dplyr y tidyverse, son muy útiles para explorar y resumir estadísticas.

Por otro lado, las operaciones geoespaciales pueden mejorar nuestra comprensión de varios problemas que afectan a las regiones geográficas. Por ejemplo, desea averiguar cuál es la ubicación de aquellos municipios cuyos rendimientos de cosecha son sobresalientes (o, alternativamente, más bajos que el promedio). Para realizar dicha exploración, necesitamos unir datos no espaciales con datos espaciales. Además, también podríamos explorar uniones espaciales. Estas operaciones se basan en la intersección entre dos objetos espaciales, a menudo puntos y los polígonos.

Comencemos por eliminar el contenido de la memoria:

rm(list=ls())

Ahora, instalemos las librerías que necesitamos (en el siguiente fragmento, cualquier paquete se instala solo si no se ha instalado previamente.):

list.of.packages <- c("here", "tidyverse", "rgeos", "maptools", "raster", "sf",  "viridis", "rnaturalearth", "GSODR", "ggrepel", "cowplot")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

Ahora, carguemos las librerías:

library(here)
## here() starts at C:/Users/Brian/Desktop/Daniela/Agricultura
library(tidyverse)
## -- Attaching packages ------------- tidyverse 1.3.0 --
## v ggplot2 3.3.0     v purrr   0.3.3
## v tibble  2.1.3     v dplyr   0.8.5
## v tidyr   1.0.2     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0
## -- Conflicts ---------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(rgeos)
## Loading required package: sp
## rgeos version: 0.5-2, (SVN revision 621)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 
##  Linking to sp version: 1.4-1 
##  Polygon checking: TRUE
library(maptools)
## Checking rgeos availability: TRUE
library(raster)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:tidyr':
## 
##     extract
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(viridis)
## Loading required package: viridisLite
library(rnaturalearth)
library(GSODR)
library(ggrepel)
library(cowplot)
## 
## ********************************************************
## Note: As of version 1.0.0, cowplot does not change the
##   default ggplot2 theme anymore. To recover the previous
##   behavior, execute:
##   theme_set(theme_cowplot())
## ********************************************************

3. Explorando estadísticas agrícolas en Boyacá

Con los datos estadísticos ya descargados en formato csv, como Estadisticas Municipales Agropecuarias en mi computadora, he usado Excel para eliminar las filas correspondientes a los departamentos diferentes de Boyacá, guardando el archivo correspondiente como un archivo local con el nombre EVA_Boyacá.csv. Hecho esto, leamos el archivo csv “estadisticas municipales agropecuarias” para Boyacá, así:

datos <- read_csv2("C:/Users/Brian/Desktop/Daniela/Agricultura/EVA_Boyacá.csv")
## Using ',' as decimal and '.' as grouping mark. Use read_delim() for more control.
## Parsed with column specification:
## cols(
##   COD_DEP = col_double(),
##   DEPARTAMENTO = col_character(),
##   COD_MUN = col_double(),
##   MUNICIPIO = col_character(),
##   GRUPO_DE_CULTIVO = col_character(),
##   SUBGRUPO_DE_CULTIVO = col_character(),
##   CULTIVO = col_character(),
##   YEAR = col_double(),
##   `Area_Sembrada_(ha)` = col_double(),
##   `Area_Cosechada_(ha)` = col_double(),
##   `Produccion_(t)` = col_double(),
##   `Rendimiento_(t/ha)` = col_double(),
##   ESTADO_FISICO_PRODUCCION = col_character(),
##   CICLO_DE_CULTIVO = col_character()
## )

Veamos cuáles son los atributos de los datos:

head(datos)
## # A tibble: 6 x 14
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO_DE_CULTIVO SUBGRUPO_DE_CUL~
##     <dbl> <chr>          <dbl> <chr>     <chr>            <chr>           
## 1      15 BOYACA         15114 BUSBANZA  HORTALIZAS       ACELGA          
## 2      15 BOYACA         15516 PAIPA     HORTALIZAS       ACELGA          
## 3      15 BOYACA         15491 NOBSA     HORTALIZAS       ACELGA          
## 4      15 BOYACA         15500 OICATA    HORTALIZAS       ACELGA          
## 5      15 BOYACA         15516 PAIPA     HORTALIZAS       ACELGA          
## 6      15 BOYACA         15500 OICATA    HORTALIZAS       ACELGA          
## # ... with 8 more variables: CULTIVO <chr>, YEAR <dbl>,
## #   `Area_Sembrada_(ha)` <dbl>, `Area_Cosechada_(ha)` <dbl>,
## #   `Produccion_(t)` <dbl>, `Rendimiento_(t/ha)` <dbl>,
## #   ESTADO_FISICO_PRODUCCION <chr>, CICLO_DE_CULTIVO <chr>
tail(datos)
## # A tibble: 6 x 14
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO_DE_CULTIVO SUBGRUPO_DE_CUL~
##     <dbl> <chr>          <dbl> <chr>     <chr>            <chr>           
## 1      15 BOYACA         15514 PAEZ      FRUTALES         ZAPOTE          
## 2      15 BOYACA         15753 SOATA     FRUTALES         ZAPOTE          
## 3      15 BOYACA         15533 PAYA      FRUTALES         ZAPOTE          
## 4      15 BOYACA         15514 PAEZ      FRUTALES         ZAPOTE          
## 5      15 BOYACA         15753 SOATA     FRUTALES         ZAPOTE          
## 6      15 BOYACA         15533 PAYA      FRUTALES         ZAPOTE          
## # ... with 8 more variables: CULTIVO <chr>, YEAR <dbl>,
## #   `Area_Sembrada_(ha)` <dbl>, `Area_Cosechada_(ha)` <dbl>,
## #   `Produccion_(t)` <dbl>, `Rendimiento_(t/ha)` <dbl>,
## #   ESTADO_FISICO_PRODUCCION <chr>, CICLO_DE_CULTIVO <chr>

Cada municipio tiene estadísticas sobre el área sembrada, área cosechada, y rendimiento para diferentes cultivos en diferentes años. Los atributos SUBGRUPO y CULTIVO parecen referirse a la misma cosa (es decir, un cultivo). Los cultivos se clasifican además en un GRUPO dado. En esta tabla, tenemos unidades y si verificamos el archivo csv original, encontramos que las unidades de área son hectáreas y que las unidades de rendimiento son Ton / ha, entre otros.

Utilizaremos la librería dplyr para explorar el contenido del objeto de datos. Primero, obtengamos un resumen del rendimiento (es decir, el rendimiento promedio durante varios años) por grupo y municipio:

datos %>%
  group_by(MUNICIPIO, GRUPO_DE_CULTIVO) %>%
  summarise(rend_prom = mean(`Rendimiento_(t/ha)`, na.rm = TRUE)) -> rend_resumen
### Let's visualize only the first six records
head(rend_resumen)
## # A tibble: 6 x 3
## # Groups:   MUNICIPIO [1]
##   MUNICIPIO GRUPO_DE_CULTIVO      rend_prom
##   <chr>     <chr>                     <dbl>
## 1 ALMEIDA   CEREALES                   36.9
## 2 ALMEIDA   FRUTALES                   45.2
## 3 ALMEIDA   HORTALIZAS                 80.4
## 4 ALMEIDA   LEGUMINOSAS                61.4
## 5 ALMEIDA   OTROS PERMANENTES          70.7
## 6 ALMEIDA   TUBERCULOS Y PLATANOS      70.9

También podemos calcular el rendimiento promedio por GRUPO en el departamento de Boyacá:

datos %>%
  group_by(GRUPO_DE_CULTIVO) %>%
  summarise(rend_dep = mean(`Rendimiento_(t/ha)`, na.rm = TRUE)) -> rend_Boyaca

rend_Boyaca
## # A tibble: 12 x 2
##    GRUPO_DE_CULTIVO                                 rend_dep
##    <chr>                                               <dbl>
##  1 CEREALES                                            33.1 
##  2 FIBRAS                                              24.5 
##  3 FLORES Y FOLLAJES                                   32.0 
##  4 FRUTALES                                            52.9 
##  5 HONGOS                                               1   
##  6 HORTALIZAS                                         164.  
##  7 LEGUMINOSAS                                         43.2 
##  8 OLEAGINOSAS                                         26   
##  9 OTROS PERMANENTES                                   49.6 
## 10 OTROS TRANSITORIOS                                  32.4 
## 11 PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDICINALES     8.66
## 12 TUBERCULOS Y PLATANOS                              157.

Se observa que los rendimientos más altos corresponden a HORTALIZAS, FRUTALES y TUBERCULOS Y PLATANOS. Ahora, veamos cuáles son los municipios con mayor rendimiento para cada grupo de cultivos en 2018:

datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRUPO_DE_CULTIVO, MUNICIPIO) %>%
  summarize(max_rend = max(`Rendimiento_(t/ha)`, na.rm = TRUE)) %>%
    slice(which.max(max_rend)) -> rend_max_18

rend_max_18
## # A tibble: 11 x 3
## # Groups:   GRUPO_DE_CULTIVO [11]
##    GRUPO_DE_CULTIVO                                 MUNICIPIO     max_rend
##    <chr>                                            <chr>            <dbl>
##  1 CEREALES                                         CIENEGA            311
##  2 FIBRAS                                           RONDON               6
##  3 FLORES Y FOLLAJES                                FIRAVITOBA         105
##  4 FRUTALES                                         TUNUNGUA          1653
##  5 HORTALIZAS                                       GUATEQUE          6105
##  6 LEGUMINOSAS                                      GARAGOA           1067
##  7 OLEAGINOSAS                                      PUERTO BOYACA       35
##  8 OTROS PERMANENTES                                BOAVITA            256
##  9 OTROS TRANSITORIOS                               SUSACON            205
## 10 PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDICINALES SOATA               15
## 11 TUBERCULOS Y PLATANOS                            PAZ DE RIO        3011

Ahora, veamos cuáles son los municipios con mayor área cosechada para cada grupo de cultivos en 2018:

datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRUPO_DE_CULTIVO, MUNICIPIO) %>%
  summarize(max_area_cosecha = max(`Area_Cosechada_(ha)`, na.rm = TRUE)) %>%
    slice(which.max(max_area_cosecha)) -> area_cosecha_max

area_cosecha_max
## # A tibble: 11 x 3
## # Groups:   GRUPO_DE_CULTIVO [11]
##    GRUPO_DE_CULTIVO                                MUNICIPIO    max_area_cosecha
##    <chr>                                           <chr>                   <dbl>
##  1 CEREALES                                        LA VICTORIA               440
##  2 FIBRAS                                          RAQUIRA                    20
##  3 FLORES Y FOLLAJES                               TOCA                       32
##  4 FRUTALES                                        MONIQUIRA                 700
##  5 HORTALIZAS                                      AQUITANIA                1210
##  6 LEGUMINOSAS                                     SORACA                   1600
##  7 OLEAGINOSAS                                     PUERTO BOYA~               15
##  8 OTROS PERMANENTES                               CHITARAQUE               7528
##  9 OTROS TRANSITORIOS                              SOATA                      41
## 10 PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDICINAL~ SOATA                      25
## 11 TUBERCULOS Y PLATANOS                           SABOYA                   3000

Ahora bien, el área más grande cosechada en 2018 correspondió a OTROS PERMANENTES en CHITARAQUE. Buscando información se sabe que Chitaraque es un municipio eminentemente agropecuario, dedicado de manera preferencial al cultivo de la caña panelera y a la fabricación y comercialización de panela, siendo el primer cultivador departamental de éste producto. Se destaca a continuación el cultivo de café, plátano, yuca y algunos frutales.

Dicho esto, seleccionemos la producción (t) de caña en Chitaraque para cada año (se modificó en el excel caña por cana para evitar errores en RStudio):

datos %>% 
  filter(MUNICIPIO=="CHITARAQUE" & SUBGRUPO_DE_CULTIVO=="CANA") %>% 
  group_by(YEAR, CULTIVO) ->  chitaraque_cana

chitaraque_cana
## # A tibble: 13 x 14
## # Groups:   YEAR, CULTIVO [13]
##    COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO_DE_CULTIVO SUBGRUPO_DE_CUL~
##      <dbl> <chr>          <dbl> <chr>     <chr>            <chr>           
##  1      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
##  2      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
##  3      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
##  4      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
##  5      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
##  6      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
##  7      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
##  8      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
##  9      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
## 10      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
## 11      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
## 12      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
## 13      15 BOYACA         15185 CHITARAQ~ OTROS PERMANENT~ CANA            
## # ... with 8 more variables: CULTIVO <chr>, YEAR <dbl>,
## #   `Area_Sembrada_(ha)` <dbl>, `Area_Cosechada_(ha)` <dbl>,
## #   `Produccion_(t)` <dbl>, `Rendimiento_(t/ha)` <dbl>,
## #   ESTADO_FISICO_PRODUCCION <chr>, CICLO_DE_CULTIVO <chr>

Hagamos una exploración gráfica básica:

g <- ggplot(aes(x=YEAR, y=`Produccion_(t)`/1000), data = chitaraque_cana) + geom_bar(stat='identity') + labs(y='Producción de caña [Ton x 1000]')
g + ggtitle("Evoluvión de la producción de caña en Chitaraque de 2007 a 2018") + labs(caption= "Basado en datos de EMA (DANE, 2018)")

En esta gráfica es posible observar que la mayor producción de caña en este municipio fue en el año 2007 con ~76000 toneladas, así mismo la producción más baja dió lugar en el año 2017 con ~30500 toneladas de caña.

Ahora, investiguemos qué cultivos tuvieron la mayor área cosechada en 2018:

datos %>% 
  filter(YEAR==2018) %>% 
  group_by(GRUPO_DE_CULTIVO) %>%
  summarize(sum_area_cosecha = sum(`Area_Cosechada_(ha)`, na.rm = TRUE)) %>%
     arrange(desc(sum_area_cosecha)) -> total_area_cosecha

total_area_cosecha
## # A tibble: 11 x 2
##    GRUPO_DE_CULTIVO                                 sum_area_cosecha
##    <chr>                                                       <dbl>
##  1 OTROS PERMANENTES                                           37166
##  2 TUBERCULOS Y PLATANOS                                       35849
##  3 FRUTALES                                                    15898
##  4 CEREALES                                                     9226
##  5 HORTALIZAS                                                   9100
##  6 LEGUMINOSAS                                                  8924
##  7 FLORES Y FOLLAJES                                              88
##  8 OTROS TRANSITORIOS                                             80
##  9 PLANTAS AROMATICAS, CONDIMENTARIAS Y MEDICINALES               48
## 10 FIBRAS                                                         35
## 11 OLEAGINOSAS                                                    25

Podemos ver que otros cultivos permanentes tuvieron la mayor parte del área cosechada en 2018 para Boyacá. ¿Qué cultivos pertenecen a este GRUPO? Puede buscar más información en DANE, la agencia gubernamental a cargo de las estadísticas nacionales en Colombia. También podemos buscar esta información en los datos mismos:

datos %>%
  filter(GRUPO_DE_CULTIVO=="OTROS PERMANENTES" & YEAR==2018) %>%
  group_by(CULTIVO) %>%
  summarize(sum_cosecha = sum(`Area_Cosechada_(ha)`, na.rm = TRUE)) %>%
     arrange(desc(sum_cosecha)) -> total_cosecha


total_cosecha
## # A tibble: 5 x 2
##   CULTIVO       sum_cosecha
##   <chr>               <dbl>
## 1 CANA PANELERA       19021
## 2 CAFE                 9653
## 3 CACAO                4840
## 4 CANA MIEL            3454
## 5 TABACO NEGRO          198

De los cultivos permantentes, el que tuvo la mayor área cosechada en el año 2018 fue la caña panelera, pues municipios como Buenavista, Coper, San Pablo de Borbur, Muzo, Briceño, Maripí, chitaraque y Miraflores, son destacados por ser productores de panela, posesionando al departamento como el tercer productor de Colombia.

Ahora, veamos cuáles son los municipios con mayor área cosechada para cada cultivo permanente en 2018:

datos %>% 
  filter(YEAR==2018 & GRUPO_DE_CULTIVO=="OTROS PERMANENTES") %>% 
  group_by(CULTIVO, MUNICIPIO) %>%
  summarize(max_area2 = max(`Area_Cosechada_(ha)`, na.rm = TRUE)) %>%
    slice(which.max(max_area2)) -> area_cosecha2

area_cosecha2
## # A tibble: 5 x 3
## # Groups:   CULTIVO [5]
##   CULTIVO       MUNICIPIO  max_area2
##   <chr>         <chr>          <dbl>
## 1 CACAO         MARIPI           780
## 2 CAFE          MONIQUIRA       1356
## 3 CANA MIEL     MARIPI           406
## 4 CANA PANELERA CHITARAQUE      7528
## 5 TABACO NEGRO  COVARACHIA       181

Teniendo en cuenta lo anterior, se verifica que CHITARAQUE se dedica preferencialmente al cultivo de la caña panelera, pues se considera el municipio con mayor área cosechada para este cultivo en el 2018. Por otra parte se resalta que el municipio de Maripi sobresale para los cultivos de cacao y caña miel.

Volvamos a los datos de cultivos permanentes. Antes de trazar, necesitaremos agregar, al objeto total_area_cosecha, un nuevo campo con abreviaturas para cada GRUPO de cultivos. De lo contrario, el plot se verá desordenado.

total_area_cosecha$CROP <- abbreviate(total_area_cosecha$GRUPO_DE_CULTIVO, 3)

Ahora sigue el plot:

g <- ggplot(aes(x=CROP, y=sum_area_cosecha), data = total_area_cosecha) + geom_bar(stat='identity') + labs(y='Área total cosechada [Ha]')
g+ ggtitle("Área total cosechada por grupos de cultivos en 2018 para Boyacá") + theme(plot.title = element_text(hjust = 0.5)) +
   labs(caption= "Basado en datos de EMA (DANE, 2018)")

Explorando más las estadísticas correspondientes a Boyacá. Intentaré obtener información sobre cuál es el municipio que presenta el rendimiento máximo de la papa en el año 2018, recordar que la papa es un cultivo escencial para la economía del departamento.

datos %>% 
  filter(YEAR=="2018", SUBGRUPO_DE_CULTIVO=="PAPA") %>% 
  group_by(CULTIVO, MUNICIPIO) %>%
  summarise(max_rend_18 = max(`Rendimiento_(t/ha)`, na.rm = TRUE)) %>%
  slice(which.max(max_rend_18)) -> Rend_max
Rend_max
## # A tibble: 1 x 3
## # Groups:   CULTIVO [1]
##   CULTIVO MUNICIPIO  max_rend_18
##   <chr>   <chr>            <dbl>
## 1 PAPA    PAZ DE RIO        3011

El municipio de Paz de Río para el año 2018 presentó el maximo rendimiento en el cultivo de papa con 3011(t/ha).

4. Unir las estadísticas agrícolas a los municipios.

En este caso usaremos el archivo que corresponde a Marco Geoestadistico Departamental que está disponible en el Geoportal DANE. Comenzaré leyendo los datos usando la biblioteca sf:

ant_munic <- sf::st_read("C:/Users/Brian/Desktop/Daniela/Boyacá/Marco Geoestadístico Boyacá/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\Brian\Desktop\Daniela\Boyacá\Marco Geoestadístico Boyacá\15_BOYACA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 123 features and 9 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

¿Qué es ant_munic?

ant_munic
## Simple feature collection with 123 features and 9 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
##    DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                           MPIO_CRSLC
## 1          15      15001      TUNJA                                 1541
## 2          15      15022    ALMEIDA                                 1908
## 3          15      15047  AQUITANIA                                 1789
## 4          15      15051   ARCABUCO                                 1856
## 5          15      15087      BELÉN                                 1756
## 6          15      15090     BERBEO      Ordenanza 28 de Abril 9 de 1913
## 7          15      15092  BETÉITIVA                                 1754
## 8          15      15097    BOAVITA                                 1613
## 9          15      15104     BOYACÁ                                 1537
## 10         15      15106    BRICEÑO Ordenanza 14 del 25 de Julio de 1890
##    MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng  Shape_Area
## 1   119.68957      2017     BOYACÁ  0.5723744 0.009766301
## 2    57.67212      2017     BOYACÁ  0.3484692 0.004701759
## 3   942.14660      2017     BOYACÁ  1.8003115 0.076843504
## 4   137.89859      2017     BOYACÁ  0.7527090 0.011256738
## 5   163.08822      2017     BOYACÁ  0.6293489 0.013314920
## 6    58.01301      2017     BOYACÁ  0.4291743 0.004730850
## 7   101.89955      2017     BOYACÁ  0.4738184 0.008317810
## 8   145.30529      2017     BOYACÁ  0.6597822 0.011867743
## 9    48.02287      2017     BOYACÁ  0.3256140 0.003918022
## 10   64.59970      2017     BOYACÁ  0.4849753 0.005273255
##                          geometry
## 1  POLYGON ((-73.34014 5.58308...
## 2  POLYGON ((-73.36793 5.01349...
## 3  POLYGON ((-72.76242 5.63856...
## 4  POLYGON ((-73.50487 5.84347...
## 5  POLYGON ((-72.91692 6.08612...
## 6  POLYGON ((-73.0677 5.27048,...
## 7  POLYGON ((-72.81796 5.97422...
## 8  POLYGON ((-72.64907 6.43640...
## 9  POLYGON ((-73.34806 5.47411...
## 10 POLYGON ((-73.89118 5.73749...

ant_munic es una colección de características simples. Las características simples o el acceso a características simples se refieren a un estándar formal (ISO 19125-1: 2004) que describe cómo los objetos en el mundo real se pueden representar en las computadoras, con énfasis en la geometría espacial de estos objetos. También describe cómo dichos objetos pueden almacenarse y recuperarse de las bases de datos, y qué operaciones geométricas deben definirse para ellos.

Por otra parte, los datos utilizan el sistema de referencia de coordenadas geográficas WGS1984 (es decir, código 4326 epsg). Podemos usar la función left_join para unir los municipios y las estadísticas agrícolas seleccionadas.

Necesitamos un atributo común (o variable compartida) para basar la unión. El mejor atributo es una identificación. En ant_munic, el atributo MPIO_CCDGO parece estar bien (se lee 15001 para Tunja). En datos, el atributo correspondiente es COD_MUN (lee 15001 para Tunja). Verifiquemos la última declaración:

datos %>% filter (MUNICIPIO =="TUNJA") ->  tunja_datos
tunja_datos
## # A tibble: 189 x 14
##    COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO_DE_CULTIVO SUBGRUPO_DE_CUL~
##      <dbl> <chr>          <dbl> <chr>     <chr>            <chr>           
##  1      15 BOYACA         15001 TUNJA     TUBERCULOS Y PL~ ARRACACHA       
##  2      15 BOYACA         15001 TUNJA     TUBERCULOS Y PL~ ARRACACHA       
##  3      15 BOYACA         15001 TUNJA     TUBERCULOS Y PL~ ARRACACHA       
##  4      15 BOYACA         15001 TUNJA     TUBERCULOS Y PL~ ARRACACHA       
##  5      15 BOYACA         15001 TUNJA     TUBERCULOS Y PL~ ARRACACHA       
##  6      15 BOYACA         15001 TUNJA     TUBERCULOS Y PL~ ARRACACHA       
##  7      15 BOYACA         15001 TUNJA     TUBERCULOS Y PL~ ARRACACHA       
##  8      15 BOYACA         15001 TUNJA     LEGUMINOSAS      ARVEJA          
##  9      15 BOYACA         15001 TUNJA     LEGUMINOSAS      ARVEJA          
## 10      15 BOYACA         15001 TUNJA     LEGUMINOSAS      ARVEJA          
## # ... with 179 more rows, and 8 more variables: CULTIVO <chr>, YEAR <dbl>,
## #   `Area_Sembrada_(ha)` <dbl>, `Area_Cosechada_(ha)` <dbl>,
## #   `Produccion_(t)` <dbl>, `Rendimiento_(t/ha)` <dbl>,
## #   ESTADO_FISICO_PRODUCCION <chr>, CICLO_DE_CULTIVO <chr>
class(tunja_datos$COD_MUN)
## [1] "numeric"

Para poder hacer la unión, necesitamos cambiar tanto el tipo de datos. Para esta tarea, es una buena idea crear una copia de los datos estadísticos originales. Con este enfoque, cualquier movimiento falso no estropeará los datos originales. Avancemos paso a paso:

datos2 <- datos
datos2$TEMP <- as.character(datos2$COD_MUN)
datos2$MPIO_CCDGO <- as.factor(paste(datos2$TEMP, sep=""))
head(datos2)
## # A tibble: 6 x 16
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO_DE_CULTIVO SUBGRUPO_DE_CUL~
##     <dbl> <chr>          <dbl> <chr>     <chr>            <chr>           
## 1      15 BOYACA         15114 BUSBANZA  HORTALIZAS       ACELGA          
## 2      15 BOYACA         15516 PAIPA     HORTALIZAS       ACELGA          
## 3      15 BOYACA         15491 NOBSA     HORTALIZAS       ACELGA          
## 4      15 BOYACA         15500 OICATA    HORTALIZAS       ACELGA          
## 5      15 BOYACA         15516 PAIPA     HORTALIZAS       ACELGA          
## 6      15 BOYACA         15500 OICATA    HORTALIZAS       ACELGA          
## # ... with 10 more variables: CULTIVO <chr>, YEAR <dbl>,
## #   `Area_Sembrada_(ha)` <dbl>, `Area_Cosechada_(ha)` <dbl>,
## #   `Produccion_(t)` <dbl>, `Rendimiento_(t/ha)` <dbl>,
## #   ESTADO_FISICO_PRODUCCION <chr>, CICLO_DE_CULTIVO <chr>, TEMP <chr>,
## #   MPIO_CCDGO <fct>

Ahora, filtremos un solo cultivo y seleccionemos solo dos atributos relevantes:

datos2 %>% filter(YEAR==2018, CULTIVO == "MAIZ")  -> datos3
head(datos3)
## # A tibble: 6 x 16
##   COD_DEP DEPARTAMENTO COD_MUN MUNICIPIO GRUPO_DE_CULTIVO SUBGRUPO_DE_CUL~
##     <dbl> <chr>          <dbl> <chr>     <chr>            <chr>           
## 1      15 BOYACA         15810 TIPACOQUE CEREALES         MAIZ            
## 2      15 BOYACA         15673 SAN MATEO CEREALES         MAIZ            
## 3      15 BOYACA         15660 SAN EDUA~ CEREALES         MAIZ            
## 4      15 BOYACA         15401 LA VICTO~ CEREALES         MAIZ            
## 5      15 BOYACA         15764 SORACA    CEREALES         MAIZ            
## 6      15 BOYACA         15469 MONIQUIRA CEREALES         MAIZ            
## # ... with 10 more variables: CULTIVO <chr>, YEAR <dbl>,
## #   `Area_Sembrada_(ha)` <dbl>, `Area_Cosechada_(ha)` <dbl>,
## #   `Produccion_(t)` <dbl>, `Rendimiento_(t/ha)` <dbl>,
## #   ESTADO_FISICO_PRODUCCION <chr>, CICLO_DE_CULTIVO <chr>, TEMP <chr>,
## #   MPIO_CCDGO <fct>
class(datos3)
## [1] "spec_tbl_df" "tbl_df"      "tbl"         "data.frame"
datos4 <- datos3 %>% dplyr::select(MUNICIPIO, MPIO_CCDGO, YEAR, `Produccion_(t)`, `Rendimiento_(t/ha)`) 
datos4
## # A tibble: 112 x 5
##    MUNICIPIO          MPIO_CCDGO  YEAR `Produccion_(t)` `Rendimiento_(t/ha)`
##    <chr>              <fct>      <dbl>            <dbl>                <dbl>
##  1 TIPACOQUE          15810       2018              198                   99
##  2 SAN MATEO          15673       2018              270                    2
##  3 SAN EDUARDO        15660       2018               19                   16
##  4 LA VICTORIA        15401       2018              440                    1
##  5 SORACA             15764       2018              300                   75
##  6 MONIQUIRA          15469       2018              260                    1
##  7 SUTATENZA          15778       2018              417                  183
##  8 TIPACOQUE          15810       2018              199                   99
##  9 SOATA              15753       2018              240                   15
## 10 SAN LUIS DE GACENO 15667       2018              330                  216
## # ... with 102 more rows
datos4 %>% 
  gather("YEAR", "Produccion_(t)", "Rendimiento_(t/ha)" , key = variable, value = number)
## # A tibble: 336 x 4
##    MUNICIPIO          MPIO_CCDGO variable number
##    <chr>              <fct>      <chr>     <dbl>
##  1 TIPACOQUE          15810      YEAR       2018
##  2 SAN MATEO          15673      YEAR       2018
##  3 SAN EDUARDO        15660      YEAR       2018
##  4 LA VICTORIA        15401      YEAR       2018
##  5 SORACA             15764      YEAR       2018
##  6 MONIQUIRA          15469      YEAR       2018
##  7 SUTATENZA          15778      YEAR       2018
##  8 TIPACOQUE          15810      YEAR       2018
##  9 SOATA              15753      YEAR       2018
## 10 SAN LUIS DE GACENO 15667      YEAR       2018
## # ... with 326 more rows
head(datos4)
## # A tibble: 6 x 5
##   MUNICIPIO   MPIO_CCDGO  YEAR `Produccion_(t)` `Rendimiento_(t/ha)`
##   <chr>       <fct>      <dbl>            <dbl>                <dbl>
## 1 TIPACOQUE   15810       2018              198                   99
## 2 SAN MATEO   15673       2018              270                    2
## 3 SAN EDUARDO 15660       2018               19                   16
## 4 LA VICTORIA 15401       2018              440                    1
## 5 SORACA      15764       2018              300                   75
## 6 MONIQUIRA   15469       2018              260                    1

Esta es una tarea clave. Implica varios pasos para poder convertir la tabla de atributos de formato largo a formato ancho.

datos4 %>% 
  group_by(MPIO_CCDGO) %>% 
  mutate(Visit = 1:n()) %>% 
  gather("YEAR", "Produccion_(t)", "Rendimiento_(t/ha)", key = variable, value = number) %>% 
  unite(combi, variable, Visit) %>%
  spread(combi, number) -> datos5
head(datos5)
## # A tibble: 6 x 8
## # Groups:   MPIO_CCDGO [123]
##   MUNICIPIO MPIO_CCDGO `Produccion_(t)~ `Produccion_(t)~ `Rendimiento_(t~
##   <chr>     <fct>                 <dbl>            <dbl>            <dbl>
## 1 ALMEIDA   15022                   114               NA              175
## 2 ARCABUCO  15051                    11               NA               14
## 3 BELEN     15087                    32               NA              202
## 4 BERBEO    15090                    17               NA                1
## 5 BETEITIVA 15092                    19               NA                8
## 6 BOAVITA   15097                   105               NA                7
## # ... with 3 more variables: `Rendimiento_(t/ha)_2` <dbl>, YEAR_1 <dbl>,
## #   YEAR_2 <dbl>
tail(datos5)
## # A tibble: 6 x 8
## # Groups:   MPIO_CCDGO [123]
##   MUNICIPIO MPIO_CCDGO `Produccion_(t)~ `Produccion_(t)~ `Rendimiento_(t~
##   <chr>     <fct>                 <dbl>            <dbl>            <dbl>
## 1 TURMEQUE  15835                   120               NA               15
## 2 TUTA      15837                    40               NA                5
## 3 TUTAZA    15839                     3               NA               13
## 4 UMBITA    15842                    93               NA              185
## 5 VILLA DE~ 15407                    30               NA                1
## 6 VIRACACHA 15879                    89               NA              124
## # ... with 3 more variables: `Rendimiento_(t/ha)_2` <dbl>, YEAR_1 <dbl>,
## #   YEAR_2 <dbl>

También haremos una copia de la colección de características simples (de nuevo, solo en caso de un movimiento falso):

ant_munic2 <- ant_munic

Ahora. sigue la unión:

ant_munic_stat = left_join(ant_munic2, datos5, by="MPIO_CCDGO")
summary(ant_munic_stat)
##  DPTO_CCDGO   MPIO_CCDGO      MPIO_CNMBR 
##  15:123     15001  :  1   ALMEIDA  :  1  
##             15022  :  1   AQUITANIA:  1  
##             15047  :  1   ARCABUCO :  1  
##             15051  :  1   BELÉN    :  1  
##             15087  :  1   BERBEO   :  1  
##             15090  :  1   BETÉITIVA:  1  
##             (Other):117   (Other)  :117  
##                                 MPIO_CRSLC   MPIO_NAREA        MPIO_NANO   
##  1754                                :10   Min.   :  25.35   Min.   :2017  
##  1777                                : 6   1st Qu.:  64.26   1st Qu.:2017  
##  Ordenanza 41 de Diciembre 14 de 1978: 6   Median : 122.70   Median :2017  
##  1799                                : 4   Mean   : 188.11   Mean   :2017  
##  1776                                : 3   3rd Qu.: 203.16   3rd Qu.:2017  
##  1780                                : 3   Max.   :1513.60   Max.   :2017  
##  (Other)                             :91                                   
##   DPTO_CNMBR    Shape_Leng       Shape_Area        MUNICIPIO        
##  BOYACÁ:123   Min.   :0.2045   Min.   :0.002069   Length:123        
##               1st Qu.:0.3858   1st Qu.:0.005245   Class :character  
##               Median :0.5350   Median :0.010017   Mode  :character  
##               Mean   :0.6365   Mean   :0.015353                     
##               3rd Qu.:0.7505   3rd Qu.:0.016568                     
##               Max.   :2.4039   Max.   :0.123610                     
##                                                                     
##  Produccion_(t)_1 Produccion_(t)_2 Rendimiento_(t/ha)_1 Rendimiento_(t/ha)_2
##  Min.   :  3.00   Min.   : 21.00   Min.   :  1.00       Min.   : 13.0       
##  1st Qu.: 24.00   1st Qu.: 34.50   1st Qu.:  7.00       1st Qu.: 56.0       
##  Median : 50.00   Median : 48.00   Median : 18.00       Median : 99.0       
##  Mean   : 81.94   Mean   : 89.33   Mean   : 64.61       Mean   : 72.0       
##  3rd Qu.:105.00   3rd Qu.:123.50   3rd Qu.:116.00       3rd Qu.:101.5       
##  Max.   :440.00   Max.   :199.00   Max.   :257.00       Max.   :104.0       
##  NA's   :14       NA's   :120      NA's   :14           NA's   :120         
##      YEAR_1         YEAR_2              geometry  
##  Min.   :2018   Min.   :2018   POLYGON      :123  
##  1st Qu.:2018   1st Qu.:2018   epsg:4326    :  0  
##  Median :2018   Median :2018   +proj=long...:  0  
##  Mean   :2018   Mean   :2018                      
##  3rd Qu.:2018   3rd Qu.:2018                      
##  Max.   :2018   Max.   :2018                      
##  NA's   :14     NA's   :120

5. Plotting

#install.packages("RColorBrewer")
library(RColorBrewer)
library(leaflet)

Tracemos los municipios con su producción de maíz correspondiente para un solo año:

bins <- c(0, 50, 100, 150, 200, 300, 400, 600)
pal <- colorBin("YlOrRd", domain = ant_munic_stat$`Produccion_(t)_1`, bins = bins)

  mapa <- leaflet(data = ant_munic_stat) %>%
  addTiles() %>%
  addPolygons(label = ~`Produccion_(t)_1`,
              popup = ~MPIO_CNMBR,
              fillColor = ~pal(`Produccion_(t)_1`),
              color = "#444444",
              weight = 1,
              smoothFactor = 0.5,
              opacity = 1.0,
              fillOpacity = 0.5,
              highlightOptions = highlightOptions(color = "white", weight = 2, bringToFront = TRUE)
              ) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend("topleft", pal = pal, values = ~`Produccion_(t)_1`,
    title = "Producción de maíz en Boyacá [Ton] (2018)",
    opacity = 2
  )
mapa

Observando el mapa, se obtiene que el municipio con mayor producción de maíz es la victoria seguido de Sutatenza con mas de 400 toneladas en el año 2018.

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252  LC_CTYPE=Spanish_Colombia.1252   
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] leaflet_2.0.3       RColorBrewer_1.1-2  cowplot_1.0.0      
##  [4] ggrepel_0.8.2       GSODR_2.0.1         rnaturalearth_0.1.0
##  [7] viridis_0.5.1       viridisLite_0.3.0   sf_0.8-1           
## [10] raster_3.0-12       maptools_0.9-9      rgeos_0.5-2        
## [13] sp_1.4-1            forcats_0.5.0       stringr_1.4.0      
## [16] dplyr_0.8.5         purrr_0.3.3         readr_1.3.1        
## [19] tidyr_1.0.2         tibble_2.1.3        ggplot2_3.3.0      
## [22] tidyverse_1.3.0     here_0.1           
## 
## loaded via a namespace (and not attached):
##  [1] httr_1.4.1              jsonlite_1.6.1          modelr_0.1.6           
##  [4] assertthat_0.2.1        cellranger_1.1.0        yaml_2.2.1             
##  [7] pillar_1.4.3            backports_1.1.5         lattice_0.20-38        
## [10] glue_1.3.2              digest_0.6.25           rvest_0.3.5            
## [13] leaflet.providers_1.9.0 colorspace_1.4-1        htmltools_0.4.0        
## [16] pkgconfig_2.0.3         broom_0.5.5             haven_2.2.0            
## [19] scales_1.1.0            farver_2.0.3            generics_0.0.2         
## [22] ellipsis_0.3.0          withr_2.1.2             cli_2.0.2              
## [25] magrittr_1.5            crayon_1.3.4            readxl_1.3.1           
## [28] evaluate_0.14           fs_1.3.2                fansi_0.4.1            
## [31] nlme_3.1-144            xml2_1.2.5              foreign_0.8-75         
## [34] class_7.3-15            tools_3.6.3             data.table_1.12.8      
## [37] hms_0.5.3               lifecycle_0.2.0         munsell_0.5.0          
## [40] reprex_0.3.0            compiler_3.6.3          e1071_1.7-3            
## [43] rlang_0.4.5             classInt_0.4-2          units_0.6-5            
## [46] grid_3.6.3              rstudioapi_0.11         htmlwidgets_1.5.1      
## [49] crosstalk_1.1.0.1       labeling_0.3            rmarkdown_2.1          
## [52] gtable_0.3.0            codetools_0.2-16        DBI_1.1.0              
## [55] R6_2.4.1                gridExtra_2.3           lubridate_1.7.4        
## [58] knitr_1.28              utf8_1.1.4              rprojroot_1.3-2        
## [61] KernSmooth_2.23-16      stringi_1.4.6           Rcpp_1.0.3             
## [64] vctrs_0.2.4             dbplyr_1.4.2            tidyselect_1.0.0       
## [67] xfun_0.12