02_spatial-data.qmd

::: {.content-hidden unless-format="html"}
$$
\newcommand{\tr}{\mathrm{tr}}
\newcommand{\rank}{\mathrm{rank}}
\newcommand{\plim}{\operatornamewithlimits{plim}}
\newcommand{\diag}{\mathrm{diag}}
\newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}}
\newcommand{\Var}{\mathrm{Var}}
\newcommand{\Exp}{\mathrm{E}}
\newcommand{\Cov}{\mathrm{Cov}}
\newcommand\given[1][]{\:#1\vert\:}
\newcommand{\irow}[1]{%
\begin{pmatrix}#1\end{pmatrix}
}
$$
:::

# Data Manipulation & Visualization

### Required packages {.unnumbered}

```{r, message = FALSE, warning = FALSE, results = 'hide'}
pkgs <- c("sf", "gstat", "mapview", "nngeo", "rnaturalearth", "dplyr", "areal",
          "nomisr", "osmdata", "OpenStreetMap", "tidyr", "texreg", "downlit", "xml2") 
lapply(pkgs, require, character.only = TRUE)

```

For mapping

```{r, message = FALSE, warning = FALSE, results = 'hide'}
pkgs <- c("tmap", "tmaptools", "viridis", "viridisLite", 
          "ggplot2", "ggthemes", "rmapshaper", "cowplot", "gridExtra") 
lapply(pkgs, require, character.only = TRUE)

```

### Session info {.unnumbered}

```{r}
sessionInfo()

```

### Reload data from pervious session {.unnumbered}

```{r}
load("_data/msoa_spatial.RData")
load("_data/ulez_spatial.RData")
load("_data/pollution_spatial.RData")
load("_data/pubs_spatial.RData")
```

## Manipulation and linkage

Having data with geo-spatial information allows to perform a variety of methods to manipulate and link different data sources. Commonly used methods include 1) subsetting, 2) point-in-polygon operations, 3) distance measures, 4) intersections or buffer methods.

The [online Vignettes of the sf package](https://r-spatial.github.io/sf/articles/) provide a comprehensive overview of the multiple ways of spatial manipulations.

#### Check if data is on common projection

```{r}
st_crs(msoa.spdf) == st_crs(pol.spdf)
st_crs(msoa.spdf) == st_crs(pubs.spdf)
st_crs(msoa.spdf) == st_crs(ulez.spdf)
```

The spatial data files are on different projections. Before we can do any spatial operations with them, we have to transform them into a common projection.

```{r}
# MSOA in different crs --> transform
pol.spdf <- st_transform(pol.spdf, crs = st_crs(msoa.spdf))
pubs.spdf <- st_transform(pubs.spdf, crs = st_crs(msoa.spdf))
ulez.spdf <- st_transform(ulez.spdf, crs = st_crs(msoa.spdf))


# Check if all geometries are valid, and make valid if needed
msoa.spdf <- st_make_valid(msoa.spdf)

```

The `st_make_valid()` function can help if the spatial geometries have some problems such as holes or points that don't match exactly.

### Subsetting

We can subset spatial data in a similar way as we subset conventional data.frames or matrices. For instance, below we simply reduce the pollution grid across the UK to observations in London only.

```{r}

# Subset to pollution estimates in London
pol_sub.spdf <- pol.spdf[msoa.spdf, ] # or:
pol_sub.spdf <- st_filter(pol.spdf, msoa.spdf)
mapview(pol_sub.spdf)

```

Or we can reverse the above and exclude all intersecting units by specifying `st_disjoint` as alternative spatial operation using the `op =` option (note the empty space for column selection). `st_filter()` with the `.predicate` option does the same job. See the [sf Vignette](https://cran.r-project.org/web/packages/sf/vignettes/sf3.html) for more operations.

```{r}
# Subset pubs to pubs not in the ulez area
sub2.spdf <- pubs.spdf[ulez.spdf, , op = st_disjoint] # or:
sub2.spdf <- st_filter(pubs.spdf, ulez.spdf, .predicate = st_disjoint)
mapview(sub2.spdf)
```

We can easily create indicators of whether an MSOA is within ulez or not.

```{r}
msoa.spdf$ulez <- 0

# intersecting lsoas
within <- msoa.spdf[ulez.spdf,]

# use their ids to create binary indicator 
msoa.spdf$ulez[which(msoa.spdf$MSOA11CD %in% within$MSOA11CD)] <- 1
table(msoa.spdf$ulez)
```

### Point in polygon

We are interested in the number of pubs in each MSOA. So, we count the number of points in each polygon.

```{r}
# Assign MSOA to each point
pubs_msoa.join <- st_join(pubs.spdf, msoa.spdf, join = st_within)

# Count N by MSOA code (drop geometry to speed up)
pubs_msoa.join <- dplyr::count(st_drop_geometry(pubs_msoa.join),
                               MSOA11CD = pubs_msoa.join$MSOA11CD,
                               name = "pubs_count")
sum(pubs_msoa.join$pubs_count)

# Merge and replace NAs with zero (no matches, no pubs)
msoa.spdf <- merge(msoa.spdf, pubs_msoa.join,
                   by = "MSOA11CD", all.x = TRUE)
msoa.spdf$pubs_count[is.na(msoa.spdf$pubs_count)] <- 0

```


### Distance measures

We might be interested in the distance to the nearest pub. Here, we use the package `nngeo` to find k nearest neighbours with the respective distance.

```{r}
# Use geometric centroid of each MSOA
cent.sp <- st_centroid(msoa.spdf[, "MSOA11CD"])

# Get K nearest neighbour with distance
knb.dist <- st_nn(cent.sp, 
                  pubs.spdf,
                  k = 1,             # number of nearest neighbours
                  returnDist = TRUE, # we also want the distance
                  progress = FALSE)
msoa.spdf$dist_pubs <- unlist(knb.dist$dist)
summary(msoa.spdf$dist_pubs)

```

### Intersections + Buffers

We may also want the average pollution within 1 km radius around each MSOA centroid. Note that it is usually better to use a ego-centric method where you calculate the average within a distance rather than using the characteristic of the intersecting cells only [@Lee.2008; @Mohai.2007].

Therefore, we first create a buffer with `st_buffer()` around each midpoint and subsequently use `st_intersetion()` to calculate the overlap.

```{r}
# Create buffer (1km radius)
cent.buf <- st_buffer(cent.sp, 
                      dist = 1000) # dist in meters
mapview(cent.buf)


### New version (using areal package)
# We use area-weighted interpolation from the areal package
int.spdf <- aw_interpolate(
  cent.buf,                         # interpolate to
  tid = "MSOA11CD",                 # id of target
  source = pol.spdf,                # the source to be interpolated from
  sid = "ukgridcode",               # source id
  weight = "sum",                   # function for interpolation
  output = "sf",                    # output object
  intensive = "no22011"             # variables to be interpolated
)

# # ### Old version ("by hand")
# # Add area of each buffer (in this constant)
# cent.buf$area <- as.numeric(st_area(cent.buf))
# 
# # Calculate intersection of pollution grid and buffer
# int.df <- st_intersection(cent.buf, pol.spdf)
# int.df$int_area <- as.numeric(st_area(int.df)) # area of intersection
# 
# # And we use the percent overlap areas as the weights to calculate a weighted mean.
# 
# # Area of intersection as share of buffer
# int.df$area_per <- int.df$int_area / int.df$area
# 
# # Aggregate as weighted mean
# int.df <- st_drop_geometry(int.df)
# int.df$no2_weighted <- int.df$no22011 * int.df$area_per
# int.df <- aggregate(list(no2 = int.df[, "no2_weighted"]), 
#                     by = list(MSOA11CD = int.df$MSOA11CD),
#                     sum)
```

And we merge that back to the LSOA data.

```{r}
# Back to non-spatial df
int.df <- st_drop_geometry(int.spdf)
names(int.df)[2] <- "no2"

# Merge back to spatial data.frame
msoa.spdf <- merge(msoa.spdf, int.df, by = "MSOA11CD", all.x = TRUE)

mapview(msoa.spdf, zcol = "no2")
```

Note: for buffer related methods, it often makes sense to use population weighted centroids instead of geographic centroids (see [here](https://geoportal.statistics.gov.uk/datasets/ons::middle-layer-super-output-areas-december-2011-ew-population-weighted-centroids/about) for MSOA population weighted centroids). However, often this information is not available.

### and more

There are more spatial operation possible using sf. Have a look at the [sf Cheatsheet](fig/sf.pdf).

![](fig/sf_1.png)


### Air pollution and ethnic minorities

With a few lines of code, we have compiled an original dataset containing demographic information, air pollution, and some infrastructural information.

Let's see what we can do with it.

```{r}
# Define ethnic group shares
msoa.spdf$per_mixed <- msoa.spdf$KS201EW_200 / msoa.spdf$KS201EW0001 * 100
msoa.spdf$per_asian <- msoa.spdf$KS201EW_300 / msoa.spdf$KS201EW0001 * 100
msoa.spdf$per_black <- msoa.spdf$KS201EW_400 / msoa.spdf$KS201EW0001 * 100
msoa.spdf$per_other <- msoa.spdf$KS201EW_500 / msoa.spdf$KS201EW0001 * 100

# Define tenure
msoa.spdf$per_owner <- msoa.spdf$KS402EW_100 / msoa.spdf$KS402EW0001 * 100
msoa.spdf$per_social <- msoa.spdf$KS402EW_200 / msoa.spdf$KS402EW0001 * 100

# Non British passport
msoa.spdf$per_nonUK <- (msoa.spdf$KS205EW0001 - msoa.spdf$KS205EW0003)/ msoa.spdf$KS205EW0001 * 100
msoa.spdf$per_nonEU <- (msoa.spdf$KS205EW0001 - msoa.spdf$KS205EW0003 -
                          msoa.spdf$KS205EW0004 - msoa.spdf$KS205EW0005  - 
                          msoa.spdf$KS205EW0006)/ msoa.spdf$KS205EW0001 * 100
msoa.spdf$per_nonUK_EU <- (msoa.spdf$KS205EW0005  + msoa.spdf$KS205EW0006)/ msoa.spdf$KS205EW0001 * 100


# Run regression
mod1.lm <- lm(no2 ~ per_mixed + per_asian + per_black + per_other +
                per_owner + per_social + pubs_count + POPDEN + ulez,
              data = msoa.spdf)

# summary
screenreg(list(mod1.lm), digits = 3)

```

For some examples later, we also add data on house prices. We use the median house prices in 2017 from the [London Datastore](https://data.london.gov.uk/dataset/average-house-prices).

```{r house-prices, cache=TRUE}
# Download
hp.link <- "https://data.london.gov.uk/download/average-house-prices/bdf8eee7-41e1-4d24-90ce-93fe5cf040ae/land-registry-house-prices-MSOA.csv"
hp.df <- read.csv(hp.link)
hp.df <- hp.df[which(hp.df$Measure == "Median" &
                       grepl("2011", hp.df$Year)), ]
table(hp.df$Year)

# Aggregate across 2011 values
hp.df$med_house_price <- as.numeric(hp.df$Value)
hp.df <- aggregate(hp.df[, "med_house_price", drop = FALSE],
                   by = list(MSOA11CD = hp.df$Code),
                   FUN = function(x) mean(x, na.rm = TRUE))

# Merge spdf and housing prices
msoa.spdf <- merge(msoa.spdf, hp.df,
                   by = "MSOA11CD",
                   all.x = TRUE, all.y = FALSE)
hist(log(msoa.spdf$med_house_price))
```

### Save spatial data

```{r}
# Save
save(msoa.spdf, file = "_data/msoa2_spatial.RData")
```




## Visualization

A large advantage of spatial data is that different data sources can be connected and combined. Another nice advantage is: you can create very nice maps. And it's quite easy to do! [Stefan Jünger](https://stefanjuenger.github.io/) & [Anne-Kathrin Stroppe](https://www.gesis.org/institut/mitarbeitendenverzeichnis/person/Anne-Kathrin.Stroppe) provide more comprehensive materials on mapping in their [GESIS workshop on geospatial techniques in R](https://github.com/StefanJuenger/gesis-workshop-geospatial-techniques-R-2023).

Many packages and functions can be used to plot maps of spatial data. For instance, ggplot as a function to plot spatial data using `geom_sf()`. I am personally a fan of `tmap`, which makes many steps easier (but sometimes is less flexible).

A great tool for choosing coulour is for instance [Colorbrewer](https://colorbrewer2.org/). `viridis` provides another great resource to chose colours.



### ggplot

Many of you will be familiar with ggplot. Obviously, you can also create very nice maps with ggplot and the `geom_sf()` layer. Particularly for continuous scales which you don't want to cut into bins, ggplot is probably the best choice. The options `coord_sf(datum = NA)` and `theme_map()` help to improve the layout (e.g. omitting the scales and the ticks).


```{r}
gp <- ggplot(msoa.spdf)+
    geom_sf(aes(fill = no2))+
    scale_fill_viridis_c(option = "B")+
    coord_sf(datum = NA)+
    theme_map()+
    theme(legend.position = c(.9, .6))
gp
```

As always with ggplot, we can add different layers of data on top of each other. Below, we combine different layers of `geom_sf()`.

```{r}
# Get some larger scale boundaries
borough.spdf <- st_read(
  dsn = paste0("_data", "/statistical-gis-boundaries-london/ESRI"),
  layer = "London_Borough_Excluding_MHW" # Note: no file ending
)

# transform to only inner lines
borough_inner <- ms_innerlines(borough.spdf)

# Plot with inner lines
gp <- ggplot(msoa.spdf) +
  geom_sf(aes(fill = no2), color = NA) +
  scale_fill_viridis_c(option = "A") +
  geom_sf(data = borough_inner, color = "gray92") +
  geom_sf(
    data = ulez.spdf,
    color = "red",
    fill = NA,
    linewidth = 1.5
  ) +
  coord_sf(datum = NA) +
  theme_map() +
  labs(fill = "NO2 in mg/m3", title = "Air pollution") +
  theme(legend.position = c(.7, .3),
        plot.title = element_text(hjust = 0.5, size = 18),
        legend.background = element_rect(fill = "white", color = NA))
gp
```


And we can combine single maps, e.g. with the gridExtra package and the `grid.arrange()` function. We use exactly the same function, and only change the colouring option in `scale_fill_viridis_c()` which is the viridis colour option for continuous scales.

```{r}
# Plot with inner lines
gp2 <- ggplot(msoa.spdf) +
  geom_sf(aes(fill = med_house_price/1000), color = NA) +
  scale_fill_viridis_c(option = "D") +
  geom_sf(data = borough_inner, color = "gray92") +
  geom_sf(
    data = ulez.spdf,
    color = "red",
    fill = NA,
    linewidth = 1.5
  ) +
  coord_sf(datum = NA) +
  theme_map() +
  labs(fill = "Median/1000", title = "House price") +
  theme(legend.position = c(.7, .3),
        plot.title = element_text(hjust = 0.5, size = 18),
        legend.background = element_rect(fill = "white", color = NA))

gpc <- grid.arrange(gp, gp2, ncol = 2, nrow = 1)

plot(gpc)
```

And we can use the standard functions to export a png is good resolution.


```{r}
png(file =  "fig/London_combined.png", 
    width = 13, height = 6, units = "in",
    res = 300, bg = "white")
plot(gpc)
dev.off()
```


### Tmaps

Another nice package for data visualisation is tmp. It's particularly nice, if you want to cut the colour into different brackets / bins. The `tm_scale_intervals()` function provides very helpful algorithms to cut the variable range into different bins for colouring.

For instance, lets plot the NO2 estimates using tmap + `tm_polygon()` and the fisher algorithm for the colour intervals.

```{r}

# Define colours
cols <- viridis(n = 7, direction = 1, option = "C")

mp1 <- tm_shape(msoa.spdf) +
  tm_polygons(
    fill = "no2",                      # variable for fill colouring
    fill_alpha = 1,                    # transparency 
    fill.scale = tm_scale_intervals(
      style = "fisher",                # algorithm to def cut points
      n = 7,                           # Number of requested cut points
      values = cols                    # colouring
    ),
    fill.legend = tm_legend(
      title = "NO2",
      hist = FALSE
    )
  ) +
  tm_borders(
    col = "white",
    lwd = 0.5,
    fill_alpha = 0.5
  )


mp1

```

Tmap allows to easily combine different objects by defining a new object via `tm_shape()`.

```{r}

# Define colours
cols <- viridis(n = 7, direction = 1, option = "C")

mp1 <- tm_shape(msoa.spdf) +
  tm_polygons(
    fill = "no2",
    fill_alpha = 1,
    fill.scale = tm_scale_intervals(
      style = "fisher",
      n = 7,
      values = cols
    ),
    fill.legend = tm_legend(
      title = "NO2",
      hist = FALSE
    )
  ) +
  tm_borders(
    col = "white",
    lwd = 0.5,
    fill_alpha = 0.5
  ) +
  tm_shape(ulez.spdf) +
  tm_borders(
    col = "red",
    lwd = 1,
    fill_alpha = 1
  )


mp1

```

And it is easy to change the layout.

```{r}

# Define colours
cols <- viridis(n = 7, direction = 1, option = "C")

mp1 <- tm_shape(msoa.spdf) +
  tm_polygons(
    fill = "no2",
    fill_alpha = 1,
    fill.scale = tm_scale_intervals(
      style = "fisher",
      n = 7,
      values = cols
    ),
    fill.legend = tm_legend(
      title = expression('in'~mu*'g'/m^{3}),
      hist = FALSE
    )
  ) +
  tm_borders(
    col = "white",
    lwd = 0.5,
    fill_alpha = 0.5
  ) +
  tm_shape(ulez.spdf) +
  tm_borders(
    col = "red",
    lwd = 1,
    fill_alpha = 1
  ) +
  tm_layout(
    frame = FALSE,
    legend.frame = TRUE,
    legend.bg.color = "white",
    legend.position = c("right", "bottom"),
    legend.outside = FALSE,
    legend.title.size = 0.8,
    legend.text.size = 0.8
  ) +
  tm_title_out(
    text = "NO2",
    position = tm_pos_out("center", "top",      # top cell
                          pos.h = "center"),    # position in cell
    size = 1.6
  )

mp1


```


We can also add some map information from OSM. However, it's sometimes a bit tricky with the projection. That's why we switch into the OSM projection here. Note that this osm query is build on retiring packages.

```{r, warning=FALSE}
# Save old projection
crs_orig <- st_crs(msoa.spdf)

# Change projection
ulez.spdf <- st_transform(ulez.spdf, 4326)
msoa.spdf <- st_transform(msoa.spdf, 4326)

# Get OSM data for background
osm_tmp <- read_osm(st_bbox(msoa.spdf), ext = 1.1, type = "osm-german") 

# Define colours
cols <- viridis(n = 7, direction = 1, option = "C")

mp1 <- tm_shape(osm_tmp) + 
  tm_rgb() +
  
  tm_shape(msoa.spdf) + 
  tm_polygons(
    fill = "no2",
    fill_alpha = 0.8,
    fill.scale = tm_scale_intervals(
      style = "fisher",
      n = 7,
      values = cols
    ),
    fill.legend = tm_legend(
      title = expression('in'~mu*'g'/m^{3}),
      hist = FALSE
    )
  ) +
  tm_shape(ulez.spdf) +
  tm_borders(
    col = "red",
    lwd = 1,
    fill_alpha = 1
  ) +
  tm_layout(
    frame = FALSE,
    legend.frame = TRUE,
    legend.bg.color = "white",
    legend.position = c("right", "bottom"),
    legend.outside = FALSE,
    legend.title.size = 0.8,
    legend.text.size = 0.8
  ) +
  tm_title_out(
    text = "NO2",
    position = tm_pos_out("center", "top", pos.h = "center"),
    size = 1.6
  )

mp1


```

Tmap also makes it easy to combine single maps

```{r}
# Define colours
cols1 <- viridis(n = 7, direction = 1, option = "C")

# Define colours
cols2 <- viridis(n = 7, direction = 1, option = "D")

# First map: NO2
mp1 <- tm_shape(osm_tmp) + 
  tm_rgb() +
  tm_shape(msoa.spdf) + 
  tm_polygons(
    fill = "no2",
    fill_alpha = 0.8,
    fill.scale = tm_scale_intervals(
      style = "fisher",
      n = 7,
      values = cols1
    ),
    fill.legend = tm_legend(
      title = expression('in'~mu*'g'/m^{3}),
      hist = FALSE
    )
  ) +
  tm_shape(ulez.spdf) +
  tm_borders(
    col = "red",
    lwd = 1,
    fill_alpha = 1
  ) +
  tm_layout(
    frame = FALSE,
    legend.frame = TRUE,
    legend.bg.color = "white",
    legend.position = c("right", "bottom"),
    legend.outside = FALSE,
    legend.title.size = 0.8,
    legend.text.size = 0.8
  ) +
  tm_title_out(
    text = "NO2",
    position = tm_pos_out("center", "top", pos.h = "center"),
    size = 1.4
  )

# Second map: % Black
mp2 <- tm_shape(osm_tmp) + 
  tm_rgb() +
  tm_shape(msoa.spdf) + 
  tm_polygons(
    fill = "per_black",
    fill_alpha = 0.8,
    fill.scale = tm_scale_intervals(
      style = "fisher",
      n = 7,
      values = cols2
    ),
    fill.legend = tm_legend(
      title = "% black",
      hist = FALSE
    )
  ) +
  tm_shape(ulez.spdf) +
  tm_borders(
    col = "red",
    lwd = 1,
    fill_alpha = 1
  ) +
  tm_layout(
    frame = FALSE,
    legend.frame = TRUE,
    legend.bg.color = "white",
    legend.position = c("right", "bottom"),
    legend.outside = FALSE,
    legend.title.size = 0.8,
    legend.text.size = 0.8
  ) +
  tm_title_out(
    text = "Ethnic Black inhabitants",
    position = tm_pos_out("center", "top", pos.h = "center"),
    size = 1.4
  )

# Arrange side by side
tmap_arrange(mp1, mp2, ncol = 2, nrow = 1)
```

And you can easily export those to png or pdf

```{r}
png(file = paste("London.png", sep = ""), width = 14, height = 7, units = "in", 
    res = 100, bg = "white")
par(mar=c(0,0,3,0))
par(mfrow=c(1,1),oma=c(0,0,0,0))
tmap_arrange(mp1, mp2, ncol = 2, nrow = 1)
dev.off()
```



## Exercises

1) What is the difference between a spatial "sf" object and a conventional "data.frame"? What's the purpose of the function `st_drop_geometry()`?

It's the same. A spatial "sf" object just has an additional column containing the spatial coordinates.

2) Using msoa.spdf, please create a spatial data frame that contains only the MSOA areas that are within the ulez zone.

```{r}
sub4.spdf <- msoa.spdf[ulez.spdf, ]
```


3) Please create a map for London (or only the msoa-ulez subset) which shows the share of Asian residents (or any other ethnic group).

```{r}
gp <- ggplot(msoa.spdf)+
    geom_sf(aes(fill = per_asian))+
    scale_fill_viridis_c(option = "E")+
    coord_sf(datum = NA)+
    theme_map()+
    theme(legend.position = c(.9, .6))
gp
```


4) Please calculate the distance of each MSOA to the London city centre 
  a) use google maps to get lon and lat, 
  b) use `st_as_sf()` to create the spatial point
  c) use `st_distance()` to calculate the distance
  
```{r}
### Distance to city center
# Define centre
centre <- st_as_sf(data.frame(lon = -0.128120855701165, 
                              lat = 51.50725909644806),
                   coords = c("lon", "lat"), 
                   crs = 4326)
# Reproject
centre <- st_transform(centre, crs = st_crs(msoa.spdf))
# Calculate distance
msoa.spdf$dist_centre <- as.numeric(st_distance(msoa.spdf, centre)) / 1000
# hist(msoa.spdf$dist_centre)
```
  
  
5) Can you create a plot with the distance to the city centre and pub counts next to each other?  


```{r}

# Define colours
cols <- viridis(n = 10, direction = 1, option = "B")
cols2 <- viridis(n = 10, direction = 1, option = "E")


library(tmap)

# Map 1: Distance - Fisher
mp1 <- tm_shape(msoa.spdf) + 
  tm_polygons(
    fill = "dist_centre",
    fill_alpha = 1,
    fill.scale = tm_scale_intervals(
      style = "fisher",
      n = 10,
      values = cols
    ),
    fill.legend = tm_legend(
      title = "Distance",
      hist = FALSE
    )
  ) +
  tm_borders(col = "white", lwd = 0.5, fill_alpha = 0.5) +
  tm_layout(
    frame = FALSE,
    legend.frame = TRUE,
    legend.bg.color = "white",
    legend.position = c("right", "bottom"),
    legend.outside = FALSE,
    legend.title.size = 0.8,
    legend.text.size = 0.8
  ) +
  tm_title_out(
    text = "Dist centre",
    position = tm_pos_out("center", "top", pos.h = "center"),
    size = 1.6
  )

# Map 2: Distance - Quantile
mp2 <- tm_shape(msoa.spdf) + 
  tm_polygons(
    fill = "dist_centre",
    fill_alpha = 1,
    fill.scale = tm_scale_intervals(
      style = "quantile",
      n = 10,
      values = cols
    ),
    fill.legend = tm_legend(
      title = "Distance",
      hist = FALSE
    )
  ) +
  tm_borders(col = "white", lwd = 0.5, fill_alpha = 0.5) +
  tm_layout(
    frame = FALSE,
    legend.frame = TRUE,
    legend.bg.color = "white",
    legend.position = c("right", "bottom"),
    legend.outside = FALSE,
    legend.title.size = 0.8,
    legend.text.size = 0.8
  ) +
  tm_title_out(
    text = "Dist centre",
    position = tm_pos_out("center", "top", pos.h = "center"),
    size = 1.6
  )

# Map 3: Pubs - Fisher
mp3 <- tm_shape(msoa.spdf) + 
  tm_polygons(
    fill = "pubs_count",
    fill_alpha = 1,
    fill.scale = tm_scale_intervals(
      style = "fisher",
      n = 10,
      values = cols
    ),
    fill.legend = tm_legend(
      title = "Count",
      hist = FALSE
    )
  ) +
  tm_borders(col = "white", lwd = 0.5, fill_alpha = 0.5) +
  tm_layout(
    frame = FALSE,
    legend.frame = TRUE,
    legend.bg.color = "white",
    legend.position = c("right", "bottom"),
    legend.outside = FALSE,
    legend.title.size = 0.8,
    legend.text.size = 0.8
  ) +
  tm_title_out(
    text = "Pubs",
    position = tm_pos_out("center", "top", pos.h = "center"),
    size = 1.6
  )

# Map 4: Pubs - Quantile
mp4 <- tm_shape(msoa.spdf) + 
  tm_polygons(
    fill = "pubs_count",
    fill_alpha = 1,
    fill.scale = tm_scale_intervals(
      style = "quantile",
      n = 10,
      values = cols
    ),
    fill.legend = tm_legend(
      title = "Count",
      hist = FALSE
    )
  ) +
  tm_borders(col = "white", lwd = 0.5, fill_alpha = 0.5) +
  tm_layout(
    frame = FALSE,
    legend.frame = TRUE,
    legend.bg.color = "white",
    legend.position = c("right", "bottom"),
    legend.outside = FALSE,
    legend.title.size = 0.8,
    legend.text.size = 0.8
  ) +
  tm_title_out(
    text = "Pubs",
    position = tm_pos_out("center", "top", pos.h = "center"),
    size = 1.6
  )



tmap_arrange(mp1, mp2, mp3, mp4, ncol = 2, nrow = 2)
```