1. ¿Por qué este cuaderno?

Este es un cuaderno realizado en RStudio Cloud que ilustra la cartografía temática para el departamento del Valle del Cauca y tiene como objetivo ayudar a los estudiantes de la Universidad Nacional de Colombia que se encuentran cursando geomática básica a aprender como realizar un mapa temático en R

2. Cartografía temática

Según el Instituto Geográfico Agustín Codazzi (2018) los mapas temáticos son elaborados con un propósito especial, según su contenido se pueden clasificar en: geológicos, catastrales, turísticos, de suelos, entre otros. Los mapas temáticos requieren para su elaboración de un mapa topográfico como base.

Paa Olaya, la cartografía temática se centra en la representación de un tema concreto (una variable espacial dada), de cualquier índole: física, social, política, cultural, etc. Se excluyen de la lista de esos temas posibles a los puramente topográficos, que constituyen el objeto de la cartografía base. Las entidades encargadas de realizar mapas temáticos en Colombia son entidades como el Ministerio de Agricultura, IDEAM, DANE, etc.

Un mapa temático se compone, así pues, de dos partes bien diferenciadas:

2.1 Tipo de mapas cartográficos

Existen diferentes alternativas para representar un mapa, por lo que resulta pertienente entender la forma en la que se presenta. A continuación se explicará brevemente los mapas de símbolos proporcionales, mapas de puntos, mapas de isolineas y mapas de coropletas.

- Mapas de símbolos proporcionales

Representa variables cuantitativas mediante símbolos cuyo tamaño está en relación con el valor a representar de dicha variable. Estos mapas utilizan la variable visual tamaño, la cual es la única que presenta la propiedad cuantitativa. La forma de los diferentes símbolos es siempre la misma, y la más frecuente por cuestiones de simplicidad es el círculo, aunque se puede usar otros como lineas. Para estos mapas el tamaño es el elemento que diferencia a los distintos símbolos y el que transmite la información cuantitativa. La elección de un tamaño implica elegir uno mínimo y uno máximo, correspondientes a los valores mínimo y máximo de la variable en el mapa. Entre estos se situarán los distintos tamaños correspondientes al resto de posible valores que toma la variable

- Mapas de puntos

Se emplean para la representación de variables que simbolicen algún tipo de cantidad, como la población o producción de algún cultivo en particular. Estas cantidades se representan a través de la repetición de puntos en un número proporcional a su magnitud. Cada uno de esos puntos representa un valor unitario, y el conjunto de ellos sobre la zona en cuestión suma la cantidad total a representar. Los puntos tienen todos la misma forma y tamaño, a diferencia del tipo de mapa anterior en donde sí varía el tamaño de los símbolos.

- Mapas de isolíneas

Son unos de los más usados para la representación de la información cuantitativa. Son utilizados comunmente para representar campos escalares. Un mapa de isolíneas está formado por un conjunto de líneas, cada una de las cuales une puntos que presentan el mismo valor de la variable. Estas líneas no pueden cruzarse, ya que ello significaría que en un punto se presentan dos valores. El caso más típico de mapa de isolíneas son las curvas de nivel que aparecen el un mapa topográfico, indicando la elevación del terreno

- Mapas de coropletas

Habitualmente son utilizados para representar la información geográfica en un SIG. En un mapa de coropletas se tiene una serie de áreas definidas, cada una de las cuales tiene un valor de una variable. Este valor, afecta a todo el área y es el que se representa por medio de alguna variable visual, normalmente el color a través de su componente valor.

3. Datos

Para llevar a cabo este cuaderno de R se utilizarán datos de las Necesidades Básicas Insatisfechas del Censo Nacional de Población y Vivienda, los cuales están disponibles siguiendo este link del geoportal oficial del DANE https://www.dane.gov.co/index.php/estadisticas-por-tema/pobreza-y-condiciones-de-vida/necesidades-basicas-insatisfechas-nbi

Este tipo de indicador que mide la condición de pobreza estructural se expresa como el porcentaje de personas vs. hogares sobre la población total vs. total de hogares que tiene al menos una necesidad básica insatisfecha (NBI). En Colombia y según el DANE, se tienen en cuenta las siguientes NBI: viviendas con hacinamiento crítico, con condiciones físicas impropias para el alojamiento humano, servicios inadecuados, alta dependencia económica o niños en edad escolar que no asisten a la escuela.

Se descargó el archivo de Necesidades Basicas Insatisfechas en formato xlsx del geoportal del Dane, posteriormente se usó Excel para remover aquellas columnas con imagenes institucionales o sin información. También se eliminó información que no se viera relacionada con el departamento de estudio, es decir, del Valle del Cauca e información de necesidades basicas insatisfechas correspondientes a zonas de cabecera y rural. Después se cambio el título de algunas columnas que tenían espacios por guiones para que el programa R los leyera con facilidad y se concatenaron los atributos codigo de departamento con código de municipio. Finalmente se subío el archivo xlsx a RStudio Cloud como NBI_Valle.xlsx

4. Preparación

Removamos todas las variables que hay en la memoria para partir desde 0.

rm(list=ls())

A continuación instalamos las librerías que vamos a utilizar

list.of.packages <- c("tidyverse", "rgeos", "sf", "raster", "cartography", "SpatialPosition")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

Cargamos las librerías

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.3     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(readxl) # la cual es parte de tidyverse
library(rgeos)
## Loading required package: sp
## rgeos version: 0.5-3, (SVN revision 634)
##  GEOS runtime version: 3.5.1-CAPI-1.9.1 
##  Linking to sp version: 1.4-1 
##  Polygon checking: 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.5.1, GDAL 2.2.2, PROJ 4.9.2
library(cartography)
library(SpatialPosition)

5. Leer datos de Necesidades Básicas Insatisfechas

Leamos el archivo xlsx de Necesidades Básicas Insatisfechas para el departamento del Valle del Cauca.

nbi <- read_excel("./NBI_Valle.xlsx")

Miremos cuales son los atributos de los datos de NBI, viendo las primeras 6 filas del archivo:

head(nbi)

Ahora, miremos cual es el municipio del Valle del Cauca que posee el porcentaje más alto de necesidades básicas insatisfechas por medio del codigo que se mostrará a continiación el cual selecciona aquella fila en la cual el valor de NBI es el máximo.

nbi %>% 
    slice(which.max(NBI)) -> max_nbi

max_nbi

Según los registros del DANE, la ciudad con mayor porcentaje de NBI es Buenaventura, con un porcentaje de 16.57241.

Esto se debe a que Buenaventura, el puerto más importante de Colombia, ubicado al occidente del departamento, presenta altos indices de pobreza, las condiciones de salud de la ciudad son perversas, pues se cuenta con un hospital de tercel nivel el cual no alcanza para cubrira toda la población, los habitantes de la ciudad expresan que “no hay cama pa’ tanta gente”. Además no hay agua potable y la tasa de desempleo sobrepasa el 62%.

Miremos cúal es la ciudad del departamento con el porcentaje más bajo de NBI.

nbi %>% 
    slice(which.min(NBI)) -> min_nbi

min_nbi

Según el DANE, la ciudad del departamento del Valle del Cauca con menor porcentaje de NBI es Guadalajara de Buga. En Buga el 98,9% de las viviendas tiene conexión a Energía eléctrica . El 49,3% tiene conexión a Gas Natural, el 83,6% de las viviendas de BUGA son casas. Los anteriores reportes dan argumento al bajo porcentaje de NBI en esta ciudad.

Ahora, ordenaremos los municipios en función del NBI de manera descendente con la funcion arrange que ordenan las filas de un data frame en función de los valores de una o más columnas.

nbi %>% 
  arrange(desc(NBI))  -> desc_nbi

desc_nbi

De igual manera miraremos cúal es el municipio que reporta mayor hacinamiento en el departamento:

hacinamiento <- read_excel("./NBI_Valle.xlsx")
head(hacinamiento)
hacinamiento %>% 
    slice(which.max(HACINAMIENTO)) -> max_hacinamiento

max_hacinamiento

Según el DANE el municipio con mayor hacinamiento es El Dovio, municipio ubicado en el norte del Valle.

Miraremos cúal es el municipio con menor hacinamiento.

hacinamiento %>% 
    slice(which.min(HACINAMIENTO)) -> min_hacinamiento

min_hacinamiento

El municipio que reporta menor hacinamiento es Palmira, en donde el número de personas por hogar es de 3,7 (DANE).

6. Uniendo datos de Necesidades Básicas Insatisfechas a municipios.

Se suben los datos de los municipios del Valle del Cauca a RStudio Cloud. Leeremos los datos usando la libreía st

munic <- st_read("./ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source `/cloud/project/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 42 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
## CRS:            4326

Ahora miremos lo que se encuentra dentro del atributo MPIO_CCDGO:

head(munic$MPIO_CCDGO)
## [1] 76001 76036 76041 76054 76111 76113
## 42 Levels: 76001 76020 76036 76041 76054 76100 76109 76111 76113 ... 76895

Podemos usar la función left_join para unir los datos de los municipios con los de NBI.

nbi_munic = left_join(munic, nbi, by=c("MPIO_CCDGO"="CODIGO"))
## Warning: Column `MPIO_CCDGO`/`CODIGO` joining factor and character vector,
## coercing into character vector
nbi_munic %>%
  dplyr::select(MUNICIPIO, MPIO_CCDGO, NBI)  ->  check_nbi_munic

head(check_nbi_munic)

Esto se hace con el fin de verificar que las operaciones conducen a datos que tienen valores. Notemos que a Gudalajara de Buga, la ciudad que analizamos previamente, le corresponde el valor de NBI de 4.01, lo cual se puede verificar con los datos que se analizaron de NBI.

Ahora lo que haremos es reproyectar los municipios, teniendo en cuenta que la libreria cartography requiere coordenadas planas y los datos originales estan en coordenadas geograficas, por lo tanto se realiza una una transformación mediante en siguiente código:

nbi_munic_new <- st_transform(nbi_munic, crs = 3116)

El sistema de referencia de coordenadas que se usa para el caso de Colombia es el sistema MAGNA-SIRGAS / Colombia Bogota zone y el codigo EPSG es el 3116

7. Ejemplos de mapas temáticos

Se usará el paquete cartography debido a que la calidad de los mapas es alta, se puede incluso comparar con la calidad de mapas que brindan softwares como QGIS. Este paquete usa objetos sf o sp para producir gráficos base. Como la mayoría de los componentes internos del paquete se basan en funcionalidades sf, el formato preferido para los objetos espaciales es sf.

7.1. Mapa de simbolos profesionales usando como mapa básico OpenStreetMap.

Primero se descarga el paquete osmdata y se lee la librería con el objetivo de trabajar con el mapa de OpenStreetMap.

Las funciones getTiles() y tilesLayer() descargan y muestran mosaicos de OpenStreetMap. propSymbolsLayer() muestra símbolos con áreas proporcionales a una variable cuantitativa (como Necesidades Básicas Insatisfechas). Símbolos como cuadrados, barras y círculos están disponibles y el elemento inches se utiliza para personalizar los tamaños de estos símbolos.

Miremos como se distribuye el hacinamiento en el Valle del Cauca. Con este indicador se busca captar los niveles críticos de ocupación de los recursos de la vivienda por el grupo que la habita. Se consideran en esta situación las viviendas con más de tres personas por cuarto ( excluyendo cocina, baño y garaje).

mun.osm <- getTiles(
x = nbi_munic_new, 
type = "OpenStreetMap", 
zoom = 8,
cachedir = TRUE,
crop = FALSE
)
# se establecen márgenes
opar <- par(mar = c(1.2,1.2,1.2,1.2))
# trazar azulejos osm
tilesLayer(x = mun.osm)
# trazar unicamente las fronteras de los municipios
plot(st_geometry(nbi_munic_new), col = NA, border = "black", add=TRUE)
# plot HACINAMIENTO
propSymbolsLayer(
  x = nbi_munic_new, 
  var = "HACINAMIENTO", 
  inches = 0.12, 
  col = "deepskyblue",
  legend.pos = "topright",  
  legend.title.txt = "Hacinamiento"
)
# diseño
layoutLayer(title = " Distribución del hacinamiento en el Valle del Cauca",
            sources = "Fuente: DANE, 2018\n© OpenStreetMap",
            author = " Nicolás Cifuentes",
            frame = TRUE, north = FALSE, tabtitle = TRUE)
# flecha norte
north(pos = "topleft")

Ahora miremos como se distribuye la miseria en el Valle del Cauca

# se establecen márgenes
opar <- par(mar = c(1.2,1.2,1.2,1.2))
# trazar azulejos osm
tilesLayer(x = mun.osm)
# trazar unicamente las fronteras de los municipios
plot(st_geometry(nbi_munic_new), col = NA, border = "black", add=TRUE)
# plot HACINAMIENTO
propSymbolsLayer(
  x = nbi_munic_new, 
  var = "MISERIA", 
  inches = 0.2, 
  col = "red",
  legend.pos = "topright",  
  legend.title.txt = "Miseria"
)
# diseño
layoutLayer(title = " Distribución de la miseria en el Valle del Cauca",
            sources = "Fuente: DANE, 2018\n© OpenStreetMap",
            author = " Nicolás Cifuentes",
            frame = TRUE, north = FALSE, tabtitle = TRUE)
# flecha norte
north(pos = "topleft")

Se puede ver que los municipios con mayor porcentaje de miseria, en su mayoeía, se encuentran unicados al norte del país

Comprobemos la afirmación anterior ordenando los municipios de forma descendente con función a la miseria.

miseria <- read_excel("./NBI_Valle.xlsx")
miseria %>% 
  arrange(desc(MISERIA))  -> desc_miseria

desc_miseria

Observamos que efectivamente los municipios ubicados al norte del departamento como Bolívar, El Dovio, Trujillo, Argelia, El Cairo, Ansermanuevo, Toro, Versalles y Alcalá, presentan altos valores de miseria según el reporte.

7.2. Mapas de Coropletas

choroLayer() muestra mapas de coropletas. Los argumentos nclass, method y breaks permiten personalizar la clasificación variable

getBreaks() permite clasificar fuera de la función misma. Las paletas de colores se definen con color y se puede crear un conjunto de colores con carto.pal().

Miremos como se distribuyen las NBI en los municipios del Valle del Cauca, teniendo en cuenta que el tono más oscuro corresponde a una variable de mayor porcentaje de NBI y el tono más claro corresponde a una variable de menor porcentaje de NBI.

# se establecen márgenes
opar <- par(mar = c(1.2,1.2,1.2,1.2))
# establecer el color de fondo de la figura
par(bg="red")
# trazar municipios
plot(st_geometry(nbi_munic_new), col = NA, border = NA, bg = "#aadaff")
# plot NBI
choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "red.pal", n1 = 5),
  border = "gray4", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "NBI",
  add = TRUE
) 
# diseño
layoutLayer(title = "Distribución de NBI en el Valle del Cauca ", 
            sources = "Fuente: DANE, 2018",
            author = "Nicolás Cifuentes", 
            frame = TRUE, north = TRUE, tabtitle = TRUE, col="black") 
# north arrow
north(pos = "topleft")

Se puede comprobar que el municipio con mayor porcentaje de NBI es Buenaventura. Esto coincide con lo que estudiamos anteriormente.

7.3. Símbolos proporcionales en un mapa de tipología.

La función propSymbolsTypoLayer() crea un mapa de símbolos que son proporcionales a los valores de una primera variable y coloreados para reflejar las modalidades de una segunda variable cualitativa. Se utiliza una combinación de argumentos propSymbolsLayer() y typoLayer().

Primero necesitamos crear una variable cualitativa. Para esta tarea, usaremos la funcion mutate.

nbi_munic_2 <- dplyr::mutate(nbi_munic_new, poverty = ifelse(MISERIA > 5, "Extreme", 
                                                         ifelse(HACINAMIENTO > 3, "High", "Intermediate")))

Tener en cuenta que el nuevo atributo se llama pobreza y su valor depende del umbral de los datos de MISERIA y HACINAMIENTO.

head(nbi_munic_2)
library(sf)
library(cartography)
# establecer márgenes
opar <- par(mar = c(0,0,1.2,0))
# trazar los municipios
plot(st_geometry(nbi_munic_2), col="#f2efe9", border="#b38e43", bg = "#aad3df", 
     lwd = 0.5)
# trazar símbolos con coloración coropleta
propSymbolsTypoLayer(
  x = nbi_munic_2, 
  var = "NBI", 
  inches = 0.3,
  symbols = "circle",
  border = "white",
  lwd = .5,
  legend.var.pos = "topright", 
  legend.var.title.txt = "NBI",
  var2 = "poverty",
  legend.var2.values.order = c("Extreme", "High", "Intermediate"),
  col = carto.pal(pal1 = "multi.pal", n1 = 3),
  legend.var2.pos = "top",
  legend.var2.title.txt = "Pobreza"
) 
# diseño
layoutLayer(title="Distribución de NBI en el Valle del Cauca", 
            author = "Nicolas Cifuentes", 
            sources = "Source: DANE, 2018", 
            scale = 1, tabtitle = TRUE, frame = TRUE)
# flecha nortte
north(pos = "topleft")

7.4 Mapas de etiqueta

Intentemos combinar las funciones de choroLayer() y labelLayer:

Miremos como se distribuye en factor economía, el cual es un indicador indirecto sobre los niveles de ingreso. Se clasifican aquí, las viviendas en los cuales haya más de tres personas por miembro ocupado y el jefe tenga, como máximo, dos años de educación primaria aprobados.

library(sf)
library(cartography)
# establecer márgenes
opar <- par(mar = c(0,0,1.2,0))
# establecer color de fondo de figura
par(bg="grey25")
# trazar municipios
plot(st_geometry(nbi_munic_2), col = "#e4e9de", border = "darkseagreen4", 
     bg = "grey75", lwd = 0.5)
# plot NBI
choroLayer(
  x = nbi_munic_new, 
  var = "ECONOMIA",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "red.pal", n1 = 5),
  border = "black", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "NBI",
  add = TRUE
) 
# trazar etiquetas
labelLayer(
  x = nbi_munic_2, 
  txt = "MUNICIPIO", 
  col= "black", 
  cex = 0.5, 
  font = 4,
  halo = TRUE, 
  bg = "gray", 
  r = 0.2, 
  overlap = FALSE, 
  show.lines = FALSE
)
# diseño de mapa
layoutLayer(
  title = "Municipios del Valle del Cauca", 
  sources = "Fuente: DANE, 2018",  
  author = " Nicolas Cifuentes", 
  frame = TRUE,
  north = TRUE, 
  tabtitle = TRUE, 
  theme = "taupe.minimal"
)

Los municipios con mayor dependencia económica se encuentran en el norte del departamento.

7.5 Mapas de Isolíneas

Permiten una representación espacial del fenómeno independiente de la heterogeneidad inicial de la división territorial smoothLayer() depende en gran medida del paquete SpatialPosition. La función utiliza una capa de puntos marcados y un conjunto de parámetros (una función de interacción espacial y sus parámetros) y muestra una capa de mapa isopleta. Usaremos otro conjunto de datos para hacer un mapa isopleta. Se subirán otro conjunto de datos para hacer un mapa ispleta, en este caso se subirán datos estadísticos sobre la producción de café en el Valle del Cauca, ya que es un cultivo importante para la economía del departamento.

Se lee el conjunto de datos.

crops2018 <- read_excel("./valle/Evaluaciones_Valle.xlsx")

Se analiza que hay en las primeras filas de los datos

head(crops2018)

Filtraremos las filas para obtener unicamente datos que corresponden a Café

crops2018 %>%
  filter(CULTIVO == "CAFE") -> cafe2018

Se cambiará el tipo de dato que tiene COD_MUN, el cual es el codigo del municipio, donde inicialmente era de tipo numerico, pero con la funcion que se ejecutará a continuacion pasara a ser de tipo caracter, con el fin de poder realizar correctamdente el join.

cafe2018$TEMP <-  as.character(cafe2018$COD_MUN)
cafe2018$MPIO_CCDGO <- as.factor(paste(cafe2018$TEMP, sep=""))

Realicemos el join usando la función left_join

cafe_munic = left_join(munic, cafe2018, by="MPIO_CCDGO")
## Warning: Column `MPIO_CCDGO` joining factors with different levels, coercing to
## character vector

Miremos los datos que arroja

head(cafe_munic)

Ahora vamos a reproyectar los municipios

rep_cafe <- st_transform(cafe_munic, crs = 3116)

A mapear:

NO SE PUDO MAPEAR DEBIDO A QUE EL PAQUETE “LWGEON” NO ESTABA DISPONIBLE PARA RStudio Cloud

8. Guardando mapas

Realizaremos otro mapa de producción de café en 2018. Esta vez usaremos símbolos proporcionales y mapas de coropletas. La salida se guardará como un archivo .png.

Primero, algunas observaciones: propSymbolsChoroLayer() crea un mapa de símbolos que son proporcionales a los valores de una primera variable y coloreados para reflejar la clasificación de una segunda variable. Se utiliza una combinación de argumentos propSymbolsLayer() y choroLayer().

El siguiente fragmento no muestra un mapa. En cambio, escribe el mapa en el nombre de archivo cafe_2018.png debajo del directorio de trabajo.

### open the plot
png("./valle/cafe_2018.png", width = 2048, height = 1526)
# set margins
opar <- par(mar = c(0,0,5,5))
# Plot the municipalities
plot(st_geometry(rep_cafe), col="darkseagreen3", border="darkseagreen4",  
     bg = "grey", lwd = 0.6)
# Plot symbols with choropleth coloration
propSymbolsChoroLayer(x = rep_cafe, var = "Produccion", var2 = "Rendimiento",
                      col = carto.pal(pal1 = "green.pal", n1 = 3,
                                      pal2 = "red.pal", n2 = 3),
                      inches = 0.8, method = "q6",
                      border = "grey50", lwd = 1,
                      legend.title.cex = 1.5,
                      legend.values.cex = 1.0,
                      legend.var.pos = "right", 
                      legend.var2.pos = "left",
                      legend.var2.values.rnd = 2,
                      legend.var2.title.txt = "Rendimiento\n(in Ton/Ha)",
                      legend.var.title.txt = "Producción de café en 2018",
                      legend.var.style = "e")
# plot labels
labelLayer(
  x = rep_cafe, 
  txt = "MPIO_CNMBR", 
  col= "white", 
  cex = 1.0, 
  font = 4,
  halo = FALSE, 
  bg = "white", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
# layout
layoutLayer(title="Producción y rendimiento de café en el Valle del Cauca, 2018",
            author = "Nicolas Cifuentes", 
            sources = "FUente MADR & DANE, 2018", 
            scale = 50, tabtitle = FALSE, frame = TRUE)
# north arrow
north(pos = "topleft")
#
title(main="Producción y rendimiento de café en el Valle del Cauca, 2018", cex.main=3,
      sub= "Fuente: MADR & DANE, 2018", cex.sub=2)
#
graticule = TRUE
#
par(opar)
### close the plot
dev.off()
## png 
##   2

Una vez que guardamos un mapa, podemos agregarlo a nuestro informe R Markdown usando la sintaxis de Markdown de la siguiente manera

image:

Se puede observar que la mayor producción de café en el valle para el año dee estudio se ubica en el Norte del departamento, donde se encuentran grandes extensiones de café e importantes ciudades productoras.

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/atlas-base/atlas/libblas.so.3.0
## LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] SpatialPosition_2.0.1 cartography_2.4.1     sf_0.9-3             
##  [4] raster_3.1-5          rgeos_0.5-3           sp_1.4-1             
##  [7] readxl_1.3.1          forcats_0.5.0         stringr_1.4.0        
## [10] dplyr_0.8.5           purrr_0.3.4           readr_1.3.1          
## [13] tidyr_1.0.3           tibble_3.0.1          ggplot2_3.3.0        
## [16] tidyverse_1.3.0      
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.4.6       lubridate_1.7.8    lattice_0.20-38    png_0.1-7         
##  [5] class_7.3-15       assertthat_0.2.1   digest_0.6.25      R6_2.4.1          
##  [9] cellranger_1.1.0   backports_1.1.7    reprex_0.3.0       evaluate_0.14     
## [13] e1071_1.7-3        httr_1.4.1         pillar_1.4.4       rlang_0.4.6       
## [17] rstudioapi_0.11    rmarkdown_2.1      rgdal_1.4-8        munsell_0.5.0     
## [21] broom_0.5.6        compiler_3.6.0     modelr_0.1.7       xfun_0.13         
## [25] pkgconfig_2.0.3    htmltools_0.4.0    tidyselect_1.1.0   codetools_0.2-16  
## [29] fansi_0.4.1        crayon_1.3.4       dbplyr_1.4.3       withr_2.2.0       
## [33] grid_3.6.0         nlme_3.1-139       jsonlite_1.6.1     gtable_0.3.0      
## [37] lifecycle_0.2.0    DBI_1.1.0          magrittr_1.5       units_0.6-6       
## [41] scales_1.1.1       KernSmooth_2.23-15 cli_2.0.2          stringi_1.4.6     
## [45] fs_1.4.1           xml2_1.3.2         ellipsis_0.3.1     generics_0.0.2    
## [49] vctrs_0.3.0        tools_3.6.0        glue_1.4.1         hms_0.5.3         
## [53] slippymath_0.3.1   yaml_2.2.1         colorspace_1.4-1   isoband_0.2.1     
## [57] classInt_0.4-3     rvest_0.3.5        knitr_1.28         haven_2.2.0