My second R Notebook: How to make thematic maps

1.Introduction

This is the second R Markdown Notebook for the Geomatica Basica 2022 course. It illustrates how to make thematic maps showing the municipal share of the two most important groups of crops for a given department. We will use as main source the csv files saved in the EVA notebook, as well as a shapefile of municipalities obtained in class.

2. Setup

In this step, we will install and load the required R libraries.

#run the following lines from the command line:
#install.packages('dplyr')
#install.packages('readxl')
#install.packages('sf)

library(sf)
library(dplyr)
library(readr)

3. Read files related to municipios, crops and cities

I will work with the Departamento de Santander. As you work in a different computer and you are interested in a different department, you may need some changes to file paths, file names or variables names.

Before starting, make sure that you have saved the following files in your computer:

• a shapefile with your department’s municipalities (as produced in class)
• a csv file with EVA 2020 statistics for the selected group of crops (as produced in the first notebook)

In addition, you will need a file with cities in Colombia. You may download it from https://simplemaps.com/data/co-cities. Then, move the downloaded co.csv file to your working directory.

Let’s start by listing the available files.

What are the csv files saved in the working directory?

list.files("./datos", pattern=c('csv'))

## [1] "co.csv"                  "stder_frutales_2020.csv"
## [3] "stder_oleag_2020.csv"

What are the shapefiles saved in the working directory?

#change the following line  to go to your data directory
list.files('./68_SANTANDER')

## [1] "STDER_MPIO.cpg"     "STDER_MPIO.dbf"     "STDER_MPIO.prj"    
## [4] "STDER_MPIO.sbn"     "STDER_MPIO.sbx"     "STDER_MPIO.shp"    
## [7] "STDER_MPIO.shp.xml" "STDER_MPIO.shx"

Now, let’s proceed to read the EVA-based files obtained from the first notebook:

# adjust the filepath acording to your data
(frutales = read_csv("./datos/stder_frutales_2020.csv",show_col_types = FALSE))

The, read the municipalities of your department:

# this shapefile must have a 4326 EPSG code
# this shapefile was obtained in class using QGIS
(mun.tmp =  st_read('./68_SANTANDER/STDER_MPIO.shp'))

## Reading layer `STDER_MPIO' from data source 
##   `/Users/ivanlizarazo/Documents/GB/68_SANTANDER/STDER_MPIO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 87 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
## Geodetic CRS:  WGS 84

Let’s select a few attributes to clean the object:

# check the objet output in the last chunk and
# change attribute names according to your own data
mun.tmp %>% select(MPIO_CCDGO, MPIO_CNMBR, MPIO_NAREA) -> municipios

municipios

Now, read the cities csv file:

# This file was downloaded from simplemaps as described above
(cities = read_csv("./datos/co.csv"))

## Rows: 1102 Columns: 9

## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): city, country, iso2, admin_name, capital
## dbl (4): lat, lng, population, population_proper

## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Note that there are 1102 rows representing cities en Colombia. Note that the cities object is a non-spatial object (a tibble, actually).

Note also that coordinates values are stored in the “lng” (latitude) and “lat” (longitude) attributes.

To convert cities to a spatial object (a simple feature object, in this case):

# Convert data frame to simple feature object
sf.cities <-  st_as_sf(x = cities, 
                        coords = c("lng", "lat"))

sf.cities

Note the missing CRS!

Let’s add the CRS by using the EPSG code:

st_crs(sf.cities) <- 4326

4. Subset data relevant to our department

As we are interested in just a single department we need to create an spatial join:

# find points (cities) within polygons (our municipalities)
sf.cities.joined <- st_join(sf.cities, municipios, join = st_within)

To display the joined object:

sf.cities.joined

Note that we obtained an sf data frame with each row of cities appended with the columns from our municipios. Cities located in a different department have a lot of NAs.

Now, let’s filter the rows which correspond to our deparment (in this case Santander):

stder.cities = dplyr::filter(sf.cities.joined, admin_name=='Santander')

stder.cities 

We obtained 87 cities located in Santander.

5. Make a map for the most important group of crops

Now, make a choropleth map of this data. We will use the tmap, ggplot2, and classInt libraries.

# run the following lines from the command window:
#install.packages("tmap")
#install.packages("ggplot2")
#install.packages("classInt")

library(tmap)

## Warning: package 'tmap' was built under R version 4.1.2

library(ggplot2)
library(ggrepel)
library(classInt)

Let’s remind that the municipios object do not have EVA attributes. In contrast, the frutales objects, which contains crop statistics, is a non-spatial object. Our next task is to join the stats object to the spatial objects to have the relevant data in a single object.

To be able to make the join, we need a shared key, i.e. a common attribute. In this case, we have it on both objects, that is the codigo municipal. However, note that both their names and their data types are different.

### check according to your own data
class(frutales$Cod_Mun)

## [1] "numeric"

class(municipios$MPIO_CCDGO)

## [1] "character"

Thus, we need to fix this:

frutales$Cod_Mun = as.character(frutales$Cod_Mun)

Now, we are ready to make a join. You may need to review your class notes for understanding what is going here:

munic_frutales = left_join(municipios, frutales, by = c("MPIO_CCDGO" = "Cod_Mun"))

What is the result?

munic_frutales

Note that munic_frutales includes now the max_prod attribute represenying the total amount of tons produced by fruit crops in 2020.

The following code customizes the colors for the municipalities. More information [here] (https://stackoverflow.com/questions/57177312/trying-to-plot-in-tmap-shapefile-with-attribute).

breaks <- classIntervals(munic_frutales$max_prod, n = 6, style = 'fisher')

#label breaks
lab_vec <- vector(length = length(breaks$brks)-1)
rounded_breaks <- round(breaks$brks,2)
lab_vec[1] <- paste0('[', rounded_breaks[1],' - ', rounded_breaks[2],']')
for(i in 2:(length(breaks$brks) - 1)){
  lab_vec[i] <- paste0('(',rounded_breaks[i], ' - ', rounded_breaks[i+1], ']')
}

Now, the static map:

munic_frutales <-  munic_frutales %>%
  mutate(faktor_class = factor(cut(max_prod, breaks$brks, include.lowest = T), labels = lab_vec))
# change attribute name
munic_frutales$Produccion = munic_frutales$faktor_class

# This code creates a new field ("mid") with centroid coordinates
munic_frutales$mid <- sf::st_centroid(munic_frutales$geometry)

# Get the longitude values
LONG = st_coordinates(munic_frutales$mid)[,1]

# Get the latitude values
LAT = st_coordinates(munic_frutales$mid)[,2]

ggplot(data = munic_frutales) +
   geom_sf(aes(fill = Produccion)) +
   geom_label_repel(aes(x = LONG, y = LAT, label = MPIO_CNMBR), 
                    label.padding =     unit(0.05,"lines"),  
                    label.r = unit(0.025, "lines"),
                    label.size = 0.05)

Another option is to produce an interactive map:

6. Make a new map for the second largest group of crops

Here, we will use EVA data for the “oleaginosas y hortalizas” group which accounts for the second largest production in 2020 in Santander.

We need to conduct similar steps to the previous ones:

(oleag = read_csv("./datos/stder_oleag_2020.csv",show_col_types = FALSE))

oleag$Cod_Mun = as.character(oleag$Cod_Mun)

#municipios$MPIO_CCDGO = as.character(municipios$MPIO_CCDGO)

munic_oleag = left_join(municipios, oleag, by = c("MPIO_CCDGO" = "Cod_Mun"))

munic_oleag

# see example at https://r-tmap.github.io/tmap-book/layout.html
facet = "max_prod"
oleag_map =  
  tm_shape(munic_oleag) + tm_polygons(facet) + tm_text(text = "MPIO_CNMBR", size = 0.7, fontfamily = "sans") +
  tm_shape(stder.cities) + tm_symbols(shape = 2, col = "red", size = 0.20) +
  tm_credits("Data source: UPRA (2020)", fontface = "bold") +
  tm_layout(main.title = "Produccion de oleaginosas en 2020",
            main.title.fontface = "bold.italic", 
            legend.title.fontfamily = "monospace") +
  tm_scale_bar(position = c("left", "bottom"))

tmap_mode("view")

## tmap mode set to interactive viewing

oleag_map

## Credits not supported in view mode.

## Symbol shapes other than circles or icons are not supported in view mode.

Now, it’s time for you to try your own maps. Several map examples are available at: https://geocompr.robinlovelace.net/adv-map.html

7. Bibliography

In your notebook, cite this work as Lizarazo, I., 2022. Getting started with thematic maps. Available at https://rpubs.com/ials2un/thematic_maps_v2.

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] classInt_0.4-3 ggrepel_0.9.1  ggplot2_3.3.5  tmap_3.3-3     readr_2.1.2   
## [6] dplyr_1.0.7    sf_1.0-7      
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.0              bit64_4.0.5             vroom_1.5.7            
##  [4] jsonlite_1.7.2          viridisLite_0.4.0       bslib_0.2.5.1          
##  [7] assertthat_0.2.1        highr_0.9               sp_1.4-5               
## [10] yaml_2.2.1              pillar_1.6.2            lattice_0.20-44        
## [13] glue_1.4.2              digest_0.6.27           RColorBrewer_1.1-2     
## [16] leaflet.providers_1.9.0 colorspace_2.0-2        htmltools_0.5.1.1      
## [19] XML_3.99-0.9            pkgconfig_2.0.3         raster_3.4-13          
## [22] stars_0.5-5             s2_1.0.7                purrr_0.3.4            
## [25] scales_1.1.1            tzdb_0.2.0              tibble_3.1.3           
## [28] proxy_0.4-26            farver_2.1.0            generics_0.1.0         
## [31] ellipsis_0.3.2          withr_2.4.2             leafsync_0.1.0         
## [34] cli_3.0.1               mime_0.11               magrittr_2.0.1         
## [37] crayon_1.4.1            evaluate_0.14           fansi_0.5.0            
## [40] lwgeom_0.2-8            class_7.3-19            tools_4.1.0            
## [43] hms_1.1.1               lifecycle_1.0.0         stringr_1.4.0          
## [46] munsell_0.5.0           compiler_4.1.0          jquerylib_0.1.4        
## [49] e1071_1.7-8             rlang_0.4.11            units_0.7-2            
## [52] grid_4.1.0              tmaptools_3.1-1         dichromat_2.0-0        
## [55] rstudioapi_0.13         htmlwidgets_1.5.3       crosstalk_1.1.1        
## [58] leafem_0.1.6            base64enc_0.1-3         rmarkdown_2.9          
## [61] wk_0.5.0                gtable_0.3.0            codetools_0.2-18       
## [64] abind_1.4-5             DBI_1.1.1               markdown_1.1           
## [67] R6_2.5.0                knitr_1.33              bit_4.0.4              
## [70] utf8_1.2.2              KernSmooth_2.23-20      stringi_1.7.3          
## [73] parallel_4.1.0          Rcpp_1.0.7              vctrs_0.3.8            
## [76] png_0.1-7               leaflet_2.0.4.1         tidyselect_1.1.1       
## [79] xfun_0.24