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('tidyverse') 
#install.packages('sf)

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

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

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 need to have a file with major cities in Colombia. You may download it from https://simplemaps.com/data/co-cities. Then, you need to uncompress it and move it under your working directory.

Let’s start by listing the available files:

list.files(pattern=c('csv'))

## [1] "Censo_Nacional_Agropecuario_Uso.csv"           
## [2] "edidiv.csv"                                    
## [3] "EVA_Boyaca_2019.csv"                           
## [4] "EVA_Boyaca.csv"                                
## [5] "EVA.csv"                                       
## [6] "Evaluaciones_Agropecuarias_Municipales_EVA.csv"
## [7] "NBI.csv"                                       
## [8] "stder_frutales_2020.csv"                       
## [9] "stder_oleag_2020.csv"

list.files(pattern=c('txt'))

## [1] "cities.txt"          "colombia_cities.txt" "test.txt"

#change the file path as needed
list.files('./shapefiles/68_STDER')

## [1] "municipios.geojson" "municipios.qmd"

Now, let’s proceed to read the files and show the output object:

(frutales = read_csv("stder_frutales_2020.csv",show_col_types = FALSE))

To read the municipalities we are interested:

(municipios =  st_read('./shapefiles/68_STDER/municipios.geojson'))

## Reading layer `municipios' from data source 
##   `/Users/ivanlizarazo/Documents/ivan/GB/datos/shapefiles/68_STDER/municipios.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 87 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
## Geodetic CRS:  WGS 84

To read the cities:

(cities = read_delim("colombia_cities.txt",col_names = FALSE, delim = ';', show_col_types = FALSE))

Note that there are 100 rows representing major 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 “X4” (latitude) and “X5” (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("X5", "X4"))

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 department we need to create an spatial join:

# find points within polygons
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 municipios .

Now, let’s select the rows whose DPTO_CNMBR are filled out:

stder.cities = dplyr::filter(sf.cities.joined, DPTO_CCDGO=='68')

stder.cities 

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. However, their data types are different.

class(frutales$Cod_Mun)

## [1] "numeric"

class(municipios$COD_MUN)

## [1] "character"

Thus, we need to fix this:

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

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

munic_frutales

This code customizes the colors for the polygons. 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')

## Warning in classIntervals(munic_frutales$max_prod, n = 6, style = "fisher"): var
## has missing values, omitted in finding classes

#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

(munic_frutales$mid <- sf::st_centroid(munic_frutales$geometry))

## Geometry set for 87 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -74.17032 ymin: 5.778074 xmax: -72.56969 ymax: 7.57697
## Geodetic CRS:  WGS 84
## First 5 geometries:

## POINT (-73.53038 6.182207)

## POINT (-73.82846 5.778074)

## POINT (-73.01163 6.716764)

## POINT (-73.62345 5.955617)

## POINT (-73.21518 6.647059)

LONG = st_coordinates(munic_frutales$mid)[,1]

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:

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

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

munic_oleag = left_join(municipios, oleag, by = c("COD_MUN" = "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.5, fontfamily = "sans") +
  tm_shape(stder.cities) + tm_symbols(shape = 2, col = "red", size = 0.25) +
  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. Reproducibility

If you use code from here please cite this work as Lizarazo, I., 2022. Getting started with thematic maps. Available at https://rpubs.com/ials2un/thematic_maps_v1.

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## 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        sp_1.4-5                highr_0.9              
## [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            generics_0.1.0          farver_2.1.0           
## [31] ellipsis_0.3.2          withr_2.4.2             leafsync_0.1.0         
## [34] mime_0.11               magrittr_2.0.1          crayon_1.4.1           
## [37] evaluate_0.14           fansi_0.5.0             lwgeom_0.2-8           
## [40] class_7.3-19            tools_4.1.0             hms_1.1.1              
## [43] lifecycle_1.0.0         stringr_1.4.0           munsell_0.5.0          
## [46] compiler_4.1.0          jquerylib_0.1.4         e1071_1.7-8            
## [49] rlang_0.4.11            units_0.7-2             grid_4.1.0             
## [52] tmaptools_3.1-1         dichromat_2.0-0         htmlwidgets_1.5.3      
## [55] crosstalk_1.1.1         leafem_0.1.6            base64enc_0.1-3        
## [58] rmarkdown_2.9           wk_0.5.0                gtable_0.3.0           
## [61] codetools_0.2-18        abind_1.4-5             DBI_1.1.1              
## [64] markdown_1.1            R6_2.5.0                knitr_1.33             
## [67] bit_4.0.4               utf8_1.2.2              KernSmooth_2.23-20     
## [70] stringi_1.7.3           parallel_4.1.0          Rcpp_1.0.7             
## [73] vctrs_0.3.8             png_0.1-7               leaflet_2.0.4.1        
## [76] tidyselect_1.1.1        xfun_0.24