6.1 Lab Goals

Chapter 4 and Chapter 5 have introduced the fundamentals of spatial data manipulation and map-making, particularly, producing single-variable thematic maps. In this chapter, we will learn to:

  • visualize time series data
  • create cartograms
  • make bivariate choropleth map

6.2 Good Practice

6.2.1 Organizing Folders & Sub-folders

Under the course folder, please create a folder called “lab6”. Next, in the lab6 folder, please create two sub-folders including a “data” folder and a “plot” folder.

6.2.2 Data

This chapter continues to explore spatial data sources other than US Census Bureau. We will use the 311 Service Requests in 2023 and Parks, Medians, and Parkway Trees.

311 data contains pertinent information related to inquiries and service requests for the City and County of Denver. It is aggregated, often publicly available, records of non-emergency service requests, questions, and complaints made by residents to local government through the 311 system.

Please follow the steps below to download data, unzip it and move the data to the required folder.

If there you have any questions about the above-mentioned steps, please refer to Chapter 3.2.3 for detailed instructions.

6.2.3 Launching R Studio

Again, we would like to start a new project from scratch with a clean R Script. Please do the following steps. If you have any questions about these steps, please refer to the relevant chapters for help.

  • Step 1: Make sure all existing R projects are properly closed.
    • If not, please close it by going to File –> Close Project –> Save changes (see Chapter 2.5).
  • Step 2: Create a New Project using Existing Directory, navigate to lab6, click open, then Create Project. (see Chapter 1.3).
  • Step 3: Create a New Script by go to File –> New File –> R Script. Save the script by giving it a proper name.

6.2.4 Before Start

Heads-Up!

All scripts are non-copyable. When writing codes, please try to understand:

  • Code logics
  • Purposes of each line and function
  • the input and output


6.3 Libraries & Input Data

Load libraries.

### Library
library(dplyr)        # data manipulation
library(tidyverse)    # tidy data 
library(reshape2)     # one-hot encoding 

library(sf)           # simple features     
library(geojsonsf)    # read geojson
library(tmap)         # mapping 
library(ggplot2)      # plot 

library(viridis)      # viz - color

library(cartogram)    # map - cartogram
library(pals)         # bivariate color palette
library(corrplot)     # correlation matrix

Input data

# Denver neighborhood (geojson --> simple features)
sf_den_ngbr <- geojson_sf("data/Denver_neighborhood.geojson")

# Denver 311 service requests (csv -- > data frame)
df_den_311 <- read.csv("data/2023_denver_311_service_request_cleaned.csv")

# Denver Trees (around 1 minute to read the data)
sf_den_tree <- geojson_sf("data/denver_trees.geojson")

6.4 Viz: 311 Service Requests by Week

In the previous chapter, we have played with AirBnb data.

This week, we will play with 311 Service Requests in Denver. 311 is a special telephone number supported in many communities in the United States that provides access to non-emergency municipal services like reporting abandoned vehicles, noise complaints, and graffiti.

6.4.1 Exploring 311 data

Before any analysis and visualization tasks, the first step is to learn about the data and check what are the attributes (or columns) in the data.

[Q1] In your opinion, which column in df_den_311 is most interesting and briefly explain why.

# open data
View(df_den_311)

Let’s have a look at the distributions of variables, e.g., the Agency. The table() builds a contingency table of the counts at each combination of factor levels

# contingency table of the attribute 'Agency'
tb_311_agency_freq <- table(df_den_311$Agency)
head(tb_311_agency_freq)
## 
##                              311                     City Council 
##                             8818                              189 
##                 Clerk & Recorder Community Planning & Development 
##                              375                            21586 
##                 Excise & License                  External Agency 
##                             9870                             3965

Create a barplot to show the count of agencies mentioned in the 311 dataset.

# plot: contingency table
barplot(tb_311_agency_freq, main="311 Requests - Agency", las=2, cex.names = 0.5)

6.4.2 Data Type

Check data type of each column. [Q2] What is the data type of the column Case_Created_dttm.

# data type
str(df_den_311)
## 'data.frame':    87321 obs. of  13 variables:
##  $ ObjectId             : chr  "ID000001" "ID000002" "ID000003" "ID000004" ...
##  $ Case_Created_dttm    : chr  "Sun, 31 Dec 2023 14:50:12 GMT" "Sun, 31 Dec 2023 13:25:09 GMT" "Sun, 31 Dec 2023 11:00:02 GMT" "Sun, 31 Dec 2023 13:30:04 GMT" ...
##  $ Case_Closed_dttm     : chr  "1/2/2024 6:55:27 AM" NA "1/8/2024 1:55:27 PM" NA ...
##  $ Incident_Address_1   : chr  "5063 W Radcliff Ave" "2929 S Elm St" "2222 N Broadway" "9132 E 34th Ave" ...
##  $ First_Call_Resolution: chr  "Y" "Y" "Y" "Y" ...
##  $ Case_Summary         : chr  "ROWE - Driveway Block" "Snow Removal Complaint/Sidewalk" "Business Complaint" "Snow Removal Complaint/Sidewalk" ...
##  $ Case_Source          : chr  "PocketGov" "PocketGov" "PocketGov" "PocketGov" ...
##  $ Council_District     : int  2 4 9 8 7 1 0 8 1 8 ...
##  $ Case_Status          : chr  "Closed - Answer Provided" "Routed to Agency" "Closed - Answer Provided" "Routed to Agency" ...
##  $ Longitude            : num  -105 -105 -105 -105 -105 ...
##  $ Latitude             : num  39.6 39.7 39.8 39.8 39.7 ...
##  $ Agency               : chr  "Public Works" "Community Planning & Development" "Excise & License" "Community Planning & Development" ...
##  $ Police_District      : int  4 3 6 5 4 1 0 5 1 5 ...

The column Case_Created_dttm records at what date and time the “311 request” was received. Let’s convert this column to a particular data type as ** date **. Read more about the datetime format in R.

# have a look at the first values in column "Case_Created_dttm"
head(df_den_311$Case_Created_dttm)  
## [1] "Sun, 31 Dec 2023 14:50:12 GMT" "Sun, 31 Dec 2023 13:25:09 GMT"
## [3] "Sun, 31 Dec 2023 11:00:02 GMT" "Sun, 31 Dec 2023 13:30:04 GMT"
## [5] "Sun, 31 Dec 2023 13:25:17 GMT" "Sun, 31 Dec 2023 10:40:02 GMT"
# convert the data type of "Case_Created_dttm" as "date"  
df_den_311$case_create_date <- as.Date(df_den_311$Case_Created_dttm, 
                                       format = "%a, %d %b %Y %H:%M:%S")

# extract week numbers 
df_den_311$case_create_wk <- strftime(df_den_311$case_create_date, 
                                      format = "%V") %>% as.numeric()

[Q3] After transformation, What is the data type of the column case_create_date?

# Check data type again. 
str(df_den_311)
## 'data.frame':    87321 obs. of  15 variables:
##  $ ObjectId             : chr  "ID000001" "ID000002" "ID000003" "ID000004" ...
##  $ Case_Created_dttm    : chr  "Sun, 31 Dec 2023 14:50:12 GMT" "Sun, 31 Dec 2023 13:25:09 GMT" "Sun, 31 Dec 2023 11:00:02 GMT" "Sun, 31 Dec 2023 13:30:04 GMT" ...
##  $ Case_Closed_dttm     : chr  "1/2/2024 6:55:27 AM" NA "1/8/2024 1:55:27 PM" NA ...
##  $ Incident_Address_1   : chr  "5063 W Radcliff Ave" "2929 S Elm St" "2222 N Broadway" "9132 E 34th Ave" ...
##  $ First_Call_Resolution: chr  "Y" "Y" "Y" "Y" ...
##  $ Case_Summary         : chr  "ROWE - Driveway Block" "Snow Removal Complaint/Sidewalk" "Business Complaint" "Snow Removal Complaint/Sidewalk" ...
##  $ Case_Source          : chr  "PocketGov" "PocketGov" "PocketGov" "PocketGov" ...
##  $ Council_District     : int  2 4 9 8 7 1 0 8 1 8 ...
##  $ Case_Status          : chr  "Closed - Answer Provided" "Routed to Agency" "Closed - Answer Provided" "Routed to Agency" ...
##  $ Longitude            : num  -105 -105 -105 -105 -105 ...
##  $ Latitude             : num  39.6 39.7 39.8 39.8 39.7 ...
##  $ Agency               : chr  "Public Works" "Community Planning & Development" "Excise & License" "Community Planning & Development" ...
##  $ Police_District      : int  4 3 6 5 4 1 0 5 1 5 ...
##  $ case_create_date     : Date, format: "2023-12-31" "2023-12-31" ...
##  $ case_create_wk       : num  52 52 52 52 52 52 52 52 52 52 ...

6.4.3 Top 5 Service Reuqests Categories

The column Case_Summary categorizes people’s 311 requests into different types. Let’s check what issues that people have reported and their frequency.

[Q4] What are the top 5 categories of 311 service requests?

# Case summary - contingency table - ordered by frequency 
df_311_summary <- table(df_den_311$Case_Summary) %>% 
  as.data.frame() %>% arrange(-Freq)
# check
View(df_311_summary)

OK, there are over 1000 categories of service requests. Let’s only focus on the top 5 frequently mentioned category and save it as a new data df_den_311_top.

# extract 311 requests in the most mentioned categories 
df_den_311_top <- df_den_311 %>% 
  filter(Case_Summary %in% df_311_summary$Var1[1:5]) %>% 
  select(case_create_wk, Case_Summary) %>% 
  mutate(n=1)

Count the number of 311 service requests in the top 5 categories by week.

#  (one-hot encoding)
df_den_311_top <- dcast(df_den_311_top, 
                        case_create_wk ~ Case_Summary, 
                        fill = 0, value.var = "n",
                        fun.aggregate = length)

head(df_den_311_top)
##   case_create_wk Application Encampment Reporting Sheriff - Abandoned vehicle
## 1              1          26                  180                          63
## 2              2          43                  305                          57
## 3              3          43                  187                          40
## 4              4          30                  170                          47
## 5              5          40                  232                          60
## 6              6          46                  279                          47
##   Snow Removal Complaint/Sidewalk Weeds/Vegetation Violation
## 1                            1732                          1
## 2                            1169                          0
## 3                            1221                          5
## 4                             847                          1
## 5                             516                          3
## 6                             177                          9

From here, we are going to change the structure of the data df_den_311_top and reshape (or tidy) it in a specific way so that:

  • Each variable must have its own column.
  • Each observation must have its own row.
  • Each value must have its own cell.

The process if tidying is illustrated in Figure 6.1. Figure 6.1: Diagram of data tidying. 

# Tidy data using the function pivot_longer()
df_den_311_top_tidy <- df_den_311_top %>% 
  pivot_longer(cols = 2:ncol(.), names_to = "category", values_to ="count")

[Q5] Please compare the two datasets before and after tidying, e.g., df_den_311_top and df_den_311_top_tidy. Try to use your own words to explain how the “tidy” process changes the data structure.

# BEFORE tidying:  "df_den_311_top_tidy"
View(df_den_311_top)

# AFTER tidying:  "df_den_311_top_tidy"
View(df_den_311_top_tidy)

6.4.4 Time Series Plot

Let’s make a basic plot exploring how the top 5 categories of 311 requests change over times (by week) by using the tidied data df_den_311_top_tidy.

# viz: time series 
p1_311 <- ggplot(df_den_311_top_tidy, 
                 aes(x=case_create_wk, y=count, fill=category)) +
  geom_area(alpha=0.6 , linewidth=.5, colour="white")

# print plot
p1_311

Improve the plot by using the viridis color palette from the viridis library.

# viz: improve 
p2_311 <- ggplot(df_den_311_top_tidy, 
                 aes(x=case_create_wk, y=count, fill=category)) +
  geom_area(alpha=0.6 , linewidth=.5, colour="white") + 
  scale_fill_viridis(discrete = T) + 
  theme_bw() +
  xlab("Weeks in 2023") + 
  ylab(" Count of 311 Service Requests")
  
# print plot
p2_311

Export

# export as .pdf 
pdf("plot/plot2_311_service_cat_week.pdf", width = 10, height = 6)
print(p2_311)
dev.off()


# export as .png
png("plot/plot2_311_service_cat_week.png")
print(p2_311)
dev.off()

Okay, while the plot looks better now, the week number is very confusing. Let’s add some labels by connecting week numbers to seasons as suggested by Figure 6.2.

[Q6] Please briefly discuss your findings from this analysis (plot3_311_service_cat_week_reference.pdf).

Figure 6.2: When four seasons start in 2024.
Figure 6.2: When four seasons start in 2024.


Add reference lines for four seasons.

# create a df showing when the four seasons start 
df_311_season_label <- data.frame(
  x_season = c(5,18,31,44),
  y_season = 2000, 
  label = c("Winter","Spring", "Summer", "Fall")
)

# viz: improve by adding annotation & reference lines  
p3_311 <- ggplot(df_den_311_top_tidy, aes(x=case_create_wk, y=count, fill=category)) +
  geom_area(alpha=0.6 , linewidth=.5, colour="white") + 
  scale_fill_viridis(discrete = T) + 
  theme_bw() +
  geom_vline(xintercept = c(10,23,36,49),colour="black", linetype = "dotted") + 
  geom_label(data=df_311_season_label, aes(x = x_season, y = y_season, label = label), fill = "white", linewidth=0.5)  + 
  xlab("Weeks in 2023") + 
  ylab(" Count of 311 Service Requests") 

p3_311

# export 
pdf("plot/plot3_311_service_cat_week_reference.pdf", width = 10, height = 6)
print(p3_311)
dev.off()

6.5 Viz: Mapping 311 Requests

This section moves from time-series to spatial visualization of the 311 data.

[Q7] Please open the data frame df_den_311 and report which two columns can be used to convert “df_den_311” to spatial data / simple features?

# open
View(df_den_311)

# list all column names 
names(df_den_311)

6.5.1 From Dataframe to Simple Features

Create simple features of 311 service requests and save the data as sf_den_311.

# Convert to spatial points 
sf_den_311 <- st_as_sf(df_den_311, coords = c("Longitude", "Latitude"), crs = 4326)

Let’s check the CRS (coordinate reference system) of all simple features we have and project them to the 2D CRS, NAD83 / Colorado North (ftUS), which has its centroid in Colorado.

The value 2231 is the unique index for the CRS NAD83 / Colorado North (ftUS). More CRS could be found here https://epsg.io/.

# four spatial data (simple features). Let's check their coordination systems 
sf_den_311$geometry
sf_den_tree$geometry
sf_den_ngbr$geometry

Project three simple features, i.e., sf_den_311, sf_den_tree, sf_den_ngbr from WGS84 to NAD83 / Colorado North (ftUS) .

# project the data to NAD83 / Colorado North (ftUS)
# https://epsg.io/2231 - NAD83 / Colorado North (ftUS)
sf_den_311 <- st_transform(sf_den_311, 2231)
sf_den_tree <- st_transform(sf_den_tree, 2231)
sf_den_ngbr <- st_transform(sf_den_ngbr, 2231)

Quick plot of three simple features sf_den_311, sf_den_tree, sf_den_ngbr

## map "sf_den_311" as symbols
tmap_mode("plot")
tm_shape(sf_den_311) + tm_symbols(size = 0.2, fill_alpha = 0.1)

## map the first 1000 records in "sf_den_tree" as dots
tmap_mode("view")
tm_shape(sf_den_tree[1:5000, ]) + tm_dots(fill = "DIAMETER", fill_alpha = 0.2)
## map "sf_den_ngbr" as polygons 
tm_shape(sf_den_ngbr) + tm_polygons(fill_alpha = 0.5)

6.5.2 Mapping 311 Calls by Neighborhood

As shown in the plot, the sf_den_311 are points. They overlap together and make patterns hard to distinguish. Thus, we will agregate 311 calls to neighborhood level.

The following chunk of code counts the number of 311 service requests in each neighborhood based on the spatial intersects relationship st_intersects() .

### Spatial Data Processing 
# Count 311 (2023) in each neighborhood  
sf_den_ngbr_311 <- sf_den_ngbr
sf_den_ngbr_311["n_311"] <- lengths(st_intersects(sf_den_ngbr_311, sf_den_311))

Quick plot of 311 calls by neighborhood.

[Q8] Have a look at the simple feature data frame sf_den_ngbr_311 and and report which neighborhood has received the highest number of complaints in 2023?

view(sf_den_ngbr_311)

Quick plot of count of 311 in neighborhood

# set mode to static mapping
tmap_mode("plot")

map1_ngbr_311 <- tm_shape(sf_den_ngbr_311) + 
  tm_polygons(fill = "n_311", 
              fill.scale = tm_scale_intervals(values = "brewer.yl_gn_bu"),
              fill.legend = tm_legend("311 Calls in 2023", position = c("bottom", "right"))) + 
  tm_compass(position = c("top", "left")) +
  tm_scalebar(width = 7, position = c("top", "left"))

map1_ngbr_311

# export 
pdf("plot/map1_ngbr_311_requests_in_2023.pdf", width = 10, height = 8)
print(map1_ngbr_311)
dev.off()

6.5.3 Cartograms

This and the following sections use the same dataset sf_den_ngbr_311 and variable n_311 to explore additional mapping techniques, i.e., cartogram.

A cartogram is a thematic map of a set of features, in which their geographic size is altered to be directly proportional to a selected variable.

# cartogram 1
# Construct a continuous area cartogram by a rubber sheet distortion algorithm
sf_den_ngbr_311_carto <- cartogram_cont(sf_den_ngbr_311, "n_311", itermax = 10)

map2_ngbr_311_c1 <- tm_shape(sf_den_ngbr_311_carto) + 
  tm_polygons(fill = "n_311", 
              fill.scale = tm_scale_intervals(values = "brewer.yl_gn_bu"),
              fill.legend = tm_legend("311 Calls in 2023", position = c("bottom", "right"))) + 
  tm_compass(position = c("top", "left")) +
  tm_scalebar(width = 7, position = c("top", "left"))

map2_ngbr_311_c1

pdf("plot/map2_ngbr_311_requests_in_2023_carto1.pdf", width = 10, height = 8)
print(map2_ngbr_311_c1)
dev.off()

6.5.4 Cartogram - Dorling

Construct a cartogram which represents each geographic region as non-overlapping circles (Dorling 1996). please try to play with a different k value, see what is happening.

# cartogram 2
# Construct a cartogram which represents each geographic region as non-overlapping circles (Dorling 1996).
sf_den_ngbr_311_carto_dorling <- cartogram_dorling(sf_den_ngbr_311, weight = "n_311", itermax = 5, k=0.7)

map3_ngbr_311_c2 <- tm_shape(sf_den_ngbr) + tm_polygons(fill_alpha = 0, col="grey") +
  tm_shape(sf_den_ngbr_311_carto_dorling) + 
  tm_polygons(fill = "n_311", 
              fill.scale = tm_scale_intervals(values = "brewer.yl_gn_bu"),
              fill.legend = tm_legend("311 Calls in 2023", position = c("bottom", "right"))) + 
  tm_compass(position = c("top", "left")) +
  tm_scalebar(width = 7, position = c("top", "left"))

map3_ngbr_311_c2

pdf("plot/map3_ngbr_311_requests_in_2023_carto2_dorling.pdf", width = 10, height = 8)
print(map3_ngbr_311_c2)
dev.off()

[Q9] Please compare the three maps: map1_ngbr_311_c1, map2_ngbr_311_c1, map3_ngbr_311_c1 and briefly discuss which map do you like the most and why?

6.6 Viz: Tree Maps

This section uses the Parks, Medians, and Parkway Trees in Denver to produce a bivariate choropleth map.

Again, have a look at the sf_den_tree data and explore what columns are in it.

[Q10] Which column in sf_den_tree do you like the most and why?

View(sf_den_tree)

6.6.1 Spatial Data Manipulation

Join the two datasets sf_den_tree and sf_den_ngbr by using st_join() , so that we can know in which neighborhood each tree record is located.

## spatial join to get neighobrhood information 
sf_den_tree_ngbr <- st_join(sf_den_tree, sf_den_ngbr["neighbourhood"], 
                            left = T, join = st_intersects)
## have a look at the data after joining
View(sf_den_tree_ngbr)

The column SPECIES_COMMON tells about the species of trees. Let’s check what are & how many trees of different species.

# count tree species 
df_tree_spec <- table(sf_den_tree_ngbr$SPECIES_COMMON) %>% 
  as.data.frame() %>% 
  arrange(-Freq)
# open data
View(df_tree_spec)
## the top 10 tree species are:
df_tree_spec$Var1[1:10]

Count trees in the top 10 species in each neighborhood using One-hot encoding.

# neighborhood tree species count 
df_ngbr_tree <- sf_den_tree_ngbr %>% 
  filter(SPECIES_COMMON %in% df_tree_spec$Var1[1:10]) %>% 
  select(neighbourhood, SPECIES_COMMON) %>% 
  mutate(n=1)

# remove geometry 
st_geometry(df_ngbr_tree) <- NULL

# one-hot encoding 
df_ngbr_tree <- dcast(df_ngbr_tree, neighbourhood ~ SPECIES_COMMON, 
                      fill = 0, value.var = "n", fun.aggregate = length)

# remove NA 
df_ngbr_tree <- na.omit(df_ngbr_tree)
# open the final data 
View(df_ngbr_tree)

6.6.2 Correlation Plot

Let’s check if the count of tree species have any correlation among themselves

## calculate correlation matrix
m <- cor(df_ngbr_tree[, 2:ncol(df_ngbr_tree)])
m
##                      Ash, Green  Ash, White Crabapple, Flowering Elm, Siberian
## Ash, Green            1.0000000  0.85493608            0.7000782   -0.10504628
## Ash, White            0.8549361  1.00000000            0.7881865   -0.07973038
## Crabapple, Flowering  0.7000782  0.78818645            1.0000000    0.10923967
## Elm, Siberian        -0.1050463 -0.07973038            0.1092397    1.00000000
## Hackberry, Common     0.8240476  0.80128051            0.6956260   -0.05773733
## Honey Locust          0.8359359  0.85990268            0.7233975   -0.16239784
## Linden, Littleleaf    0.7375544  0.79095863            0.6589631   -0.19874075
## Maple, Silver         0.3779506  0.53767817            0.4717454    0.07677355
## Pear, Flowering       0.7985263  0.85038143            0.7361297   -0.12387270
## Pine, Austrian        0.5342673  0.43598443            0.5738679   -0.01126116
##                      Hackberry, Common Honey Locust Linden, Littleleaf
## Ash, Green                  0.82404762    0.8359359          0.7375544
## Ash, White                  0.80128051    0.8599027          0.7909586
## Crabapple, Flowering        0.69562596    0.7233975          0.6589631
## Elm, Siberian              -0.05773733   -0.1623978         -0.1987408
## Hackberry, Common           1.00000000    0.8953511          0.7925347
## Honey Locust                0.89535108    1.0000000          0.7817069
## Linden, Littleleaf          0.79253469    0.7817069          1.0000000
## Maple, Silver               0.25071448    0.4208601          0.1565366
## Pear, Flowering             0.84944138    0.9089777          0.7765225
## Pine, Austrian              0.54097165    0.4406928          0.5747946
##                      Maple, Silver Pear, Flowering Pine, Austrian
## Ash, Green              0.37795057       0.7985263     0.53426728
## Ash, White              0.53767817       0.8503814     0.43598443
## Crabapple, Flowering    0.47174539       0.7361297     0.57386785
## Elm, Siberian           0.07677355      -0.1238727    -0.01126116
## Hackberry, Common       0.25071448       0.8494414     0.54097165
## Honey Locust            0.42086008       0.9089777     0.44069280
## Linden, Littleleaf      0.15653660       0.7765225     0.57479460
## Maple, Silver           1.00000000       0.3732849    -0.06196431
## Pear, Flowering         0.37328487       1.0000000     0.40626784
## Pine, Austrian         -0.06196431       0.4062678     1.00000000

Create a correlation plot.

[Q11] Please use your own words to explain your findings from the correlation matrix m and the plot. E.g., which two variables have the highest or lowest correlation coefficients?

# correlation plot. 
corrplot(m, method = "ellipse", type = "upper")

Join the data df_ngbr_tree to the spatial boundaries of neighborhood sf_den_ngbr.

# left join to simple features 
sf_den_ngbr_tree <- left_join(sf_den_ngbr, df_ngbr_tree, by = "neighbourhood")
## plot all columns together for explorations 
plot(sf_den_ngbr_tree, max.plot = 400)

6.6.3 Bivariate Mapping

Let’s say we are interested in the two species “Honey Locust” and “Hackberry, Common”, and display the counts of two species on a single map.

# two tree species 
lt_tree_specs_lbl <- c("Honey Locust", "Hackberry, Common")

At the final map, for each species, we want to divide the values into three categories using the “quantile” classification methods. You could try other classification methods that we have seen in Chapter 5.6.6.

## Bivariate mapping 
tm_shape(sf_den_ngbr_tree) + 
  tm_polygons(fill = tm_vars(lt_tree_specs_lbl, multivariate = TRUE),
              fill_alpha = 0.8, 
              fill.scale = tm_scale_bivariate(
                scale1 = tm_scale_intervals(style = "quantile", n=3, labels =  c("L", "M", "H")),
                scale2 = tm_scale_intervals(style = "quantile", n=3, labels =  c("L", "M", "H")),
                values = "stevens.pinkblue"),
              tm_legend(frame = F, position =c("left", "top")))

Improve the map by updating legends, adding scale bars and north arrow.

[Q12] Please try to change the parameters of bivariate maps and paste a screenshot of your final project.

  • Try other classification styles: sd, equal, pretty, quantile, jenks, kmeans
  • Change number of classes and their labels
  • Different bivariate color palettes (see screenshots).
Figure 6.3: Bivariate Color Palette
Figure 6.3: Bivariate Color Palette
map4_ngbr_trees_biv <- tm_shape(sf_den_ngbr_tree) + 
  tm_polygons(fill = tm_vars(lt_tree_specs_lbl, multivariate = TRUE),
              fill_alpha = 0.8, 
              fill.scale = tm_scale_bivariate(
                scale1 = tm_scale_intervals(style = "kmeans", n=3, labels =  c("L", "M", "H")),
                scale2 = tm_scale_intervals(style = "kmeans", n=3, labels =  c("L", "M", "H")),
                values = "stevens.pinkblue"),
              tm_legend(frame = F, position =c("left", "top"))) + 
  tm_compass(position = c("right", "bottom")) +
  tm_scalebar(width = 7, position =c("right", "bottom")) +
  tm_layout(frame = F)


map4_ngbr_trees_biv

Export.

# export 
pdf("plot/map4_neighborhood_bivariate_trees.pdf", width = 10, height = 7)
print(map4_ngbr_trees_biv)
dev.off()

6.7 Close & Exit

Congratulations!! You have completed the entire tutorial and learnt new spatial analysis and visualization techniques!! Excellent work.

Please go “File”–> “Close Project” – a pop window asking “Do you want to save these changes” –> “Yes”.

Don’t forget to submit the lab6 report and your script to Canvas.