Thematic Cartography
Ivan Lizarazo
6 Abril 2020
1. Why this notebook?
This is an R Markdown Notebook illustrating thematic cartography for a Departament in Colombia. It has been created, compiled and published from RStudio Cloud.
This notebook aims to help Geomatica Basica students at Universidad Nacional to learn how to make thematic maps in R. Every student should replicate this notebook BUT adapting its contents to her/his department.
The due date for publishing this notebook is on 20th May 2020 at 11:59 am. Its publication at RPubs counts as the second technical report for the course. Its mark represents 15% of the final grade.
2. Thematic cartography
As Amanda Briney states, a thematic map emphasizes a theme or topic, such as the average distribution of rainfall in an area or the population density in a municipality. They’re different from general reference maps because they don’t just show natural and human-made features such as rivers, cities, political subdivisions, and highways. If these items appear on a thematic map, they’re reference points to enhance one’s understanding of the map’s theme and purpose.
According to Briney, “the most significant factor to consider when designing thematic maps is the map’s audience, which helps determine what items should be included on the map as reference points in addition to the theme. A map being made for a political scientist, for example, would need to show political boundaries, whereas one for a biologist might need contours showing elevation”.
The sources of thematic maps’ data are very important, emphasizes Briney. “Cartographers must find accurate, recent, reliable sources of information on a wide range of subjects, from environmental features to demographic data, to make the best possible maps”.
Briney summarizes several thematic mapping techniques which are used most often:
2.1 Choropleth map
which portrays quantitative data as a color and can show density, percent, average value, or quantity of an event within a geographic area. Sequential colors represent increasing or decreasing positive or negative data values. Normally, each color also represents a range of values.
2.2 Proportional or graduated symbols map
which are used in another type of map to represent data associated with locations, such as cities. Data is displayed on these maps with proportionally sized symbols to show differences in occurrences. Circles are most often used, but squares and other geometric shapes are also suitable. The most common way to size these symbols is to make their areas proportional to the values to be depicted using mapping or drawing software.
2.3 Isarithmic or contour map
which uses isolines to depict continuous values such as precipitation levels. These maps also can display three-dimensional values, such as elevation, on topographic maps. Generally, data for isarithmic maps is gathered via measurable points (e.g. weather stations) or is collected by area (e.g. tons of corn per acre by county). Isarithmic maps also follow the basic rule that there are high and low sides in relation to the isoline. For example, in elevation, if the isoline is 500 feet, then one side must be higher than 500 feet and one side must be lower.
2.4 Dot map
which uses dots to show the presence of a theme and display a spatial pattern. A dot can represent one unit or several, depending on what is being depicted.
3. Data
We will use data on Necesidades Basicas Insatisfechas (NBI) from the Censo Nacional de Población y Vivienda 2018 which are available from the DANE Geoportal.
Thematic maps are useful to convey demographic information. You may explore several thematic maps from DANE using this link.
Previously, I downloaded the NBI file, in format .xlsx, to my computer. Then, I used excel to remove data for municipalities not located within the Antioquia Department. I also “cleaned” the data, that it, I removed several rows with either institutional images or no information. These rows were located both at the beginning and the ending of the original file. I keep the columns referring to the whole municipality (i.e. I removed the columms corresponding to “cabecera” and “zona rural”). I saved the resulting data, corresponding only to Antioquian municipalities as NBI_Antioquia.xlsx. Then, I uploaded the file to the directory “ndvi” under my project’s folder at RStudio Cloud.
4. Preparation
Let’s clean the memory:
rm(list=ls())Now, let’s install the libraries we need.
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)Then, let’s load the libraries.
library(tidyverse)## ── Attaching packages ────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.0 ✓ purrr 0.3.3
## ✓ tibble 2.1.3 ✓ dplyr 0.8.5
## ✓ tidyr 1.0.2 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0## ── Conflicts ───────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(readxl) ## this library is part of the tidyverse
library(rgeos)## Loading required package: sp## rgeos version: 0.5-2, (SVN revision 621)
## GEOS runtime version: 3.5.1-CAPI-1.9.1
## Linking to sp version: 1.4-1
## Polygon checking: TRUElibrary(raster)##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## select## The following object is masked from 'package:tidyr':
##
## extractlibrary(sf)## Linking to GEOS 3.5.1, GDAL 2.2.2, PROJ 4.9.2library(cartography)
library(SpatialPosition)5. Read NBI data
Let’s read the csv file with “estadisticas municipales agropecuarias” for Antioquia.
nbi <- read_excel("./nbi/NBI_Antioquia.xlsx")Let’s check what are the attributes of datos.
head(nbi)## # A tibble: 6 x 12
## COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO NBI MISERIA VIVIENDA SERVICIOS
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 05 ANTI… 001 05001 MEDELLÍN 5.19 0.412 0.285 0.189
## 2 05 ANTI… 002 05002 ABEJORRAL 14.0 1.38 0.376 0.495
## 3 05 ANTI… 004 05004 ABRIAQUÍ 11.9 1.21 1.73 1.59
## 4 05 ANTI… 021 05021 ALEJANDR… 12.7 0.889 0.706 0.0784
## 5 05 ANTI… 030 05030 AMAGÁ 9.66 0.948 0.528 0.277
## 6 05 ANTI… 031 05031 AMALFI 22.3 4.86 8.89 2.98
## # … with 3 more variables: HACINAMIENTO <dbl>, INASISTENCIA <dbl>,
## # ECONOMIA <dbl>Let’s find which is the muncipality with the highest percent of NBI:
nbi %>%
slice(which.max(NBI)) -> max_nbi
max_nbi## # A tibble: 1 x 12
## COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO NBI MISERIA VIVIENDA SERVICIOS
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 05 ANTI… 475 05475 MURINDÓ 81.7 34.2 29.6 60.3
## # … with 3 more variables: HACINAMIENTO <dbl>, INASISTENCIA <dbl>,
## # ECONOMIA <dbl>Let’s find which is the muncipality with the lowest percent of NBI:
nbi %>%
slice(which.min(NBI)) -> min_nbi
min_nbi## # A tibble: 1 x 12
## COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO NBI MISERIA VIVIENDA SERVICIOS
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 05 ANTI… 631 05631 SABANETA 1.59 0.0723 0.0502 0.00613
## # … with 3 more variables: HACINAMIENTO <dbl>, INASISTENCIA <dbl>,
## # ECONOMIA <dbl>Let’s sort municipalities by NBI in descending order:
nbi %>%
arrange(desc(NBI)) -> desc_nbi
desc_nbi## # A tibble: 125 x 12
## COD_DEPTO DEPTO COD_MUN CODIGO MUNICIPIO NBI MISERIA VIVIENDA SERVICIOS
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 05 ANTI… 475 05475 MURINDÓ 81.7 34.2 29.6 60.3
## 2 05 ANTI… 873 05873 VIGÍA DE… 77.8 20.1 1.74 71.8
## 3 05 ANTI… 665 05665 SAN PEDR… 66.0 33.6 61.4 21.7
## 4 05 ANTI… 051 05051 ARBOLETES 62.4 23.0 58.1 11.4
## 5 05 ANTI… 659 05659 SAN JUAN… 59.8 28.9 52.3 22.5
## 6 05 ANTI… 490 05490 NECOCLÍ 57.4 28.2 47.7 22.0
## 7 05 ANTI… 495 05495 NECHÍ 54.2 20.5 47.3 14.7
## 8 05 ANTI… 120 05120 CÁCERES 49.8 18.4 38.0 8.92
## 9 05 ANTI… 234 05234 DABEIBA 47.2 28.6 24.3 28.4
## 10 05 ANTI… 895 05895 ZARAGOZA 45.7 16.7 33.4 11.0
## # … with 115 more rows, and 3 more variables: HACINAMIENTO <dbl>,
## # INASISTENCIA <dbl>, ECONOMIA <dbl>6. Joining NBI data to municipalities
I have already uploaded to RStudio Cloud the municipalities within Antioquia.
Let’s read the data using the sf library:
munic <- st_read("./antioquia/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")## Reading layer `MGN_MPIO_POLITICO' from data source `/cloud/project/antioquia/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
## CRS: 4326Let’s check what is inside the attribute MPIO_CCDGO:
head(munic$MPIO_CNMBR)## [1] MEDELLÍN CIUDAD BOLÍVAR BRICEÑO BURITICÁ CÁCERES
## [6] BETULIA
## 125 Levels: ABEJORRAL ABRIAQUÍ ALEJANDRÍA AMAGÁ AMALFI ANDES ... ZARAGOZAWe can use the left_join function to join the municipalities and the NBI data.
nbi_munic = left_join(munic, nbi, by=c("MPIO_CCDGO"="CODIGO"))## Warning: Column `MPIO_CCDGO`/`CODIGO` joining factor and character vector,
## coercing into character vectornbi_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.09953 ymin: 5.753381 xmax: -74.99361 ymax: 7.957635
## CRS: 4326
## MUNICIPIO 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...Note that Caceres has NBI=49.84, a value that can be checked against the ordered NBI data previously obtained.
Now, let’s reproject the municipalities:
nbi_munic_new <- st_transform(nbi_munic, crs = 3116)7. Examples of thematic maps
We will use the cartography package which aims to obtain thematic maps with the visual quality of those build with a classical mapping or GIS software. Its functions are intuitive to users and ensure compatibility with common R workflows.
The cartography library uses sf or sp objects to produce base graphics. As most of the internals of the package relies on sf functionalities, the preferred format for spatial objects is sf.
7. 1 OpenStreetMap Basemap and Proportional Symbols
The functions getTiles() and tilesLayer() download and display OpenStreetMap tiles. Be careful to cite the source of the tiles appropriately.
propSymbolsLayer() displays symbols with areas proportional to a quantitative variable (e.g. NBI). Several symbols are available (circles, squares, bars). The inches argument is used to customize the symbols sizes.
# download osm tiles
mun.osm <- getTiles(
x = nbi_munic_new,
type = "OpenStreetMap",
zoom = 8,
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 Antioquia",
sources = "Sources: DANE, 2018\n© OpenStreetMap",
author = " Ivan Lizarazo ",
frame = TRUE, north = FALSE, tabtitle = TRUE)
# north arrow
north(pos = "topleft")
#
#par(opar)7. 2 Choropleth Maps
In choropleth maps, areas are shaded according to the variation of a quantitative variable. They are used to represent ratios or indices.
choroLayer() displays choropleth maps . Arguments nclass, method and breaks allow to customize the variable classification.
getBreaks() allow to classify outside of the function itself. Colors palettes are defined with col and a set of colors can be created with carto.pal(). See also display.carto.all().
# 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 = "sand.pal", n1 = 5),
border = "white",
lwd = 0.5,
legend.pos = "topright",
legend.title.txt = "NBI",
add = TRUE
)
# layout
layoutLayer(title = "NBI Distribution in Antioquia",
sources = "Source: DANE, 2018",
author = "Ivan Lizarazo",
frame = TRUE, north = TRUE, tabtitle = TRUE, col="black")
# north arrow
north(pos = "topleft")
#
#par(opar)7.3 Proportional Symbols and Typology Map
The function propSymbolsTypoLayer() creates a map of symbols that are proportional to values of a first variable and colored to reflect the modalities of a second qualitatice variable. A combination of propSymbolsLayer() and typoLayer() arguments is used.
First, we need to create a qualitative variable. Let’s use the mutate function for this task.
nbi_munic_2 <- dplyr::mutate(nbi_munic_new, poverty = ifelse(MISERIA > 20, "Extreme",
ifelse(HACINAMIENTO > 5, "High", "Intermediate")))Note that the new attribute is named poverty. Its value depends on threshold values defined by the user. Make sure to check the sintaxis of the ifelse command to understand what is going on.
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
## CRS: EPSG:3116
## DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR
## 1 05 05001 MEDELLÍN
## 2 05 05101 CIUDAD BOLÍVAR
## 3 05 05107 BRICEÑO
## 4 05 05113 BURITICÁ
## 5 05 05120 CÁCERES
## 6 05 05093 BETULIA
## MPIO_CRSLC MPIO_NAREA MPIO_NANO
## 1 1965 374.8280 2017
## 2 1869 260.4461 2017
## 3 Ordenanza 27 de Noviembre 26 de 1980 376.3468 2017
## 4 1812 355.2103 2017
## 5 Decreto departamental 160 del 16 de Marzo de 1903 1873.8033 2017
## 6 Decreto departamental 629 del 28 de Enero de 1884 262.3675 2017
## DPTO_CNMBR Shape_Leng Shape_Area COD_DEPTO DEPTO COD_MUN MUNICIPIO
## 1 ANTIOQUIA 1.0327835 0.03060723 05 ANTIOQUIA 001 MEDELLÍN
## 2 ANTIOQUIA 0.7085039 0.02124224 05 ANTIOQUIA 101 CIUDAD BOLÍVAR
## 3 ANTIOQUIA 1.0044720 0.03078496 05 ANTIOQUIA 107 BRICEÑO
## 4 ANTIOQUIA 0.9637233 0.02902757 05 ANTIOQUIA 113 BURITICÁ
## 5 ANTIOQUIA 2.9333643 0.15350440 05 ANTIOQUIA 120 CÁCERES
## 6 ANTIOQUIA 0.8476756 0.02141352 05 ANTIOQUIA 093 BETULIA
## NBI MISERIA VIVIENDA SERVICIOS HACINAMIENTO INASISTENCIA ECONOMIA
## 1 5.192345 0.4119415 0.2846426 0.1888295 1.536231 1.600220 2.039155
## 2 11.708395 1.2648358 0.4461578 0.5327904 2.395391 1.416443 8.334055
## 3 25.327952 5.1463169 3.6158762 3.3636058 7.484023 2.556340 14.278507
## 4 28.908795 6.5146580 0.4614549 3.2844734 11.644951 6.460369 14.495114
## 5 49.846861 18.4341501 38.0474732 8.9165391 7.392802 3.346095 14.873660
## 6 17.127878 2.3491937 1.1281439 0.8096091 4.532484 2.455372 10.916451
## geometry poverty
## 1 POLYGON ((823819.2 1196825,... Intermediate
## 2 POLYGON ((782138 1147634, 7... Intermediate
## 3 POLYGON ((847501.8 1290702,... High
## 4 POLYGON ((797631.6 1263319,... High
## 5 POLYGON ((875837.6 1371851,... High
## 6 POLYGON ((786889.9 1186783,... Intermediatelibrary(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 = "multi.pal", n1 = 3),
legend.var2.pos = c(1050000, 1200000),
legend.var2.title.txt = "Poverty"
)
# layout
layoutLayer(title="NBI Distribution in Antioquia",
author = "Ivan Lizarazo",
sources = "Source: DANE, 2018",
scale = 1, tabtitle = TRUE, frame = TRUE)
# north arrow
north(pos = "topleft")
#
#par(opar)7.4 Label maps
Let’s try combining the functions choroLayer and 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 Antioquia",
sources = "Source: DANE, 2018",
author = "Ivan Lizarazo",
frame = TRUE,
north = TRUE,
tabtitle = TRUE,
theme = "taupe.pal"
) 
#
#par(opar)7.5 Isopleth Maps
Isopleth maps are based on the assumption that the phenomenon to be represented has a continuous distribution. These maps use a spatial interaction modeling approach which aims to compute indicators based on stock values weighted by distance. It allows a spatial representation of the phenomenon independent from the initial heterogeneity of the territorial division. smoothLayer() heavily depends on the SpatialPosition package. The function uses a marked point layer and a set of parameters (a spatial interaction function and its parameters) and displays an isopleth map layer.
Let’s use another dataset to make an isopleth map. In this case, I will upload statistical data on 2018’s coffee production in Antioquia. We already know this dataset as it was used in a previous notebook to illustrate attribute-based joins.
Read the dataset:
crops2018 <- read_excel("./agro/EVA_Antioquia_2.xlsx")head(crops2018) ## # A tibble: 6 x 14
## COD_DEP DEPARTAMENTO C0D_MUN MUNICIPIO GRUPO SUBGRUPO CULTIVO YEAR
## <dbl> <chr> <dbl> <chr> <chr> <chr> <chr> <dbl>
## 1 5 ANTIOQUIA 5002 ABEJORRAL FRUT… AGUACATE AGUACA… 2018
## 2 5 ANTIOQUIA 5002 ABEJORRAL OTRO… CA√ëA CA√ëA … 2018
## 3 5 ANTIOQUIA 5002 ABEJORRAL OTRO… CACAO CACAO 2018
## 4 5 ANTIOQUIA 5002 ABEJORRAL OTRO… CAFE CAFE 2018
## 5 5 ANTIOQUIA 5002 ABEJORRAL LEGU… FRIJOL FRIJOL 2018
## 6 5 ANTIOQUIA 5002 ABEJORRAL FLOR… FLORES HORTEN… 2018
## # … with 6 more variables: Area_Siembra <dbl>, Area_Cosecha <dbl>,
## # Produccion <dbl>, Rendimiento <dbl>, ESTADO <chr>, CICLO <chr>Filter rows which represent only coffee data:
crops2018 %>%
filter(CULTIVO == "CAFE") -> cafe2018Create a new attribute which matches municipalities codes:
cafe2018$TEMP <- as.character(cafe2018$C0D_MUN)cafe2018$MPIO_CCDGO <- as.factor(paste(0, cafe2018$TEMP, sep=""))Now, make the join:
cafe_munic = left_join(munic, cafe2018, by="MPIO_CCDGO")## Warning: Column `MPIO_CCDGO` joining factors with different levels, coercing to
## character vectorCheck the output:
head(cafe_munic)## Simple feature collection with 6 features and 24 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -76.09953 ymin: 5.753381 xmax: -74.99361 ymax: 7.957635
## CRS: 4326
## DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR
## 1 05 05001 MEDELLÍN
## 2 05 05101 CIUDAD BOLÍVAR
## 3 05 05107 BRICEÑO
## 4 05 05113 BURITICÁ
## 5 05 05120 CÁCERES
## 6 05 05093 BETULIA
## MPIO_CRSLC MPIO_NAREA MPIO_NANO
## 1 1965 374.8280 2017
## 2 1869 260.4461 2017
## 3 Ordenanza 27 de Noviembre 26 de 1980 376.3468 2017
## 4 1812 355.2103 2017
## 5 Decreto departamental 160 del 16 de Marzo de 1903 1873.8033 2017
## 6 Decreto departamental 629 del 28 de Enero de 1884 262.3675 2017
## DPTO_CNMBR Shape_Leng Shape_Area COD_DEP DEPARTAMENTO C0D_MUN MUNICIPIO
## 1 ANTIOQUIA 1.0327835 0.03060723 5 ANTIOQUIA 5001 MEDELLIN
## 2 ANTIOQUIA 0.7085039 0.02124224 5 ANTIOQUIA 5101 CIUDAD BOLIVAR
## 3 ANTIOQUIA 1.0044720 0.03078496 5 ANTIOQUIA 5107 BRICEÑO
## 4 ANTIOQUIA 0.9637233 0.02902757 5 ANTIOQUIA 5113 BURITICA
## 5 ANTIOQUIA 2.9333643 0.15350440 NA <NA> NA <NA>
## 6 ANTIOQUIA 0.8476756 0.02141352 5 ANTIOQUIA 5093 BETULIA
## GRUPO SUBGRUPO CULTIVO YEAR Area_Siembra Area_Cosecha Produccion
## 1 OTROS PERMANENTES CAFE CAFE 2018 471 433 299
## 2 OTROS PERMANENTES CAFE CAFE 2018 10063 7374 11439
## 3 OTROS PERMANENTES CAFE CAFE 2018 589 485 586
## 4 OTROS PERMANENTES CAFE CAFE 2018 936 861 964
## 5 <NA> <NA> <NA> NA NA NA NA
## 6 OTROS PERMANENTES CAFE CAFE 2018 5734 4230 6562
## Rendimiento ESTADO CICLO TEMP
## 1 0.69 CAFE VERDE EQUIVALENTE PERMANENTE 5001
## 2 1.55 CAFE VERDE EQUIVALENTE PERMANENTE 5101
## 3 1.21 CAFE VERDE EQUIVALENTE PERMANENTE 5107
## 4 1.12 CAFE VERDE EQUIVALENTE PERMANENTE 5113
## 5 NA <NA> <NA> <NA>
## 6 1.55 CAFE VERDE EQUIVALENTE PERMANENTE 5093
## geometry
## 1 POLYGON ((-75.66974 6.37359...
## 2 POLYGON ((-76.04467 5.92774...
## 3 POLYGON ((-75.45818 7.22284...
## 4 POLYGON ((-75.90857 6.97378...
## 5 POLYGON ((-75.20358 7.95716...
## 6 POLYGON ((-76.00304 6.28171...Now, let’s reproject the municipalities:
rep_cafe <- st_transform(cafe_munic, crs = 3116)Mapping time:
# 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 = 'Produccion',
typefct = "exponential",
span = 25000,
beta = 2,
nclass = 10,
col = carto.pal(pal1 = 'green.pal', n1 = 10),
border = "grey",
lwd = 0.1,
mask = rep_cafe,
legend.values.rnd = -3,
legend.title.txt = "Production",
legend.pos = "topright",
add=TRUE
)
# 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 = "Coffee Production Distribution in Antioquia",
sources = "Sources: DANE and MADR, 2018",
author = "Ivan Lizarazo",
frame = FALSE, north = FALSE, tabtitle = TRUE, theme = "green.pal")
# north arrow
north(pos = "topleft")
#
par(opar)8. Saving maps
Let’s produce another map of coffee production in 2018. This time we will use proportional symbols and choropleth maps. The output will be saved as a .png file.
First, a few remarks: - propSymbolsChoroLayer() creates a map of symbols that are proportional to values of a first variable and colored to reflect the classification of a second variable. - a combination of propSymbolsLayer() and choroLayer() arguments is used.
The following chunk does not display a map. Instead, it writes the map to the filename cafe_2018.png under the work directory.
### open the plot
png("./agro/cafe_2018.png", width = 2048, height = 1526)
# set margins
opar <- par(mar = c(0,0,5,5))
# Plot the municipalities
plot(st_geometry(rep_cafe), col="darkseagreen3", border="darkseagreen4",
bg = "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="Coffee Production & Yield in Antioquia, 2018",
author = "Ivan Lizarazo",
sources = "Sources: MADR & DANE, 2018",
scale = 50, tabtitle = FALSE, frame = TRUE)
# north arrow
north(pos = "topleft")
#
title(main="Coffee Production & Yield in Antioquia, 2018", cex.main=3,
sub= "Source: MADR & DANE, 2018", cex.sub=2)
#
graticule = TRUE
#
par(opar)
### close the plot
dev.off()## png
## 2Once we saved a map, we can add it to our R Markdown report using markdown syntax as follows:
Coffe Production in Antioquia, 2018
9. Your turn
It’s your turn now to create your own notebook on thematic cartography. Use your department datasets to illustrate how to make five different types of thematic maps. Save at least one high-quality thematic map and display it within your notebook.
Feel free to use any municipalities data you consider relevant for your academic interests (e.g. evaluaciones agricolas, poblacion, or NBI data).
For completing your task, please have a look at the following link on cartography functionalities.
Good luck!!!
sessionInfo()## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.6 LTS
##
## Matrix products: default
## BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
## LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] SpatialPosition_2.0.0 cartography_2.4.0 sf_0.9-0
## [4] raster_3.0-12 rgeos_0.5-2 sp_1.4-1
## [7] readxl_1.3.1 forcats_0.5.0 stringr_1.4.0
## [10] dplyr_0.8.5 purrr_0.3.3 readr_1.3.1
## [13] tidyr_1.0.2 tibble_2.1.3 ggplot2_3.3.0
## [16] tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.3 lubridate_1.7.4 lattice_0.20-38 png_0.1-7
## [5] class_7.3-15 utf8_1.1.4 assertthat_0.2.1 digest_0.6.25
## [9] R6_2.4.1 cellranger_1.1.0 backports_1.1.5 reprex_0.3.0
## [13] evaluate_0.14 e1071_1.7-3 httr_1.4.1 pillar_1.4.3
## [17] rlang_0.4.5 rstudioapi_0.11 rmarkdown_2.1 rgdal_1.4-8
## [21] munsell_0.5.0 broom_0.5.5 compiler_3.6.0 modelr_0.1.6
## [25] xfun_0.12 pkgconfig_2.0.3 htmltools_0.4.0 tidyselect_1.0.0
## [29] codetools_0.2-16 fansi_0.4.1 crayon_1.3.4 dbplyr_1.4.2
## [33] withr_2.1.2 grid_3.6.0 lwgeom_0.2-1 nlme_3.1-139
## [37] jsonlite_1.6.1 gtable_0.3.0 lifecycle_0.2.0 DBI_1.1.0
## [41] magrittr_1.5 units_0.6-6 scales_1.1.0 KernSmooth_2.23-15
## [45] cli_2.0.2 stringi_1.4.6 fs_1.3.2 xml2_1.2.5
## [49] generics_0.0.2 vctrs_0.2.4 tools_3.6.0 glue_1.3.2
## [53] hms_0.5.3 slippymath_0.3.1 yaml_2.2.1 colorspace_1.4-1
## [57] isoband_0.2.0 classInt_0.4-2 rvest_0.3.5 knitr_1.28
## [61] haven_2.2.0