Cartografía temática

1. ¿Por qué este cuaderno?

Este es un cuaderno de R Markdown que ilustra cartografía temática para el departamento de Cauca en Colombia, Este cuaderno tiene como objetivo aprender a hacer mapas temáticos en R para la materia de Geomatica Básica de la Universidad Nacional.

2. Cartografía temática

Como dice Amanda Briney, un mapa temático enfatiza un tema o asunto, como la distribución promedio de las precipitaciones en un área o la densidad de población en un municipio. Se diferencian de los mapas de referencia generales porque no solo muestran características naturales y creadas por el hombre, como ríos, ciudades, subdivisiones políticas y carreteras. Si estos elementos aparecen en un mapa temático, son puntos de referencia para mejorar la comprensión del tema y el propósito del mapa.

Según Briney, “el factor más importante a considerar al diseñar mapas temáticos es la audiencia del mapa, que ayuda a determinar qué elementos deben incluirse en el mapa como puntos de referencia además del tema. Un mapa que se está haciendo para un científico político, por ejemplo, necesitaría mostrar límites políticos, mientras que uno para un biólogo podría necesitar contornos que muestren la elevación”.

Las fuentes de los datos de los mapas temáticos son muy importantes, enfatiza Briney. “Los cartógrafos deben encontrar fuentes de información precisas, recientes y confiables sobre una amplia gama de temas, desde características ambientales hasta datos demográficos, para hacer los mejores mapas posibles”.

Briney resume varias técnicas de mapeo temático que se utilizan con mayor frecuencia: ## 2.1 Mapa de coroplético que representa datos cuantitativos como un color y puede mostrar densidad, porcentaje, valor promedio o cantidad de un evento dentro de un área geográfica. Los colores secuenciales representan valores de datos positivos o negativos crecientes o decrecientes. Normalmente, cada color también representa un rango de valores.

2.2 Mapa de símbolos proporcionales o graduados

que se utilizan en otro tipo de mapa para representar datos asociados con ubicaciones, como ciudades. Los datos se muestran en estos mapas con símbolos de tamaño proporcional para mostrar diferencias en las ocurrencias. Los círculos se utilizan con mayor frecuencia, pero también son adecuados los cuadrados y otras formas geométricas. La forma más común de dimensionar estos símbolos es hacer que sus áreas sean proporcionales a los valores que se van a representar utilizando software de mapeo o dibujo.

2.3 Mapa isarítmico o de contorno

que utiliza isolíneas para representar valores continuos como los niveles de precipitación. Estos mapas también pueden mostrar valores tridimensionales, como la elevación, en mapas topográficos. Generalmente, los datos para mapas isarítmicos se recopilan a través de puntos medibles (por ejemplo, estaciones meteorológicas) o se recopilan por área (por ejemplo, toneladas de maíz por acre por condado). Los mapas isarítmicos también siguen la regla básica de que hay lados altos y bajos en relación con la isolina. Por ejemplo, en elevación, si la isolínea es de 500 pies, entonces un lado debe tener más de 500 pies y un lado debe ser más bajo.

2.4 Mapa de puntos

que utiliza puntos para mostrar la presencia de un tema y mostrar un patrón espacial. Un punto puede representar una unidad o varias, dependiendo de lo que se esté representando.

3. Datos

Usaremos datos sobre Necesidades Básicas Insatisfechas (NBI) del Censo Nacional de Población y Vivienda 2018 que están disponibles en el Geoportal del DANE.

Los mapas temáticos son útiles para transmitir información demográfica. Puede explorar varios mapas temáticos del DANE usando este enlace.

Anteriormente, descargué el archivo NBI, en formato .xlsx, a mi computadora. Luego, usé Excel para eliminar datos de municipios que no se encuentran dentro del Departamento de Cauca. También “limpié” los datos, que eliminé varias filas con imágenes institucionales o sin información. Estas filas se ubicaron tanto al principio como al final del archivo original. Conservo las columnas referentes a todo el municipio. Guardé los datos resultantes, correspondientes solo a los municipios de Cauca como NBI_Cauca.xlsx.

4. Preparación

Limpiemos la memoria:

rm(list=ls())

Ahora, instalemos las bibliotecas que necesitamos.

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)

Luego, carguemos las bibliotecas.

library(tidyverse)
library(readxl) ## this library is part of the tidyverse
library(rgeos)
library(raster)
library(sf)
library(cartography)
library(SpatialPosition)

5. Leer datos de NBI

Leamos el archivo csv con “estadisticas municipales agropecuarias” para Cauca.

nbi <- read_excel("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/informe dos/NBI_Cauca2.xlsx")

Veamos cuáles son los atributos de los datos.

head(nbi)

nbi2<-nbi
nbi2$cod_mun <-  as.character(nbi2$COD_MUN)
nbi2$MPIO_CCDGO <- as.factor(paste(19, nbi2$cod_mun, sep=""))

Busquemos cuál es el municipio con el porcentaje más alto de NBI:

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

max_nbi

Busquemos cuál es el municipio con el porcentaje más bajo de NBI:

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

min_nbi

Clasifiquemos los municipios por NBI en orden descendente:

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

desc_nbi

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

Comencemos a leer los datos administrativos de Cauca, dados por Marco Geoestadistico Departamental que se encuentra disponible en DANE Geoportal, usando la biblioteca sf:

munic <-st_read("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/Departamento/19_CAUCA/Cauca/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source `C:\Users\rojas\OneDrive\Escritorio\R estudio\Geomatica\Departamento\19_CAUCA\Cauca\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
Simple feature collection with 42 features and 9 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -77.92834 ymin: 0.9580285 xmax: -75.74782 ymax: 3.328941
geographic CRS: WGS 84

Veamos qué hay dentro del atributo MPIO_CCDGO:

head(munic$MPIO_CCDGO)

[1] "19001" "19022" "19050" "19075" "19100" "19110"

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

nbi_munic = left_join(munic, nbi2, by=c("MPIO_CCDGO"))

nbi_munic %>%
  dplyr::select(MUNICIPIO, 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: -77.40332 ymin: 1.635924 xmax: -76.38534 ymax: 3.233706
geographic CRS: WGS 84
     MUNICIPIO MPIO_CCDGO       NBI                       geometry
1      POPAYÁN      19001  8.798795 POLYGON ((-76.74145 2.57230...
2     ALMAGUER      19022 32.935908 POLYGON ((-76.80864 1.98674...
3      ARGELIA      19050 23.332497 POLYGON ((-77.3107 2.508326...
4       BALBOA      19075 20.128015 POLYGON ((-77.1431 2.156948...
5      BOLÍVAR      19100 27.177544 POLYGON ((-76.92984 2.11058...
6 BUENOS AIRES      19110 14.911238 POLYGON ((-76.76916 3.23128...

Ahora reproyectemos los datos

nbi_munic_new <- st_transform(nbi_munic, crs = 3116)

7. Ejemplos de mapas temáticos

Utilizaremos el paquete de cartografía que tiene como objetivo obtener mapas temáticos con la calidad visual de los construidos con un software clásico de cartografía o GIS. Sus funciones son intuitivas para los usuarios y garantizan la compatibilidad con los flujos de trabajo comunes de R.

La biblioteca de cartografía utiliza 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 base y símbolos proporcionales de OpenStreetMap

Las funciones getTiles () y tilesLayer () descargan y muestran mosaicos de OpenStreetMap. Tenga cuidado de citar la fuente de los mosaicos de manera adecuada.

propSymbolsLayer () muestra símbolos con áreas proporcionales a una variable cuantitativa (por ejemplo, NBI). Hay varios símbolos disponibles (círculos, cuadrados, barras). El argumento pulgadas se utiliza para personalizar los tamaños de los símbolos.

mun.osm <- getTiles(
x = nbi_munic_new, 
type = "OpenStreetMap", 
zoom = 10,
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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_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_defsDiscarded datum WGS_1984 in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

# 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 = "Black", add=TRUE)

# plot NBI
propSymbolsLayer(
  x = nbi_munic_new, 
  var = "NBI", 
  inches = 0.15, 
  col = "Red",
  legend.pos = "topright",  
  legend.title.txt = "Total NBI"
)
# layout
layoutLayer(title = " Distribución de NBI en Cauca",
            sources = "Sources: DANE, 2018\n© OpenStreetMap",
            author = " Johan Rojas",
            frame = TRUE, north = FALSE, tabtitle = TRUE)

# north arrow
north(pos = "topleft")

7. 2 mapas de coropleticos

En los mapas de coropleticos, las áreas están sombreadas de acuerdo con la variación de una variable cuantitativa. Se utilizan para representar proporciones o índices.

choroLayer () muestra mapas de coropleticos. Los argumentos nclass, método y rupturas permiten personalizar la clasificación de las variables.

getBreaks () permite clasificar fuera de la función en sí. Las paletas de colores se definen con col y se puede crear un conjunto de colores con carto.pal (). Consulte también display.carto.all ().

# establecer márgenes
opar <- par(mar = c(0,0,1.2,0))
# establecer color de fondo de las figuras 
par(bg="grey90")
# trazar municipios (solo se traza el color de fondo)
plot(st_geometry(nbi_munic_new), col = NA, border = NA, bg = "#aadaff")
# trazar NBI
choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "sand.pal", n1 = 5),
  border = "Black", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "NBI",
  add = TRUE
) 

# Diseño
layoutLayer(title = "Distribución de NBI en Cauca", 
            sources = "Source: DANE, 2018",
            author = "Johan Rojas", 
            frame = TRUE, north = TRUE, tabtitle = TRUE, col="black") 
# Flecha del norte 
north(pos = "topleft")

7.3 Mapa de símbolos proporcionales y 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 de calidad. Se utiliza una combinación de argumentos propSymbolsLayer () y typoLayer ().

Primero, necesitamos crear una variable cualitativa. Usemos la función mutate para esta tarea.

nbi_munic_2 <- dplyr::mutate(nbi_munic_new, pobreza = ifelse(MISERIA > 20, "Extrema", 
                                                         ifelse(HACINAMIENTO > 5, "Alta", "Intermedia")))

Tenga en cuenta que el nuevo atributo se llama pobreza. Su valor depende de los valores umbral definidos por el usuario. Asegúrese de verificar la sintaxis del comando ifelse para comprender lo que está sucediendo

head(nbi_munic_2)

Simple feature collection with 6 features and 22 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 629875.6 ymin: 672872.9 xmax: 743251.7 ymax: 849733.9
projected CRS:  MAGNA-SIRGAS / Colombia Bogota zone
  DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR                         MPIO_CRSLC
1         19      19001      POPAYÁN                               1537
2         19      19022     ALMAGUER                               1799
3         19      19050      ARGELIA Ordenanza 2 de Noviembre 8 de 1967
4         19      19075       BALBOA  Ordenanza 1 de Octubre 20 de 1967
5         19      19100      BOLÍVAR                               1793
6         19      19110 BUENOS AIRES                               1851
  MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area COD_DEPTO DEPTO
1   480.1798      2017      CAUCA   1.372937 0.03897000        19 CAUCA
2   236.7148      2017      CAUCA   0.710253 0.01919657        19 CAUCA
3   775.7878      2017      CAUCA   1.224153 0.06288063        19 CAUCA
4   413.2112      2017      CAUCA   1.016327 0.03349049        19 CAUCA
5   797.8203      2017      CAUCA   1.802787 0.06468337        19 CAUCA
6   435.0623      2017      CAUCA   1.633363 0.03532030        19 CAUCA
  COD_MUN    MUNICIPIO       NBI  MISERIA  VIVIENDA  SERVICIOS
1     001      POPAYÁN  8.798795 1.031876  4.481112  0.2561368
2     022     ALMAGUER 32.935908 5.555892  2.094051 24.0150094
3     050      ARGELIA 23.332497 3.723965 13.570891  1.1693852
4     075       BALBOA 20.128015 3.687050  5.062421  8.5749048
5     100      BOLÍVAR 27.177544 4.252632  2.377544 14.9726316
6     110 BUENOS AIRES 14.911238 2.092249  4.806625  3.3682042
  HACINAMIENTO INASISTENCIA ECONOMIA cod_mun
1     2.178812    0.9813082 2.067783     001
2     5.592205    1.6159293 6.233735     022
3     5.134253    3.2471769 4.577164     050
4     4.131401    1.5922556 4.914304     075
5     6.414035    2.3578947 5.959298     100
6     2.615311    1.8148676 4.632271     110
                        geometry    pobreza
1 POLYGON ((703640.6 776509.5... Intermedia
2 POLYGON ((696038 711703.1, ...       Alta
3 POLYGON ((640233.4 769570.3...       Alta
4 POLYGON ((658816 730615.4, ... Intermedia
5 POLYGON ((682562 725436.2, ...       Alta
6 POLYGON ((700729.5 849464.3... Intermedia

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="#CD853F", border="#000000", 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(845872.9, 759875.), 
  legend.var.title.txt = "NBI",
  var2 = "pobreza",
  legend.var2.values.order = c("Intermedia", "Alta", 
                               "Extrema"),
  col = carto.pal(pal1 = "orange.pal", n1 = 3),
  legend.var2.pos = c(845872.9 , 649875.6), 
  legend.var2.title.txt = "Pobreza"
) 

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

7.4 Etiquetar mapas

Intentemos combinar 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="grey25")
# plot municipalities
plot(st_geometry(nbi_munic_2), col = "#e4e9de", border = "Black", 
     bg = "grey75", lwd = 0.5)
# plot NBI
choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "wine.pal", n1 = 5),
  border = "Black", 
  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 = "Black", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
# map layout
layoutLayer(
  title = "Municipios de Cauca", 
  sources = "Source: DANE, 2018",  
  author = "Johan Rojas", 
  frame = TRUE,
  north = TRUE, 
  tabtitle = TRUE, 
  theme = "taupe.pal"
) 

7.5 Mapas de Isopleticos Los mapas de Isopleticos se basan en el supuesto de que el fenómeno a representar tiene una distribución continua. Estos mapas utilizan un enfoque de modelado de interacción espacial que tiene como objetivo calcular indicadores basados en valores de stock ponderados por distancia. Permite 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.

Usemos otro conjunto de datos para hacer un mapa de isopletas. En este caso, subiré datos estadísticos sobre la producción de café de 2018 en Antioquia. Ya conocemos este conjunto de datos, ya que se usó en un cuaderno anterior para ilustrar uniones basadas en atributos.

Leer el conjunto de datos:

cultivos2018 <-read_csv("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/trabajo7/EVA_cauca.csv")

Parsed with column specification:
cols(
  ID = col_double(),
  COD = col_character(),
  MUN = col_character(),
  GRUPO = col_character(),
  SUBGRUPO = col_character(),
  CULTIVO = col_character(),
  DESAGR = col_character(),
  YEAR = col_double(),
  PERIODO = col_character(),
  HA_SIEMBRA = col_double(),
  HA_COSECHA = col_double(),
  TON_PROD = col_double(),
  REND = col_double(),
  ESTADO = col_character(),
  NOMBRE = col_character(),
  CICLO = col_character()
)

head(cultivos2018) 

Filtrar filas que representan solo datos de café:

cultivos2018 %>%
  filter(CULTIVO == "CAFE") -> cafe2018
head(cafe2018)

Ahora, haz la unión

cafe_munic = left_join(munic, cafe2018, by=c("MPIO_CCDGO"="COD"))

Verifique la salida:

head(cafe_munic)

Simple feature collection with 6 features and 24 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -76.77328 ymin: 2.311382 xmax: -76.38534 ymax: 2.573465
geographic CRS: WGS 84
  DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CRSLC MPIO_NAREA MPIO_NANO
1         19      19001    POPAYÁN       1537   480.1798      2017
2         19      19001    POPAYÁN       1537   480.1798      2017
3         19      19001    POPAYÁN       1537   480.1798      2017
4         19      19001    POPAYÁN       1537   480.1798      2017
5         19      19001    POPAYÁN       1537   480.1798      2017
6         19      19001    POPAYÁN       1537   480.1798      2017
  DPTO_CNMBR Shape_Leng Shape_Area    ID     MUN             GRUPO
1      CAUCA   1.372937    0.03897 19000 POPAYAN OTROS PERMANENTES
2      CAUCA   1.372937    0.03897 19000 POPAYAN OTROS PERMANENTES
3      CAUCA   1.372937    0.03897 19000 POPAYAN OTROS PERMANENTES
4      CAUCA   1.372937    0.03897 19000 POPAYAN OTROS PERMANENTES
5      CAUCA   1.372937    0.03897 19000 POPAYAN OTROS PERMANENTES
6      CAUCA   1.372937    0.03897 19000 POPAYAN OTROS PERMANENTES
  SUBGRUPO CULTIVO DESAGR YEAR PERIODO HA_SIEMBRA HA_COSECHA TON_PROD REND
1     CAFE    CAFE   CAFE 2007    2007       3214       2182     2488 1.14
2     CAFE    CAFE   CAFE 2008    2008       3202       2678     2717 1.01
3     CAFE    CAFE   CAFE 2009    2009       3395       2390     2319 0.97
4     CAFE    CAFE   CAFE 2010    2010       3270       2294     2229 0.97
5     CAFE    CAFE   CAFE 2011    2011       3361       2358     2055 0.87
6     CAFE    CAFE   CAFE 2012    2012       3280       3240     2916 0.90
                  ESTADO         NOMBRE      CICLO
1 CAFE VERDE EQUIVALENTE COFFEA ARABICA PERMANENTE
2 CAFE VERDE EQUIVALENTE COFFEA ARABICA PERMANENTE
3 CAFE VERDE EQUIVALENTE COFFEA ARABICA PERMANENTE
4 CAFE VERDE EQUIVALENTE COFFEA ARABICA PERMANENTE
5 CAFE VERDE EQUIVALENTE COFFEA ARABICA PERMANENTE
6 CAFE VERDE EQUIVALENTE COFFEA ARABICA PERMANENTE
                        geometry
1 POLYGON ((-76.74145 2.57230...
2 POLYGON ((-76.74145 2.57230...
3 POLYGON ((-76.74145 2.57230...
4 POLYGON ((-76.74145 2.57230...
5 POLYGON ((-76.74145 2.57230...
6 POLYGON ((-76.74145 2.57230...

Ahora, reproyectemos los municipios:

rep_cafe <- st_transform(cafe_munic, crs = 3116)

ahora hacemos el mapeo

summary(rep_cafe$TON_PROD)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
      0     640    1276    2034    2683   11571       7 

# 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 = "grey75")
# plot isopleth map
smoothLayer(
  x = rep_cafe, 
  var = 'TON_PROD',
  typefct = "exponential",
  span = 12000,
  beta = 2,
  nclass = 10,
  col = carto.pal(pal1 = 'green.pal', n1 = 10),
  border = "Black",
  lwd = 0.1, 
  mask = rep_cafe, 
  legend.values.rnd = -3,
  legend.title.txt = "Producción",
  legend.pos = "topright", 
  add=TRUE
)

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

# annotation on the map
text(x = -76.57328, y = 2.411382, 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 la producción de cafe en Cauca",
            sources = "Sources: DANE and MADR, 2018",
            author = "Johan Rojas",
            frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "green.pal")

# north arrow
north(pos = "topleft")

8. Guardar mapas

Produzcamos otro mapa de la 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 su lugar, escribe el mapa con el nombre de archivo cafe_2018.png en el directorio de trabajo.

Por cuestiones de visualización vamos a dividir la producción de cafe en Cauca en dos partes, para esto vamos a dividir estos en dos grupos apartir de la media.

summary(cafe2018$TON_PROD)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0     640    1276    2034    2683   11571 

cafe2018 %>% 
  filter(YEAR==2018 & SUBGRUPO=="CAFE") %>% 
  group_by(MUN, CULTIVO,COD,TON_PROD,REND, HA_SIEMBRA) %>%
  summarize(sum_cosecha = sum(HA_COSECHA, na.rm = TRUE)) %>%
  filter(TON_PROD>=2034) %>%
  arrange(desc(sum_cosecha)) -> cafe20182

`summarise()` regrouping output by 'MUN', 'CULTIVO', 'COD', 'TON_PROD', 'REND' (override with `.groups` argument)

cafe20182

cafe2018 %>% 
  filter(YEAR==2018 & SUBGRUPO=="CAFE") %>% 
  group_by(MUN, CULTIVO,COD,TON_PROD,REND, HA_SIEMBRA) %>%
  summarize(sum_cosecha = sum(HA_COSECHA, na.rm = TRUE)) %>%
  filter(TON_PROD<2034) %>%
  arrange(desc(sum_cosecha)) -> cafe20182_1

`summarise()` regrouping output by 'MUN', 'CULTIVO', 'COD', 'TON_PROD', 'REND' (override with `.groups` argument)

cafe20182_1

cafe_munic2 = left_join(munic, cafe20182, by=c("MPIO_CCDGO"="COD"))

cafe_munic2_1 = left_join(munic, cafe20182_1, by=c("MPIO_CCDGO"="COD"))

rep_cafe2 <- st_transform(cafe_munic2, crs = 3116)

rep_cafe2_1 <- st_transform(cafe_munic2_1, crs = 3116)

### open the plot
png("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/informe dos/cafe_2018-2.png", width = 2048, height = 1526)
# set margins
opar <- par(mar = c(0,0,5,5))
# Plot the municipalities
plot(st_geometry(rep_cafe2), col="darkseagreen3", border="darkseagreen4",  
     bg = "white", lwd = 0.6)
# Plot symbols with choropleth coloration
propSymbolsChoroLayer(x = rep_cafe2, var = "TON_PROD", var2 = "REND",
                      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 = 3,
                      legend.values.cex = 1.0,
                      legend.var.pos = "right", 
                      legend.var2.pos = "left",
                      legend.var2.values.rnd = 2,
                      legend.var2.title.txt = "Rendimiento\n(en Ton/Ha)",
                      legend.var.title.txt = "producción de cafe (2018)",
                      legend.var.style = "e")
# plot labels
labelLayer(
  x = rep_cafe2, 
  txt = "MUN", 
  col= "Black", 
  cex = 1.0, 
  font = 5,
  halo = FALSE, 
  bg = "white", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
# layout
layoutLayer(title="Producción y rendimiento mayor a la media de café en Cauca, 2018",
            author = "Johan Rojas", 
            sources = "Sources: MADR & DANE, 2018", 
             tabtitle = FALSE, frame = TRUE,scale = "auto",)
# north arrow
north(pos = "topleft")
#
title(main="Producción y rendimiento mayor a la media de café en Cauca, 2018", cex.main=4,
      sub= "Source: MADR & DANE, 2018", cex.sub=4)
#
graticule = TRUE
#
par(opar)
### close the plot
dev.off()

null device 
          1 

png("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/informe dos/cafe_2018-2_1.png", width = 2048, height = 1526)
# set margins
opar <- par(mar = c(0,0,5,5))
# Plot the municipalities
plot(st_geometry(rep_cafe2_1), col="darkseagreen3", border="darkseagreen4",  
     bg = "white", lwd = 0.6)
# Plot symbols with choropleth coloration
propSymbolsChoroLayer(x = rep_cafe2_1, var = "TON_PROD", var2 = "REND",
                      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 = 3,
                      legend.values.cex = 1.0,
                      legend.var.pos = "right", 
                      legend.var2.pos = "left",
                      legend.var2.values.rnd = 2,
                      legend.var2.title.txt = "Rendimiento\n(en Ton/Ha)",
                      legend.var.title.txt = "producción de cafe (2018)",
                      legend.var.style = "e")
# plot labels
labelLayer(
  x = rep_cafe2_1, 
  txt = "MUN", 
  col= "Black", 
  cex = 1.0, 
  font = 5,
  halo = FALSE, 
  bg = "white", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
# layout
layoutLayer(title="Producción y rendimiento menor a la media de café en Cauca, 2018",
            author = "Johan Rojas", 
            sources = "Sources: MADR & DANE, 2018", 
             tabtitle = FALSE, frame = TRUE,scale = "auto",)
# north arrow
north(pos = "topleft")
#
title(main="Producción y rendimiento menor a la media de café en Cauca, 2018", cex.main=4,
      sub= "Source: MADR & DANE, 2018", cex.sub=5)
#
graticule = TRUE
#
par(opar)
### close the plot
dev.off()

null device 
          1 

knitr::include_graphics("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/informe dos/cafe_2018-2.png")

knitr::include_graphics("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/informe dos/cafe_2018-2_1.png")

Ahora vamos a ver la producción y área cosechada de café en Cauca en el 2018 de igual manera lo vamos ha hacer en dos mapas uno mayor a la media y otro menor a la media.

### open the plot
png("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/informe dos/cafe_2018c.png", width = 2048, height = 1526)
# set margins
opar <- par(mar = c(0,0,5,5))
# Plot the municipalities
plot(st_geometry(rep_cafe2), col="darkseagreen3", border="darkseagreen4",  
     bg = "white", lwd = 0.6)
# Plot symbols with choropleth coloration
propSymbolsChoroLayer(x = rep_cafe2, var = "TON_PROD", var2 = "sum_cosecha",
                      col = carto.pal(pal1 = "brown.pal", n1 = 3,
                                      pal2 = "orange.pal", n2 = 3),
                      inches = 0.8, method = "q6",
                      border = "grey50", lwd = 1,
                      legend.title.cex = 3,
                      legend.values.cex = 1.0,
                      legend.var.pos = "right", 
                      legend.var2.pos = "left",
                      legend.var2.values.rnd = 2,
                      legend.var2.title.txt = "Rendimiento\n(en Ton/Ha)",
                      legend.var.title.txt = "Área cosechada de cafe (2018)",
                      legend.var.style = "e")
# plot labels
labelLayer(
  x = rep_cafe2, 
  txt = "MUN", 
  col= "Black", 
  cex = 1.0, 
  font = 5,
  halo = FALSE, 
  bg = "white", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
# layout
layoutLayer(title="Producción y área cosechada mayor a la media de café en Cauca, 2018",
            author = "Johan Rojas", 
            sources = "Sources: MADR & DANE, 2018", 
             tabtitle = FALSE, frame = TRUE,scale = "auto",)
# north arrow
north(pos = "topleft")
#
title(main="Producción y área cosechada mayor a la media de café en Cauca, 2018", cex.main=4,
      sub= "Source: MADR & DANE, 2018", cex.sub=4)
#
graticule = TRUE
#
par(opar)
### close the plot
dev.off()

null device 
          1 

### open the plot
png("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/informe dos/cafe_2018c2.png", width = 2048, height = 1526)
# set margins
opar <- par(mar = c(0,0,5,5))
# Plot the municipalities
plot(st_geometry(rep_cafe2_1), col="darkseagreen3", border="darkseagreen4",  
     bg = "white", lwd = 0.6)
# Plot symbols with choropleth coloration
propSymbolsChoroLayer(x = rep_cafe2_1, var = "TON_PROD", var2 = "sum_cosecha",
                      col = carto.pal(pal1 = "brown.pal", n1 = 3,
                                      pal2 = "orange.pal", n2 = 3),
                      inches = 0.8, method = "q6",
                      border = "grey50", lwd = 1,
                      legend.title.cex = 3,
                      legend.values.cex = 1.0,
                      legend.var.pos = "right", 
                      legend.var2.pos = "left",
                      legend.var2.values.rnd = 2,
                      legend.var2.title.txt = "Rendimiento\n(en Ton/Ha)",
                      legend.var.title.txt = "Área cosechada de cafe (2018)",
                      legend.var.style = "e")
# plot labels
labelLayer(
  x = rep_cafe2_1, 
  txt = "MUN", 
  col= "Black", 
  cex = 1.0, 
  font = 5,
  halo = FALSE, 
  bg = "white", 
  r = 0.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
# layout
layoutLayer(title="Producción y área cosechada menor a la media de café en Cauca, 2018",
            author = "Johan Rojas", 
            sources = "Sources: MADR & DANE, 2018", 
             tabtitle = FALSE, frame = TRUE,scale = "auto",)
# north arrow
north(pos = "topleft")
#
title(main="Producción y área cosechada menor a la media de café en Cauca, 2018", cex.main=4,
      sub= "Source: MADR & DANE, 2018", cex.sub=4)
#
graticule = TRUE
#
par(opar)
### close the plot
dev.off()

null device 
          1 

knitr::include_graphics("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/informe dos/cafe_2018c.png")

knitr::include_graphics("C:/Users/rojas/OneDrive/Escritorio/R estudio/Geomatica/informe dos/cafe_2018c2.png")

Ahora vamos a ver hacer un mapa cloropletico con el área de siembra en Cauca para 2018

cafe2018 %>% 
  filter(YEAR==2018 & SUBGRUPO=="CAFE") %>% 
  group_by(MUN, CULTIVO,COD,TON_PROD,REND, HA_SIEMBRA) %>%
  summarize(sum_siembra = sum(HA_SIEMBRA, na.rm = TRUE)) -> cafe2018si

`summarise()` regrouping output by 'MUN', 'CULTIVO', 'COD', 'TON_PROD', 'REND' (override with `.groups` argument)

cafe2018si

cafe_municsi = left_join(munic, cafe2018si, by=c("MPIO_CCDGO"="COD"))

rep_cafesi <- st_transform(cafe_munic2, crs = 3116)

library(sf)
library(cartography)
# set margins
opar <- par(mar = c(0,0,1.2,0))
# set figure background color
par(bg="grey25")
# plot municipalities
plot(st_geometry(rep_cafesi), col = "#e4e9de", border = "Black", 
     bg = "grey75", lwd = 0.5)
# plot NBI
choroLayer(
  x = rep_cafesi, 
  var = "HA_SIEMBRA",
  method = "geom",
  nclass=7,
  col = carto.pal(pal1 = "wine.pal", n1 = 7),
  border = "Black", 
  lwd = 0.5,
  legend.pos = "topright",
  legend.title.txt = "ÁREA DE SIEMBRA",
  add = TRUE
) 

# plot labels
labelLayer(
  x =rep_cafesi, 
  txt = "MUN", 
  col= "Black", 
  cex = 0.4, 
  font = 4,
  overlap = FALSE, 
  show.lines = FALSE
)
# map layout
layoutLayer(
  title = "Área de siembra en los municipios de Cauca", 
  sources = "Source: DANE, 2018",  
  author = "Johan Rojas", 
  frame = TRUE,
  north = TRUE, 
  tabtitle = TRUE, 
  theme = "taupe.pal")

ahora miremos la distribución del área de siembra con un mapa isopletico.

# set margins
opar <- par(mar = c(0,0,1.2,0))
# plot municipalities (only the backgroung color is plotted)
plot(st_geometry(rep_cafesi), col = NA, border = "Black", bg = "grey75")
# plot isopleth map
smoothLayer(
  x = rep_cafesi, 
  var = 'HA_SIEMBRA',
  typefct = "exponential",
  span = 8428,
  beta = 2,
  nclass = 10,
  col = carto.pal(pal1 = 'green.pal', n1 = 10),
  border = "Black",
  lwd = 0.1, 
  mask = rep_cafe, 
  legend.values.rnd = -3,
  legend.title.txt = "Área de siembra",
  legend.pos = "topright", 
  add=TRUE
)

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

# annotation on the map
text(x = -76.57328, y = 2.411382, 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 del área de siembra de cafe en Cauca",
            sources = "Sources: DANE and MADR, 2018",
            author = "Johan Rojas",
            frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "green.pal")

# north arrow
north(pos = "topleft")

Por ultimo vamos a ver como se relaciona el nbi y la producción en Cauca para 2018, primero unamos los datos que necesitamos.

cafe_munic_nbi = left_join(nbi_munic, cafe2018si, by=c("MPIO_CCDGO"="COD"))

cafe_munic_nbi

Simple feature collection with 42 features and 27 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -77.92834 ymin: 0.9580285 xmax: -75.74782 ymax: 3.328941
geographic CRS: WGS 84
First 10 features:
   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR                         MPIO_CRSLC
1          19      19001      POPAYÁN                               1537
2          19      19022     ALMAGUER                               1799
3          19      19050      ARGELIA Ordenanza 2 de Noviembre 8 de 1967
4          19      19075       BALBOA  Ordenanza 1 de Octubre 20 de 1967
5          19      19100      BOLÍVAR                               1793
6          19      19110 BUENOS AIRES                               1851
7          19      19130      CAJIBÍO                               1824
8          19      19137      CALDONO                               1746
9          19      19142       CALOTO                               1543
10         19      19212      CORINTO                               1868
   MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area COD_DEPTO DEPTO
1    480.1798      2017      CAUCA  1.3729371 0.03897000        19 CAUCA
2    236.7148      2017      CAUCA  0.7102530 0.01919657        19 CAUCA
3    775.7878      2017      CAUCA  1.2241529 0.06288063        19 CAUCA
4    413.2112      2017      CAUCA  1.0163274 0.03349049        19 CAUCA
5    797.8203      2017      CAUCA  1.8027866 0.06468337        19 CAUCA
6    435.0623      2017      CAUCA  1.6333631 0.03532030        19 CAUCA
7    552.2877      2017      CAUCA  1.3822054 0.04482019        19 CAUCA
8    354.7387      2017      CAUCA  1.1029614 0.02880277        19 CAUCA
9    266.5987      2017      CAUCA  1.3621191 0.02164894        19 CAUCA
10   325.9156      2017      CAUCA  0.8276564 0.02647947        19 CAUCA
   COD_MUN    MUNICIPIO       NBI  MISERIA  VIVIENDA  SERVICIOS
1      001      POPAYÁN  8.798795 1.031876  4.481112  0.2561368
2      022     ALMAGUER 32.935908 5.555892  2.094051 24.0150094
3      050      ARGELIA 23.332497 3.723965 13.570891  1.1693852
4      075       BALBOA 20.128015 3.687050  5.062421  8.5749048
5      100      BOLÍVAR 27.177544 4.252632  2.377544 14.9726316
6      110 BUENOS AIRES 14.911238 2.092249  4.806625  3.3682042
7      130      CAJIBÍO 22.860245 4.004964  7.035831  5.7701256
8      137      CALDONO 22.278883 3.895670  4.062866  3.9346820
9      142       CALOTO 22.609246 4.678594 10.995883  7.6987017
10     212      CORINTO 18.679795 3.083252  9.721259  3.0601396
   HACINAMIENTO INASISTENCIA ECONOMIA cod_mun          MUN CULTIVO
1      2.178812    0.9813082 2.067783     001      POPAYAN    CAFE
2      5.592205    1.6159293 6.233735     022     ALMAGUER    CAFE
3      5.134253    3.2471769 4.577164     050      ARGELIA    CAFE
4      4.131401    1.5922556 4.914304     075       BALBOA    CAFE
5      6.414035    2.3578947 5.959298     100      BOLIVAR    CAFE
6      2.615311    1.8148676 4.632271     110 BUENOS AIRES    CAFE
7      3.762991    2.5779432 8.279820     130      CAJIBIO    CAFE
8      7.939029    1.7360531 8.816809     137      CALDONO    CAFE
9      4.045282    0.8391387 4.694427     142       CALOTO    CAFE
10     3.804373    1.1695095 4.701151     212      CORINTO    CAFE
   TON_PROD REND HA_SIEMBRA sum_siembra                       geometry
1      2939 0.95       3510        3510 POLYGON ((-76.74145 2.57230...
2       819 0.78       1157        1157 POLYGON ((-76.80864 1.98674...
3       830 0.62       1494        1494 POLYGON ((-77.3107 2.508326...
4      2226 0.73       3352        3352 POLYGON ((-77.1431 2.156948...
5      4598 1.87       2682        2682 POLYGON ((-76.92984 2.11058...
6      1446 0.86       1825        1825 POLYGON ((-76.76916 3.23128...
7      9500 1.55       6860        6860 POLYGON ((-76.83997 2.76373...
8      5476 1.29       4936        4936 POLYGON ((-76.36326 2.90734...
9      2265 1.29       1973        1973 POLYGON ((-76.42866 3.20990...
10     1280 0.90       1660        1660 POLYGON ((-76.28278 3.22879...

cafe_munic_nbi2 = st_join(munic, cafe_munic_nbi, by=c("MPIO_CCDGO"))

although coordinates are longitude/latitude, st_intersects assumes that they are planar

Ahora vamos a crear una escala (alta, media, baja) para la producción de cafe en 2018, para esto usaremos los datos de los cuartiles de esta variable.

summary(cafe_munic_nbi2$TON_PROD)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
    228     830    2570    3529    4635   11315      48 

cafe_nbi_munic_2 <- dplyr::mutate(cafe_munic_nbi, PROD = ifelse(TON_PROD>4635, "Alta", 
                                                         ifelse(TON_PROD > 3529, "Media", "Baja")))

cafe_nbi_munic_2

Simple feature collection with 42 features and 28 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -77.92834 ymin: 0.9580285 xmax: -75.74782 ymax: 3.328941
geographic CRS: WGS 84
First 10 features:
   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR                         MPIO_CRSLC
1          19      19001      POPAYÁN                               1537
2          19      19022     ALMAGUER                               1799
3          19      19050      ARGELIA Ordenanza 2 de Noviembre 8 de 1967
4          19      19075       BALBOA  Ordenanza 1 de Octubre 20 de 1967
5          19      19100      BOLÍVAR                               1793
6          19      19110 BUENOS AIRES                               1851
7          19      19130      CAJIBÍO                               1824
8          19      19137      CALDONO                               1746
9          19      19142       CALOTO                               1543
10         19      19212      CORINTO                               1868
   MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area COD_DEPTO DEPTO
1    480.1798      2017      CAUCA  1.3729371 0.03897000        19 CAUCA
2    236.7148      2017      CAUCA  0.7102530 0.01919657        19 CAUCA
3    775.7878      2017      CAUCA  1.2241529 0.06288063        19 CAUCA
4    413.2112      2017      CAUCA  1.0163274 0.03349049        19 CAUCA
5    797.8203      2017      CAUCA  1.8027866 0.06468337        19 CAUCA
6    435.0623      2017      CAUCA  1.6333631 0.03532030        19 CAUCA
7    552.2877      2017      CAUCA  1.3822054 0.04482019        19 CAUCA
8    354.7387      2017      CAUCA  1.1029614 0.02880277        19 CAUCA
9    266.5987      2017      CAUCA  1.3621191 0.02164894        19 CAUCA
10   325.9156      2017      CAUCA  0.8276564 0.02647947        19 CAUCA
   COD_MUN    MUNICIPIO       NBI  MISERIA  VIVIENDA  SERVICIOS
1      001      POPAYÁN  8.798795 1.031876  4.481112  0.2561368
2      022     ALMAGUER 32.935908 5.555892  2.094051 24.0150094
3      050      ARGELIA 23.332497 3.723965 13.570891  1.1693852
4      075       BALBOA 20.128015 3.687050  5.062421  8.5749048
5      100      BOLÍVAR 27.177544 4.252632  2.377544 14.9726316
6      110 BUENOS AIRES 14.911238 2.092249  4.806625  3.3682042
7      130      CAJIBÍO 22.860245 4.004964  7.035831  5.7701256
8      137      CALDONO 22.278883 3.895670  4.062866  3.9346820
9      142       CALOTO 22.609246 4.678594 10.995883  7.6987017
10     212      CORINTO 18.679795 3.083252  9.721259  3.0601396
   HACINAMIENTO INASISTENCIA ECONOMIA cod_mun          MUN CULTIVO
1      2.178812    0.9813082 2.067783     001      POPAYAN    CAFE
2      5.592205    1.6159293 6.233735     022     ALMAGUER    CAFE
3      5.134253    3.2471769 4.577164     050      ARGELIA    CAFE
4      4.131401    1.5922556 4.914304     075       BALBOA    CAFE
5      6.414035    2.3578947 5.959298     100      BOLIVAR    CAFE
6      2.615311    1.8148676 4.632271     110 BUENOS AIRES    CAFE
7      3.762991    2.5779432 8.279820     130      CAJIBIO    CAFE
8      7.939029    1.7360531 8.816809     137      CALDONO    CAFE
9      4.045282    0.8391387 4.694427     142       CALOTO    CAFE
10     3.804373    1.1695095 4.701151     212      CORINTO    CAFE
   TON_PROD REND HA_SIEMBRA sum_siembra                       geometry
1      2939 0.95       3510        3510 POLYGON ((-76.74145 2.57230...
2       819 0.78       1157        1157 POLYGON ((-76.80864 1.98674...
3       830 0.62       1494        1494 POLYGON ((-77.3107 2.508326...
4      2226 0.73       3352        3352 POLYGON ((-77.1431 2.156948...
5      4598 1.87       2682        2682 POLYGON ((-76.92984 2.11058...
6      1446 0.86       1825        1825 POLYGON ((-76.76916 3.23128...
7      9500 1.55       6860        6860 POLYGON ((-76.83997 2.76373...
8      5476 1.29       4936        4936 POLYGON ((-76.36326 2.90734...
9      2265 1.29       1973        1973 POLYGON ((-76.42866 3.20990...
10     1280 0.90       1660        1660 POLYGON ((-76.28278 3.22879...
    PROD
1   Baja
2   Baja
3   Baja
4   Baja
5  Media
6   Baja
7   Alta
8   Alta
9   Baja
10  Baja

pobrecafe <- st_transform(cafe_nbi_munic_2, crs = 3116)

Ahora hagamos al mapa

library(sf)
library(cartography)
# set margins
opar <- par(mar = c(0,0,1.2,0))
# Plot the municipalities
plot(st_geometry(pobrecafe), col="#CD853F", border="#000000", bg = "#aad3df", 
     lwd = 0.5)
# Plot symbols with choropleth coloration
propSymbolsTypoLayer(
  x = pobrecafe, 
  var = "NBI", 
  inches = 0.3,
  symbols = "square",
  border = "Black",
  lwd = .5,
  legend.var.pos = c(845872.9, 759875.), 
  legend.var.title.txt = "NBI",
  var2 = "PROD",
  legend.var2.values.order = c("Alta", "Media", 
                               "Baja"),
  col = carto.pal(pal1 = "orange.pal", n1 = 3),
  legend.var2.pos = c(845872.9 , 649875.6), 
  legend.var2.title.txt = "Producción"
) 

labelLayer(
  x =pobrecafe, 
  txt = "MPIO_CNMBR", 
  col= "Black", 
  cex = 0.4, 
  font = 4,
  overlap = FALSE, 
  show.lines = FALSE)
# layout
layoutLayer(title="Distribución de NBI y producción de cafe en Cauca", 
            author = "Johan Rojas", 
            sources = "Source: DANE & MADR, 2018", 
            scale = 1, tabtitle = TRUE, frame = TRUE)

# north arrow
north(pos = "topleft")

sessionInfo()

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Mexico.1252  LC_CTYPE=Spanish_Mexico.1252   
[3] LC_MONETARY=Spanish_Mexico.1252 LC_NUMERIC=C                   
[5] LC_TIME=Spanish_Mexico.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] SpatialPosition_2.0.1 cartography_2.4.2     sf_0.9-5             
 [4] rgeos_0.5-5           readxl_1.3.1          forcats_0.5.0        
 [7] stringr_1.4.0         dplyr_1.0.2           purrr_0.3.4          
[10] readr_1.3.1           tidyr_1.1.2           tibble_3.0.3         
[13] ggplot2_3.3.2         tidyverse_1.3.0       raster_3.3-13        
[16] sp_1.4-2             

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5         lubridate_1.7.9    lattice_0.20-41   
 [4] png_0.1-7          class_7.3-17       assertthat_0.2.1  
 [7] digest_0.6.25      R6_2.4.1           cellranger_1.1.0  
[10] backports_1.1.7    reprex_0.3.0       evaluate_0.14     
[13] e1071_1.7-3        httr_1.4.2         pillar_1.4.6      
[16] rlang_0.4.7        rstudioapi_0.11    blob_1.2.1        
[19] rmarkdown_2.3      rgdal_1.5-16       munsell_0.5.0     
[22] broom_0.7.0        compiler_4.0.2     modelr_0.1.8      
[25] xfun_0.16          base64enc_0.1-3    pkgconfig_2.0.3   
[28] htmltools_0.5.0    tidyselect_1.1.0   codetools_0.2-16  
[31] fansi_0.4.1        crayon_1.3.4       dbplyr_1.4.4      
[34] withr_2.2.0        grid_4.0.2         lwgeom_0.2-5      
[37] jsonlite_1.7.0     gtable_0.3.0       lifecycle_0.2.0   
[40] DBI_1.1.0          magrittr_1.5       units_0.6-7       
[43] scales_1.1.1       KernSmooth_2.23-17 cli_2.0.2         
[46] stringi_1.4.6      fs_1.5.0           xml2_1.3.2        
[49] ellipsis_0.3.1     generics_0.0.2     vctrs_0.3.3       
[52] tools_4.0.2        glue_1.4.2         hms_0.5.3         
[55] slippymath_0.3.1   yaml_2.2.1         colorspace_1.4-1  
[58] classInt_0.4-3     rvest_0.3.6        isoband_0.2.2     
[61] knitr_1.29         haven_2.3.1