Cartografía temática para el departamento de Antioquia

1. INTRODUCCIÓN

El presente cuaderno de R Markdown tiene como fin ilustrar la cartografía temática del departamento Antioquia, Colombia.

La cartografía temática representa cualquier fenómeno geográfico con variabiliad espacial y permite enfatizar en un tema particular, a diferencia de los mapas de referencia generales, no solo muestra características naturales como vegetación, suelos, geología, sino también conceptos abstractos como indicadores de violencia, calidad de vida, entre otros. Estos mapas están compuestos por dos elementos fundamentales, una base geográfica que proporciona información espacial para referenciar el contenido de un tema específico y una capa de contenido temático.

Un factor a tener en cuenta al diseñar este tipo de mapas es el público, ya que determina que elementos se deben incluir en el mapa como puntos de referencia además del tema.

Entre algunas técnicas de mapeo temático encontramos:

• Coropleta

Son una forma de cartografiado cuantitativo utilizada para la representación de datos de naturaleza discreta y que ocurren dentro de zonas bien definidas a las que se aplican símbolos superficiales como sombreados de acuerdo con su valor. Los datos son clasificados generalmente mediante intervalos, lo que permite tener una información general y una fácil lectura y análisis de los datos (Gutiérrez, J).

• Mapa de símbolos proporcionales

Es un tipo de mapa muy empleado en la cartografía temática por su fácil interpretación, se basa principalmente en seleccionar una forma como círculos, cuadrados, triángulos e ir variando su tamaño en proporción a las cantidades que se tengan que representar (Gutiérrez, J).

• Mapa de isolíneas

Una isolínea es una línea con un valor constante asociado a todos sus puntos como los niveles de precipitación, así como elevación y mapas topográficos. Este tipo de mapa proporciona una visión global de la configuración de la superficie estadística y muestra la distribución de la variación espacial de un fenómeno (Gutiérrez, J).

2. PREPARACIÓN

Antes de empezar a diseñar nuestro cuaderno vamos a limpiar la memoria e instalar las bibliotecas necesarias.

rm(list=ls())

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

Ahora cargamos cada libreria.

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages ----------------------------------------------------------------------- tidyverse 1.3.0 --
v ggplot2 3.3.2     v purrr   0.3.4
v tibble  3.0.4     v dplyr   1.0.2
v tidyr   1.1.2     v stringr 1.4.0
v readr   1.4.0     v forcats 0.5.0
-- Conflicts -------------------------------------------------------------------------- tidyverse_conflicts() --
x tidyr::extract() masks raster::extract()
x dplyr::filter()  masks stats::filter()
x dplyr::lag()     masks stats::lag()
x dplyr::select()  masks raster::select()

library(readxl)

library(rgeos)

rgeos version: 0.5-5, (SVN revision 640)
 GEOS runtime version: 3.8.0-CAPI-1.13.1 
 Linking to sp version: 1.4-4 
 Polygon checking: TRUE 

library(raster)

library(sf)

Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(cartography)

library(SpatialPosition)

library(rgdal)

rgdal: version: 1.5-18, (SVN revision 1082)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
Path to GDAL shared files: C:/Users/laura/OneDrive/Documentos/R/win-library/4.0/rgdal/gdal
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
Path to PROJ shared files: C:/Users/laura/OneDrive/Documentos/R/win-library/4.0/rgdal/proj
Linking to sp version:1.4-4
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

3. DATOS

Utilizaremos datos sobre las Necesidades Básicas Insatisfechas (NBI) del Censo Nacional de Población y Vivienda 2018 descargados en formato .xlsx. (Descarguelos aquí).

El archivo descargado se modificó para lograr un mejor tratamiento de los datos, dejando solamente datos de municipios del departamento seleccionado, en nuestro caso Antioquía, adicionalmente se cambió el nombre de las variables por unos mas sencillos.

Finalmente el archivo modificado se carga en la carpeta de trabajo.

4. LECTURA DE LOS DATOS NBI

Ahora realizaremos la lectura de los datos anteriomente cargados correspondienteS a las estadísticas municipales agropecuarias para Antioquia.

nbi <- read_excel("NBI_Antioquia.xlsx")

Luego observemos los atributos

head(nbi)

Ahora podemos conocer el municipio con el porcentaje más alto de NBI

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

max_nbi

Se puede observar que el municipio con NBI mas alto es Murindó con un valor de 81,68 %. De igual manera averigüemos el municipio con el porcentaje más bajo de NBI

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

min_nbi

Se observa que el municipio con NBI mas bajo es Sabaneta con un valor de 1.58 %.

Ahora ordenemos los municipios en orden descendente.

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

desc_nbi

5. Unión de los datos del NBI a los municipios

Primero debemos subir los datos municipales de Antioquia.

munic <- st_read("ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\laura\Downloads\Geomatica\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
Simple feature collection with 125 features and 9 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -77.12783 ymin: 5.418558 xmax: -73.88128 ymax: 8.873974
geographic CRS: WGS 84

Vamos a comprobar lo que hay dentro del atributo MPIO_CCDGO

head(munic$MPIO_CNMBR)

[1] "MEDELLÍN"       "CIUDAD BOLÍVAR" "BRICEÑO"        "BURITICÁ"       "CÁCERES"        "BETULIA"       

Para la unión de los municipios y los datos de NBI usaremos la función left_join

nbi_munic <- left_join(munic, nbi, by=c("MPIO_CCDGO"="Codigo"))

nbi_munic %>%
  dplyr::select(MPIO_CNMBR, MPIO_CCDGO, NBI)  ->  check_nbi_munic

head(check_nbi_munic)

Simple feature collection with 6 features and 3 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -76.09953 ymin: 5.753381 xmax: -74.99361 ymax: 7.957635
geographic CRS: WGS 84
      MPIO_CNMBR MPIO_CCDGO       NBI                       geometry
1       MEDELLÍN      05001  5.192345 POLYGON ((-75.66974 6.37359...
2 CIUDAD BOLÍVAR      05101 11.708395 POLYGON ((-76.04467 5.92774...
3        BRICEÑO      05107 25.327952 POLYGON ((-75.45818 7.22284...
4       BURITICÁ      05113 28.908795 POLYGON ((-75.90857 6.97378...
5        CÁCERES      05120 49.846861 POLYGON ((-75.20358 7.95716...
6        BETULIA      05093 17.127878 POLYGON ((-76.00304 6.28171...

Ahora, vamos a reproyeccionar los municipios

nbi_munic_new <- st_transform(nbi_munic, crs = 3116)

6. Ejemplos de mapas temáticos

Utilizaremos el paquete de cartography que permite obtener mapas temáticos con buena calidad visual. Su biblioteca utiliza objetos sf o sp para producir los gráficos base.

Mapa base OpenStreetMap y símbolos proporcionales

Las funciones getTiles() y tilesLayer() permiten tener mosaicos OpenStreetMap. Por su parte propSymbolsLayer() muestra símbolos como círculos, cuadrados, barras con áreas proporcionales a una variable cuantitativa, además los tamaños pueden ser personalizados mediante el argumento inches.

mun.osm <- getTiles(
x = nbi_munic_new, 
type = "OpenStreetMap", 
zoom = 8,
cachedir = TRUE,
crop = FALSE
)

Discarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crsDiscarded datum World Geodetic System 1984 in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition

# set margins
opar <- par(mar = c(0,0,1.2,0))
# plot osm tiles
tilesLayer(x = mun.osm)
# plot municipalities (only borders are plotted)
plot(st_geometry(nbi_munic_new), col = NA, border = "chocolate3", add=TRUE)

# plot NBI
propSymbolsLayer(
  x = nbi_munic_new, 
  var = "NBI", 
  inches = 0.15, 
  col = "chocolate3",
  legend.pos = "topright",  
  legend.title.txt = "Total NBI"
)
# layout
layoutLayer(title = "Distribución de NBI en Antioquia",
            sources = "Sources: DANE, 2018",
            author = " Laura Forero ",
            frame = TRUE, north = FALSE, tabtitle = TRUE)

# north arrow
north(pos = "topleft")

Mapas Cloropléticos

Con este tipo de mapas se podrá observar áreas sombreadas según la variación de la variable cuantitativa. Mediante la función choroLayer() se mostrará el mapa de coropletas y argumentos como nclass, method y breaks permiten personalizar la clasificación de variables, mientras que getBreaks() permite clasificar fuera de la propia función. Las paletas de colores se definen con col y se puede variar según su preferencia.

# set margins
opar <- par(mar = c(0,0,1.2,0))
# set figure background color
par(bg="grey90")
# plot municipalities (only the backgroung color is plotted)
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 = "purple.pal", n1 = 5),
  border = "white", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "NBI",
  add = TRUE
) 

# layout
layoutLayer(title = "Distribución de NBI en Antioquia",
            sources = "Sources: DANE, 2018",
            author = " Laura Forero ",
            frame = TRUE, north = FALSE, tabtitle = TRUE, col="black") 
# north arrow
north(pos = "topleft")

Mapa de símbolos proporcionales y tipología

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

Primero, debemos crear una variable cualitativa usando la función mutate.

nbi_munic_2 <- dplyr::mutate(nbi_munic_new, pobreza = ifelse(Miseria > 20, "Extreme", 
                                                         ifelse(Hacinamiento > 5, "High", "Intermediate")))

Observe que el nuevo atributo se llama pobreza y su valor depende del umbral definido por el usuario

head(nbi_munic_2)

Simple feature collection with 6 features and 21 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 776023.9 ymin: 1128331 xmax: 898859.2 ymax: 1371903
projected CRS:  MAGNA-SIRGAS / Colombia Bogota zone
  DPTO_CCDGO MPIO_CCDGO     MPIO_CNMBR                                        MPIO_CRSLC MPIO_NAREA MPIO_NANO
1         05      05001       MEDELLÍN                                              1965   374.8280      2017
2         05      05101 CIUDAD BOLÍVAR                                              1869   260.4461      2017
3         05      05107        BRICEÑO              Ordenanza 27 de Noviembre 26 de 1980   376.3468      2017
4         05      05113       BURITICÁ                                              1812   355.2103      2017
5         05      05120        CÁCERES Decreto departamental 160 del 16 de Marzo de 1903  1873.8033      2017
6         05      05093        BETULIA Decreto departamental 629 del 28 de Enero de 1884   262.3675      2017
  DPTO_CNMBR Shape_Leng Shape_Area Cod_Depto Departamento Cod_Munic      Municipio       NBI    Miseria
1  ANTIOQUIA  1.0327835 0.03060723        05    ANTIOQUIA       001       MEDELLÍN  5.192345  0.4119415
2  ANTIOQUIA  0.7085039 0.02124224        05    ANTIOQUIA       101 CIUDAD BOLÍVAR 11.708395  1.2648358
3  ANTIOQUIA  1.0044720 0.03078496        05    ANTIOQUIA       107        BRICEÑO 25.327952  5.1463169
4  ANTIOQUIA  0.9637233 0.02902757        05    ANTIOQUIA       113       BURITICÁ 28.908795  6.5146580
5  ANTIOQUIA  2.9333643 0.15350440        05    ANTIOQUIA       120        CÁCERES 49.846861 18.4341501
6  ANTIOQUIA  0.8476756 0.02141352        05    ANTIOQUIA       093        BETULIA 17.127878  2.3491937
    Vivienda Servicios Hacinamiento Inasistencia  Economia                       geometry      poverty
1  0.2846426 0.1888295     1.536231     1.600220  2.039155 POLYGON ((823819.2 1196825,... Intermediate
2  0.4461578 0.5327904     2.395391     1.416443  8.334055 POLYGON ((782138 1147634, 7... Intermediate
3  3.6158762 3.3636058     7.484023     2.556340 14.278507 POLYGON ((847501.8 1290702,...         High
4  0.4614549 3.2844734    11.644951     6.460369 14.495114 POLYGON ((797631.6 1263319,...         High
5 38.0474732 8.9165391     7.392802     3.346095 14.873660 POLYGON ((875837.6 1371851,...         High
6  1.1281439 0.8096091     4.532484     2.455372 10.916451 POLYGON ((786889.9 1186783,... Intermediate

library(sf)
library(cartography)
# set margins
opar <- par(mar = c(0,0,1.2,0))
# Plot the municipalities
plot(st_geometry(nbi_munic_2), col="#f2efe9", border="#b38e43", bg = "#aad3df", 
     lwd = 0.5)
# Plot symbols with choropleth coloration
propSymbolsTypoLayer(
  x = nbi_munic_2, 
  var = "NBI", 
  inches = 0.3,
  symbols = "square",
  border = "white",
  lwd = .5,
  legend.var.pos = c(1050000, 1350000), 
  legend.var.title.txt = "NBI",
  var2 = "poverty",
  legend.var2.values.order = c("Extreme", "High", 
                               "Intermediate"),
  col = carto.pal(pal1 = "red.pal", n1 = 3),
  legend.var2.pos = c(1050000, 1200000), 
  legend.var2.title.txt = "Poverty"
) 

# layout
layoutLayer(title="Distribución de NBI en Antioquia", 
            author = "Laura Forero", 
            sources = "Source: DANE, 2018", 
            scale = 1, tabtitle = TRUE, frame = TRUE)
# north arrow
north(pos = "topleft")

Mapas de etiquetas

Se combinan las funciones choroLayer y labelLayer

library(sf)
library(cartography)
# set margins
opar <- par(mar = c(0,0,1.2,0))
# set figure background color
par(bg="white")
# plot municipalities
plot(st_geometry(nbi_munic_2), col = "#e4e9de", border = "darkseagreen4", 
     bg = "gray96", lwd = 0.5)
# plot NBI
choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "purple.pal", n1 = 5),
  border = "white", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "NBI",
  add = TRUE
) 

# plot labels
labelLayer(
  x = nbi_munic_2, 
  txt = "Municipio", 
  col= "white", 
  cex = 0.4, 
  font = 4,
  halo = TRUE, 
  bg = "grey25", 
  r = 0.06, 
  overlap = FALSE, 
  show.lines = FALSE
)
# map layout
layoutLayer(
  title = "Municipios de Antioquia", 
  sources = "Source: DANE, 2018",  
  author = "Laura Forero", 
  frame = TRUE,
  tabtitle = TRUE, 
  theme = "purple.pal"
) 

# north arrow
north(pos = "topleft")

Mapas de isopleta

Permite una representación espacial del fenómeno independiente de la heterogeneidad inicial de la división territorial. La función smoothLayer() utiliza una capa de puntos marcados y un conjunto de parámetros y muestra una capa de mapa de isóptica.

Para diseñar este mapa usaremos otro conjunto de datos correspondientes a estadísticas de producción de café para el 2018 en Antioquia (revisar este link)

Ahora leamos los datos

datos <- read_csv("Evaluaciones_Agropecuarias_Municipales_EVA.csv")

-- Column specification ----------------------------------------------------------------------------------------
cols(
  `C搼㸳D. 
DEP.` = col_double(),
  DEPARTAMENTO = col_character(),
  `C搼㸳D. MUN.` = col_double(),
  MUNICIPIO = col_character(),
  `GRUPO 
DE CULTIVO` = col_character(),
  `SUBGRUPO 
DE CULTIVO` = col_character(),
  CULTIVO = col_character(),
  `DESAGREGACI搼㸳N REGIONAL Y/O SISTEMA PRODUCTIVO` = col_character(),
  A搼㸱O = col_double(),
  PERIODO = col_character(),
  `挼㸱rea Sembrada
(ha)` = col_double(),
  `挼㸱rea Cosechada
(ha)` = col_double(),
  `Producci昼㸳n
(t)` = col_double(),
  `Rendimiento
(t/ha)` = col_double(),
  `ESTADO FISICO PRODUCCION` = col_character(),
  `NOMBRE 
CIENTIFICO` = col_character(),
  `CICLO DE CULTIVO` = col_character()
)

Ahora ejecutemos las siguientes líneas de código hasta obtener los datos de la producción de café para el año 2018

names(datos) <-c("Cod", "Departamento", "Cod_municipio", "Municipio", 
                 "Grupo", "Subgrupo", "Cultivo", "Sistema","Year", "Periodo", "Area_sembrada", "Area_cosechada",
                 "Produccion", "Rendimiento", "Estado", "Nombre",
                 "Ciclo")
datos

datos %>%
  filter(Departamento == "ANTIOQUIA") -> Antioquia
Antioquia

Antioquia %>% 
  filter(Year==2018) -> datos2
head(datos2)

datos2 %>% filter(Cultivo == "CAFE")  -> cafe2018

Ahora cree un nuevo atributo que coincida con los códigos de municipios

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

Realicemos la unión de los datos e inspeccionemos los mismos

cafe_munic = left_join(munic, cafe2018, by="MPIO_CCDGO")
head(cafe_munic)

Simple feature collection with 6 features and 27 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -76.09953 ymin: 5.753381 xmax: -74.99361 ymax: 7.957635
geographic CRS: WGS 84
  DPTO_CCDGO MPIO_CCDGO     MPIO_CNMBR                                        MPIO_CRSLC MPIO_NAREA MPIO_NANO
1         05      05001       MEDELLÍN                                              1965   374.8280      2017
2         05      05101 CIUDAD BOLÍVAR                                              1869   260.4461      2017
3         05      05107        BRICEÑO              Ordenanza 27 de Noviembre 26 de 1980   376.3468      2017
4         05      05113       BURITICÁ                                              1812   355.2103      2017
5         05      05120        CÁCERES Decreto departamental 160 del 16 de Marzo de 1903  1873.8033      2017
6         05      05093        BETULIA Decreto departamental 629 del 28 de Enero de 1884   262.3675      2017
  DPTO_CNMBR Shape_Leng Shape_Area Cod Departamento Cod_municipio      Municipio             Grupo Subgrupo
1  ANTIOQUIA  1.0327835 0.03060723   5    ANTIOQUIA          5001       MEDELLIN OTROS PERMANENTES     CAFE
2  ANTIOQUIA  0.7085039 0.02124224   5    ANTIOQUIA          5101 CIUDAD BOLIVAR OTROS PERMANENTES     CAFE
3  ANTIOQUIA  1.0044720 0.03078496   5    ANTIOQUIA          5107        BRICEÑO OTROS PERMANENTES     CAFE
4  ANTIOQUIA  0.9637233 0.02902757   5    ANTIOQUIA          5113       BURITICA OTROS PERMANENTES     CAFE
5  ANTIOQUIA  2.9333643 0.15350440  NA         <NA>            NA           <NA>              <NA>     <NA>
6  ANTIOQUIA  0.8476756 0.02141352   5    ANTIOQUIA          5093        BETULIA OTROS PERMANENTES     CAFE
  Cultivo Sistema Year Periodo Area_sembrada Area_cosechada Produccion Rendimiento                 Estado
1    CAFE    CAFE 2018    2018           471            433        299        0.69 CAFE VERDE EQUIVALENTE
2    CAFE    CAFE 2018    2018         10063           7374      11439        1.55 CAFE VERDE EQUIVALENTE
3    CAFE    CAFE 2018    2018           589            485        586        1.21 CAFE VERDE EQUIVALENTE
4    CAFE    CAFE 2018    2018           936            861        964        1.12 CAFE VERDE EQUIVALENTE
5    <NA>    <NA>   NA    <NA>            NA             NA         NA          NA                   <NA>
6    CAFE    CAFE 2018    2018          5734           4230       6562        1.55 CAFE VERDE EQUIVALENTE
          Nombre      Ciclo TEMP                       geometry
1 COFFEA ARABICA PERMANENTE 5001 POLYGON ((-75.66974 6.37359...
2 COFFEA ARABICA PERMANENTE 5101 POLYGON ((-76.04467 5.92774...
3 COFFEA ARABICA PERMANENTE 5107 POLYGON ((-75.45818 7.22284...
4 COFFEA ARABICA PERMANENTE 5113 POLYGON ((-75.90857 6.97378...
5           <NA>       <NA> <NA> POLYGON ((-75.20358 7.95716...
6 COFFEA ARABICA PERMANENTE 5093 POLYGON ((-76.00304 6.28171...

Ahora, reproyectemos los municipios:

rep_cafe <- st_transform(cafe_munic, crs = 3116)

Visualicemos el mapa

# set margins
opar <- par(mar = c(0,0,1.2,0))
# plot municipalities (only the backgroung color is plotted)
plot(st_geometry(rep_cafe), col = NA, border = "black", bg = "grey96")
# plot isopleth map
smoothLayer(
  x = rep_cafe, 
  var = 'Produccion',
  typefct = "exponential",
  span = 25000,
  beta = 2,
  nclass = 10,
  col = carto.pal(pal1 = 'red.pal', n1 = 10),
  border = "grey",
  lwd = 0.1, 
  mask = rep_cafe, 
  legend.values.rnd = -3,
  legend.title.txt = "Production",
  legend.pos = "topright", 
  add=TRUE
)

Discarded datum Unknown based on GRS80 ellipsoid in CRS definitionDiscarded datum Marco Geocentrico Nacional de Referencia in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definitionDiscarded datum Marco Geocentrico Nacional de Referencia in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definitionDiscarded datum Marco Geocentrico Nacional de Referencia in CRS definition

# annotation on the map
text(x = 650000, y = 1200000, cex = 0.6, adj = 0, font = 3,  labels = 
       "Distance function:\n- type = exponential\n- beta = 2\n- span = 20 km")
# layout
layoutLayer(title = "Distribución de producción de café en Antioquia",
            sources = "Sources: DANE and MADR, 2018",
            author = "Laura Forero",
            frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "red.pal")

# north arrow
north(pos = "topleft")

Guardar mapas

Vamos a crear otro mapa de la producción de café durante el año 2018, usando símbolos proporcionales y mapas de coropletas. El mapa final se guardará como un archivo .png.

El siguiente fragmento no visualiza un mapa, para ello escribe el mapa en el nombre de archivo cafe_2018.png en el directorio de trabajo.

### open the plot
png("cafe.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 = "white", 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 = "Coffe Production in 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 ANTIOQUIA, 2018",
            author = "Laura Forero", 
            sources = "Sources: MADR & DANE, 2018", 
            scale = 50, tabtitle = FALSE, frame = TRUE)
# north arrow
north(pos = "topleft")
#
title(main="Producción y rendimiento de café en Antioquia, 2018", cex.main=3,
      sub= "Source: MADR & DANE, 2018", cex.sub=2)
#
graticule = TRUE
#
par(opar)

knitr::include_graphics("cafe.png")