Informe 2. Cartografia tematica

1.Thematic Cartography for Caqueta Department

Este es un cuaderno de notas de R Markdown que ilustra la cartografía temática del departamento de caqueta en Colombia. Ha sido creado, compilado y publicado a partir de RStudio Cloud.

2.1. Cartografía temática

Un mapa temático hace hincapié en un tema o asunto, como la distribución media de las precipitaciones en un área o la densidad de población en un municipio. Se diferencian de los mapas de referencia general porque no sólo muestran características naturales y artificiales 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.

Briney resume varias técnicas de cartografía temática que se utilizan con mayor frecuencia:

2.2. Mapa de Choropleth

que retrata los datos cuantitativos como un color y puede mostrar la densidad, el porcentaje, el valor medio o la 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.3. Mapa de símbolos proporcionales o graduados

que se utilizan en otro tipo de mapa para representar datos asociados a lugares, como las ciudades. Los datos se muestran en estos mapas con símbolos de tamaño proporcional para mostrar las diferencias en las ocurrencias. Los círculos son los más utilizados, pero los cuadrados y otras formas geométricas también son adecuados. 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 usando software de mapeo o dibujo.

2.4. 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 los mapas topográficos. Por lo general, los datos de los mapas isarítmicos se reúnen por medio de puntos medibles (por ejemplo, las estaciones meteorológicas) o se recogen por zonas (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 isolínea. Por ejemplo, en elevación, si la isolínea está a 500 pies, entonces un lado debe ser más alto que 500 pies y un lado debe ser más bajo.

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

Utilizaremos los 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.

4. Limpiemos la memoria

rm(list=ls())

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)

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.3     v dplyr   1.0.2
v tidyr   1.1.2     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(readxl) 
library(rgeos)

Loading required package: sp
rgeos version: 0.5-5, (SVN revision 640)
 GEOS runtime version: 3.8.0-CAPI-1.13.1 
 Linking to sp version: 1.4-2 
 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.8.0, GDAL 3.0.4, PROJ 6.3.1

library(cartography)
library(SpatialPosition)

5. Ahora cargamos los datos, que fueron previamente filtrados y copiados en un nuevo archivo para el departamento de caqueta a partir de los datos de Necesidades Básicas Insatisfechas (NBI) del Censo Nacional de Población y Vivienda 2018 que pueden ser descargados aquí.

nbi <- read_excel("d:/Users/Janus/Documents/Geomatica basica/18_CAQUETA/Caqueta_NBI.xlsx")

nbi

Encontremos que el municipio de Milan es el mayor porcentaje de necesidades basicas insatisfechas (NBI):

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

max_nbi

Encontremos cque Florencia es el municipio con el porcentaje más bajo de necesidades basicas insatisfechas NBI:

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

min_nbi

Clasificacion de municipios por NBI en orden descendente:

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

desc_nbi

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

Vamos a leer los datos usando la biblioteca sf:

Caq_MUN <-sf::st_read("d:/Users/Janus/Documents/Basic Geomatics/18_CAQUETA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")

Reading layer `MGN_MPIO_POLITICO' from data source `D:\Users\Janus\Documents\Basic Geomatics\18_CAQUETA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
Simple feature collection with 16 features and 9 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -76.30622 ymin: -0.70584 xmax: -71.25385 ymax: 2.964148
geographic CRS: WGS 84

Comprobemos lo que hay dentro del atributo MPIO_CCDGO:

(Caq_MUN$MPIO_CNMBR)

 [1] "FLORENCIA"              "ALBANIA"               
 [3] "BELÉN DE LOS ANDAQUÍES" "EL DONCELLO"           
 [5] "EL PAUJIL"              "LA MONTAÑITA"          
 [7] "MILÁN"                  "MORELIA"               
 [9] "SAN JOSÉ DEL FRAGUA"    "VALPARAÍSO"            
[11] "CARTAGENA DEL CHAIRÁ"   "SAN VICENTE DEL CAGUÁN"
[13] "SOLANO"                 "SOLITA"                
[15] "CURILLO"                "PUERTO RICO"           

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

nbi_munic = left_join(Caq_MUN, nbi, by=c("MPIO_CCDGO"="CODIGO"))

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: -76.1027 ymin: 0.9764735 xmax: -74.89527 ymax: 2.326755
geographic CRS: WGS 84
               Municipio MPIO_CCDGO      NBI
1              FLORENCIA      18001 14.44902
2                ALBANIA      18029 19.80289
3 BELÉN DE LOS ANDAQUIES      18094 28.28503
4            EL DONCELLO      18247 20.37901
5              EL PAUJIL      18256 23.33489
6           LA MONTAÑITA      18410 32.51475
                        geometry
1 POLYGON ((-75.42074 2.19413...
2 POLYGON ((-75.89506 1.36569...
3 POLYGON ((-75.78705 1.74982...
4 POLYGON ((-75.36167 2.32142...
5 POLYGON ((-75.36691 2.21234...
6 POLYGON ((-75.40404 1.76944...

check_nbi_munic

Simple feature collection with 16 features and 3 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -76.30622 ymin: -0.70584 xmax: -71.25385 ymax: 2.964148
geographic CRS: WGS 84
First 10 features:
                Municipio MPIO_CCDGO      NBI
1               FLORENCIA      18001 14.44902
2                 ALBANIA      18029 19.80289
3  BELÉN DE LOS ANDAQUIES      18094 28.28503
4             EL DONCELLO      18247 20.37901
5               EL PAUJIL      18256 23.33489
6            LA MONTAÑITA      18410 32.51475
7                   MILÁN      18460 46.96275
8                 MORELIA      18479 20.17910
9     SAN JOSÉ DEL FRAGUA      18610 26.27451
10             VALPARAÍSO      18860 27.38640
                         geometry
1  POLYGON ((-75.42074 2.19413...
2  POLYGON ((-75.89506 1.36569...
3  POLYGON ((-75.78705 1.74982...
4  POLYGON ((-75.36167 2.32142...
5  POLYGON ((-75.36691 2.21234...
6  POLYGON ((-75.40404 1.76944...
7  POLYGON ((-75.39362 1.35738...
8  POLYGON ((-75.77185 1.57991...
9  POLYGON ((-76.16722 1.58752...
10 POLYGON ((-75.73128 1.32740...

Reproyectemos los municipios:

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 que se construyen con un software de cartografía clásica o SIG.

Como la mayoría de los elementos internos del paquete se basan en las funciones sf, el formato preferido para los objetos espaciales es sf.

7.1. Símbolos de base y proporcionales de OpenStreetMap

Las funciones getTiles() y tilesLayer() descargan y muestran los azulejos de OpenStreetMap. Tengan cuidado de citar la fuente de los azulejos apropiadamente.

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 de las pulgadas se utiliza para personalizar el tamaño de los símbolos.

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

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

# plot NBI
propSymbolsLayer(
  x = nbi_munic_new, 
  var = "NBI", 
  inches = 0.15, 
  col = "brown4",
  legend.pos = "topright",  
  legend.title.txt = "Total NBI"
)
# layout
layoutLayer(title = " NBI Distribution in Caqueta",
            sources = "Sources: DANE, 2018\n© OpenStreetMap",
            author = " Luis Castillo ",
            frame = TRUE, north = FALSE, tabtitle = TRUE)

# north arrow
north(pos = "topleft")

7.2 mapas de Choropleth

Ahora vamos a usar el fondo de OpenStreetMap y trazar el NBI con diferentes intensidades de color en su lugar, lo que también se conoce como mapeo coropleto.

Usaremos la función choroLayer(), mientras que la función getBreaks() permite clasificar fuera de la función en sí. Además, las paletas de colores se definen con col y se puede crear un conjunto de colores con carto.pal().

opar <- par(mar = c(0,0,1.2,0))

par(bg="grey90")

plot(st_geometry(nbi_munic_new), col = NA, border = NA, bg = "#aadaff")

choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "sand.pal", n1 = 5),
  border = "white", 
  lwd = 0.5,
  legend.pos = "topright", 
  legend.title.txt = "NBI",
  add = TRUE
) 

layoutLayer(title = "NBI Distribution in caqueta", 
            sources = "Source: DANE, 2018",
            author = "Luis castillo", 
            frame = TRUE, north = TRUE, tabtitle = TRUE, col="black") 

north(pos = "topleft")

7.3. Símbolos proporcionales y mapa tipológico

Vamos a crear una nueva variable creando una nueva categoría llamada Pobreza usando las categorías de Miseria y Hacinamiento.

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

Vamos a hacer un mapa usando las categorías de NBI y Pobreza juntas.

head(nbi_munic_2)

nbi_munic_2

Simple feature collection with 16 features and 21 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 751894.5 ymin: 413692 xmax: 1314457 ymax: 819562.9
projected CRS:  MAGNA-SIRGAS / Colombia Bogota zone
First 10 features:
   DPTO_CCDGO MPIO_CCDGO             MPIO_CNMBR
1          18      18001              FLORENCIA
2          18      18029                ALBANIA
3          18      18094 BELÉN DE LOS ANDAQUÍES
4          18      18247            EL DONCELLO
5          18      18256              EL PAUJIL
6          18      18410           LA MONTAÑITA
7          18      18460                  MILÁN
8          18      18479                MORELIA
9          18      18610    SAN JOSÉ DEL FRAGUA
10         18      18860             VALPARAÍSO
                              MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR
1        Decreto 642 de Junio 17 de 1912  2547.6384      2017    CAQUETÁ
2    Ordenanza 3 de Noviembre 12 de 1985   414.1220      2017    CAQUETÁ
3        Decreto 963 de Marzo 14 de 1950  1191.6187      2017    CAQUETÁ
4   Decreto 1678 de Septiembre 7 de 1967  1105.8029      2017    CAQUETÁ
5  Decreto 1678 de Septiembret 7 de 1967  1234.7427      2017    CAQUETÁ
6          Decreto 83 de Julio 6 de 1955  1701.0522      2017    CAQUETÁ
7    Ordenanza 3 de Noviembre 12 de 1985  1220.5726      2017    CAQUETÁ
8    Ordenanza 3 de Noviembre 12 de 1985   462.4799      2017    CAQUETÁ
9    Ordenanza 3 de Noviembre 12 de 1985  1304.7690      2017    CAQUETÁ
10   Ordenanza 3 de Noviembre 12 de 1985  1330.2126      2017    CAQUETÁ
   Shape_Leng Shape_Area Dept_cod Departamento Mun_cod
1    2.942508 0.20692777       18      CAQUETÁ     001
2    1.112829 0.03361758       18      CAQUETÁ     029
3    2.234657 0.09674460       18      CAQUETÁ     094
4    3.154370 0.08986744       18      CAQUETÁ     247
5    3.529316 0.10030928       18      CAQUETÁ     256
6    3.402939 0.13817351       18      CAQUETÁ     410
7    1.863197 0.09912782       18      CAQUETÁ     460
8    1.518688 0.03755356       18      CAQUETÁ     479
9    2.040837 0.10589313       18      CAQUETÁ     610
10   2.313848 0.10800551       18      CAQUETÁ     860
                Municipio      NBI   Miseria  Vivienda Servicios
1               FLORENCIA 14.44902  2.774899  5.224497  1.488778
2                 ALBANIA 19.80289  4.400642  9.832684  3.208801
3  BELÉN DE LOS ANDAQUIES 28.28503  7.091305 12.143430  3.906519
4             EL DONCELLO 20.37901  3.898247 10.300478  1.473936
5               EL PAUJIL 23.33489  5.170800 13.718609  2.823273
6            LA MONTAÑITA 32.51475  8.972579 17.086081  7.150295
7                   MILÁN 46.96275 25.214900 26.934097 32.865330
8                 MORELIA 20.17910  4.477612 12.000000  3.074627
9     SAN JOSÉ DEL FRAGUA 26.27451  5.763794 13.807570  5.389877
10             VALPARAÍSO 27.38640  6.005636 12.134554  7.467418
   Hacinamiento Inasistencia Dep_económica                       geometry
1      4.434199     2.077731      4.529297 POLYGON ((850567.7 734450.4...
2      4.606922     1.260601      6.119642 POLYGON ((797693.6 642855.1...
3      5.602016     4.926108     10.654141 POLYGON ((809754.1 685340, ...
4      3.704758     2.025950      7.608696 POLYGON ((857152.4 748524, ...
5      5.186398     2.737482      5.357978 POLYGON ((856558 736458.9, ...
6      5.128428     2.724748     11.315515 POLYGON ((852387.3 687475.8...
7      6.432665     3.309456      8.896848 POLYGON ((853518.9 641899, ...
8      3.104478     1.880597      6.238806 POLYGON ((811429.8 666542.7...
9      3.748290     1.842225      8.189694 POLYGON ((767411.4 667425.3...
10     1.990137     4.209229      8.365622 POLYGON ((815926 638605.9, ...
        poverty
1  Intermediate
2  Intermediate
3          High
4  Intermediate
5          High
6          High
7       Extreme
8  Intermediate
9  Intermediate
10 Intermediate

La función propSymbolsTypoLayer() crea mapas de símbolos que son proporcionales a los variables de la primera variable y el color refleja la modamilidad de la segunda variable cualitativa.

library(sf)
library(cartography)

opar <- par(mar = c(0,0,1.2,0))

plot(st_geometry(nbi_munic_2), col="#f2efe9", border="#b38e43", bg = "#aad3df", 
     lwd = 0.5)

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 = "multi.pal", n1 = 3),
  legend.var2.pos = c(1050000, 1200000), 
  legend.var2.title.txt = "Poverty"
) 

layoutLayer(title="NBI Distribution in Caqueta", 
            author = "Luis Castillo", 
            sources = "Source: DANE, 2018", 
            scale = 1, tabtitle = TRUE, frame = TRUE)

north(pos = "topleft")

7.4 Label maps

Se combinarán 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 = "darkseagreen4", 
     bg = "grey75", lwd = 0.5)
# plot NBI
choroLayer(
  x = nbi_munic_new, 
  var = "NBI",
  method = "geom",
  nclass=5,
  col = carto.pal(pal1 = "sand.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.1, 
  overlap = FALSE, 
  show.lines = FALSE
)
# map layout
layoutLayer(
  title = "Municipalities of Caqueta", 
  sources = "Source: DANE, 2018",  
  author = "Luis Castillo", 
  frame = TRUE,
  north = TRUE, 
  tabtitle = TRUE, 
  theme = "taupe.pal"
) 

7.5. Mapas de isopletos

Los mapas isopletos se basan en la hipotesis de que el fenómeno a representar tiene una distribución continua. Estos mapas utilizan un enfoque de modelización de la interacción espacial que tiene por objeto calcular indicadores basados en valores de existencias ponderados por la distancia. Permite una representación espacial del fenómeno independiente de la heterogeneidad inicial de la división territorial.

7.6. Para caña panelera

Usemos otro conjunto de datos para hacer un mapa de la isla. En este caso, subiré datos estadísticos sobre la producción de caña panelera en caqueta en 2018. Ya conocemos este conjunto de datos ya que fue utilizado en un cuaderno anterior para ilustrar uniones basadas en atributos.

crops2018 <- read_excel("d:/Users/Janus/Documents/Geomatica basica/18_CAQUETA/Caqueta_EVA.xlsx")

crops2018

crops2018 %>%
  filter(CULTIVO == "CANA PANELERA",YEAR == 2018) -> Cana18
head(Cana18)

Cana18$MPIO_CCDGO <- as.character(Cana18$COD_MUN)
Cana_Mun = left_join(Caq_MUN, Cana18, by="MPIO_CCDGO")
head(Cana_Mun)

Simple feature collection with 6 features and 24 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -76.1027 ymin: 0.9764735 xmax: -74.89527 ymax: 2.326755
geographic CRS: WGS 84
  DPTO_CCDGO MPIO_CCDGO             MPIO_CNMBR
1         18      18001              FLORENCIA
2         18      18029                ALBANIA
3         18      18094 BELÉN DE LOS ANDAQUÍES
4         18      18247            EL DONCELLO
5         18      18256              EL PAUJIL
6         18      18410           LA MONTAÑITA
                             MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR
1       Decreto 642 de Junio 17 de 1912   2547.638      2017    CAQUETÁ
2   Ordenanza 3 de Noviembre 12 de 1985    414.122      2017    CAQUETÁ
3       Decreto 963 de Marzo 14 de 1950   1191.619      2017    CAQUETÁ
4  Decreto 1678 de Septiembre 7 de 1967   1105.803      2017    CAQUETÁ
5 Decreto 1678 de Septiembret 7 de 1967   1234.743      2017    CAQUETÁ
6         Decreto 83 de Julio 6 de 1955   1701.052      2017    CAQUETÁ
  Shape_Leng Shape_Area COD_DEP DEPARTAMENTO COD_MUN
1   2.942508 0.20692777      18      CAQUETA   18001
2   1.112829 0.03361758      18      CAQUETA   18029
3   2.234657 0.09674460      18      CAQUETA   18094
4   3.154370 0.08986744      18      CAQUETA   18247
5   3.529316 0.10030928      18      CAQUETA   18256
6   3.402939 0.13817351      18      CAQUETA   18410
               MUNICIPIO             GRUPO SUBGRUPO       CULTIVO YEAR
1              FLORENCIA OTROS PERMANENTES     CANA CANA PANELERA 2018
2                ALBANIA OTROS PERMANENTES     CANA CANA PANELERA 2018
3 BELEN DE LOS ANDAQUIES OTROS PERMANENTES     CANA CANA PANELERA 2018
4            EL DONCELLO OTROS PERMANENTES     CANA CANA PANELERA 2018
5              EL PAUJIL OTROS PERMANENTES     CANA CANA PANELERA 2018
6              MONTANITA OTROS PERMANENTES     CANA CANA PANELERA 2018
  PERIODO H_Sembrada\r\n H_Cosechada PROD REND ESTADO FISICO PRODUCCION
1    2018            258         188  752  4.0                   PANELA
2    2018            597         577 3462  6.0                   PANELA
3    2018            422         372 1488  4.0                   PANELA
4    2018            363         348 2784  8.0                   PANELA
5    2018            120          90  360  4.0                   PANELA
6    2018            319         239 1291  5.4                   PANELA
  CICLO_ CULTIVO                       geometry
1     PERMANENTE POLYGON ((-75.42074 2.19413...
2     PERMANENTE POLYGON ((-75.89506 1.36569...
3     PERMANENTE POLYGON ((-75.78705 1.74982...
4     PERMANENTE POLYGON ((-75.36167 2.32142...
5     PERMANENTE POLYGON ((-75.36691 2.21234...
6     PERMANENTE POLYGON ((-75.40404 1.76944...

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

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Marco_Geocentrico_Nacional_de_Referencia in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Marco_Geocentrico_Nacional_de_Referencia in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded 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 = "Caña panalera Production Distribution in Caqueta",
            sources = "Sources: DANE and MADR, 2018",
            author = "Luis Castillo",
            frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "orange.pal")

# north arrow
north(pos = "topleft")

par(opar)

7.7. Para cafe

Ahora vamos realizar el mismo ejercicio pero para el cultivo de cafe en caqueta para el año 2018

crops2018 %>%
  filter(CULTIVO == "CAFE",YEAR == 2018) -> Cafe18
head(Cafe18)

Cafe18$MPIO_CCDGO <- as.character(Cafe18$COD_MUN)
Cafe_Mun = left_join(Caq_MUN, Cafe18, by="MPIO_CCDGO")
head(Cafe_Mun)

Simple feature collection with 6 features and 24 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: -76.1027 ymin: 0.9764735 xmax: -74.89527 ymax: 2.326755
geographic CRS: WGS 84
  DPTO_CCDGO MPIO_CCDGO             MPIO_CNMBR
1         18      18001              FLORENCIA
2         18      18029                ALBANIA
3         18      18094 BELÉN DE LOS ANDAQUÍES
4         18      18247            EL DONCELLO
5         18      18256              EL PAUJIL
6         18      18410           LA MONTAÑITA
                             MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR
1       Decreto 642 de Junio 17 de 1912   2547.638      2017    CAQUETÁ
2   Ordenanza 3 de Noviembre 12 de 1985    414.122      2017    CAQUETÁ
3       Decreto 963 de Marzo 14 de 1950   1191.619      2017    CAQUETÁ
4  Decreto 1678 de Septiembre 7 de 1967   1105.803      2017    CAQUETÁ
5 Decreto 1678 de Septiembret 7 de 1967   1234.743      2017    CAQUETÁ
6         Decreto 83 de Julio 6 de 1955   1701.052      2017    CAQUETÁ
  Shape_Leng Shape_Area COD_DEP DEPARTAMENTO COD_MUN   MUNICIPIO
1   2.942508 0.20692777      18      CAQUETA   18001   FLORENCIA
2   1.112829 0.03361758      NA         <NA>      NA        <NA>
3   2.234657 0.09674460      NA         <NA>      NA        <NA>
4   3.154370 0.08986744      18      CAQUETA   18247 EL DONCELLO
5   3.529316 0.10030928      18      CAQUETA   18256   EL PAUJIL
6   3.402939 0.13817351      18      CAQUETA   18410   MONTANITA
              GRUPO SUBGRUPO CULTIVO YEAR PERIODO H_Sembrada\r\n
1 OTROS PERMANENTES     CAFE    CAFE 2018    2018          1.122
2              <NA>     <NA>    <NA>   NA    <NA>             NA
3              <NA>     <NA>    <NA>   NA    <NA>             NA
4 OTROS PERMANENTES     CAFE    CAFE 2018    2018        322.000
5 OTROS PERMANENTES     CAFE    CAFE 2018    2018        231.000
6 OTROS PERMANENTES     CAFE    CAFE 2018    2018        187.000
  H_Cosechada PROD REND ESTADO FISICO PRODUCCION CICLO_ CULTIVO
1       1.004 1038 1.03   CAFE VERDE EQUIVALENTE     PERMANENTE
2          NA   NA   NA                     <NA>           <NA>
3          NA   NA   NA                     <NA>           <NA>
4     304.000  472 1.55   CAFE VERDE EQUIVALENTE     PERMANENTE
5     197.000  254 1.29   CAFE VERDE EQUIVALENTE     PERMANENTE
6     177.000  153 0.86   CAFE VERDE EQUIVALENTE     PERMANENTE
                        geometry
1 POLYGON ((-75.42074 2.19413...
2 POLYGON ((-75.89506 1.36569...
3 POLYGON ((-75.78705 1.74982...
4 POLYGON ((-75.36167 2.32142...
5 POLYGON ((-75.36691 2.21234...
6 POLYGON ((-75.40404 1.76944...

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

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Marco_Geocentrico_Nacional_de_Referencia in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded datum Marco_Geocentrico_Nacional_de_Referencia in CRS definitionDiscarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preservedDiscarded 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 = "Cafe Production Distribution in Caqueta",
            sources = "Sources: DANE and MADR, 2018",
            author = "Luis Castillo",
            frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "orange.pal")

# north arrow
north(pos = "topleft")

8. Guardar los mapas

Produzcamos otro mapa de la producción de cana panelera en 201Produzcamos otro mapa de la producción de cana panelera en 2018. Esta vez usaremos símbolos proporcionales y mapas de coropletes. La salida se guardará como un archivo .png.

Para caña panelera

### open the plot
png("d:/Users/Janus/Documents/Geomatica basica/18_CAQUETA/cana_2018.png", width = 2048, height = 1526)
# set margins
opar <- par(mar = c(0,0,5,5))
# Plot the municipalities
plot(st_geometry(Cana_Proj), col="darkseagreen3", border="darkseagreen4",  
     bg = "white", lwd = 0.6)
# Plot symbols with choropleth coloration
propSymbolsChoroLayer(x = Cana_Proj, var = "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 = 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 = Cana_Proj, 
  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="Caña panelera Production & Yield in Caqueta, 2018",
            author = "Luis Castillo", 
            sources = "Sources: MADR & DANE, 2018", 
            scale = 50, tabtitle = FALSE, frame = TRUE)
# north arrow
north(pos = "topleft")
#
title(main="Caña panelera Production & Yield in caqueta, 2018", cex.main=3,
      sub= "Source: MADR & DANE, 2018", cex.sub=2)
#
graticule = TRUE
#
par(opar)
### close the plot
dev.off()

null device 
          1 

knitr::include_graphics("d:/Users/Janus/Documents/Geomatica basica/18_CAQUETA/cana_2018.png")

8.1. para cafe

Produzcamos otro mapa de la producción de cana panelera en 2018.

### open the plot
png("d:/Users/Janus/Documents/Geomatica basica/18_CAQUETA/cafe_2018.png", width = 2048, height = 1526)
# set margins
opar <- par(mar = c(0,0,5,5))
# Plot the municipalities
plot(st_geometry(Cafe_Proj), col="darkseagreen3", border="darkseagreen4",  
     bg = "white", lwd = 0.6)
# Plot symbols with choropleth coloration
propSymbolsChoroLayer(x = Cafe_Proj, var = "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 = 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 = Cafe_Proj, 
  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="Cafe Production & Yield in Caqueta, 2018",
            author = "Luis Castillo", 
            sources = "Sources: MADR & DANE, 2018", 
            scale = 50, tabtitle = FALSE, frame = TRUE)
# north arrow
north(pos = "topleft")
#
title(main="Cafe Production & Yield in caqueta, 2018", cex.main=3,
      sub= "Source: MADR & DANE, 2018", cex.sub=2)
#
graticule = TRUE
#
par(opar)
### close the plot
dev.off()

null device 
          1 

knitr::include_graphics("d:/Users/Janus/Documents/Geomatica basica/18_CAQUETA/cafe_2018.png")

9. BIBLIOGRAFÍA

CDC.2018. Thypes of Thematic Maps. Diponible en CDC

Bonilla, J. 2011. Sistemas de Información Geográfica. Disponible en Volaya

Rodríguez, L. 2019. Capítulo 1. Introducción a la Cartografía Temática. Disponible en Red Geomática