---
title: "In-Class Exercise 10"
date: "20 March 2023"
date-modified: "`r Sys.Date()`"
author: "Xinyee How"
format: html
editor: visual
execute: 
  echo: true
  eval: true
  warning: false
---

Installing packages

```{r}
pacman::p_load(tmap, SpatialAcc, sf, 
               ggstatsplot, reshape2,
               tidyverse)
```

Importing Geospatial data

```{r}
mpsz <- st_read(dsn = "data/geospatial", layer = "MP14_SUBZONE_NO_SEA_PL")
```

```{r}
hexagons <- st_read(dsn = "data/geospatial", layer = "hexagons") 
```

```{r}
eldercare <- st_read(dsn = "data/geospatial", layer = "ELDERCARE") 
```

Updating CRS information

```{r}
mpsz <- st_transform(mpsz, 3414)
eldercare <- st_transform(eldercare, 3414)
hexagons <- st_transform(hexagons, 3414)
```

Verify crs

```{r}
st_crs(mpsz)
```

Cleaning data

```{r}
eldercare <- eldercare %>%
  select(fid, ADDRESSPOS) %>%
  mutate(capacity = 100)
```

```{r}
hexagons <- hexagons %>%
  select(fid) %>%
  mutate(demand = 100)
```

Importing aspatial data

```{r}
ODMatrix <- read_csv("data/aspatial/OD_Matrix.csv", skip = 0)
```

Tidying OD matrix

```{r}
distmat <- ODMatrix %>%
  select(origin_id, destination_id, total_cost) %>%
  spread(destination_id, total_cost)%>%
  select(c(-c('origin_id')))
```

Changing metres to kilometres

```{r}
distmat_km <- as.matrix(distmat/1000)
```

Computing distance matrix

```{r}
eldercare_coord <- st_coordinates(eldercare)
hexagon_coord <- st_coordinates(hexagons)
```

```{r}
EucMatrix <- SpatialAcc::distance(hexagon_coord,
                                  eldercare_coord,
                                  type = "euclidean")
```

```{r}
EucMatrix_km <- as.matrix(EucMatrix/1000)
```

Modelling and Visualising Accessibility using Hansen Method

Computing Hansen's Method

```{r}
acc_Hansen <- data.frame(ac(hexagons$demand,
                            eldercare$capacity,
                            distmat_km, 
                            #d0 = 50,
                            power = 2, 
                            family = "Hansen"))
```

Tidying output - Dataframe

```{r}
colnames(acc_Hansen) <- "accHansen"
```

```{r}
acc_Hansen <- as_tibble(acc_Hansen)
```

Binding accHansen with hexagons

```{r}
hexagon_Hansen <- bind_cols(hexagons, acc_Hansen)
```

Alternatively, use this combined code chunk

```{r}
#| eval: false
acc_Hansen <- data.frame(ac(hexagons$demand,
                            eldercare$capacity,
                            distmat_km, 
                            #d0 = 50,
                            power = 0.5, 
                            family = "Hansen"))

colnames(acc_Hansen) <- "accHansen"
acc_Hansen <- tbl_df(acc_Hansen)
hexagon_Hansen <- bind_cols(hexagons, acc_Hansen)
```

Visualising Hansen's accessibility

Extract extend of hexagon for map plotting

```{r}
mapex <- st_bbox(hexagons)
```

Plotting map

```{r}
tmap_mode("plot")
tm_shape(hexagon_Hansen,
         bbox = mapex) + 
  tm_fill(col = "accHansen",
          n = 10,
          style = "quantile",
          border.col = "black",
          border.lwd = 1) +
tm_shape(eldercare) +
  tm_symbols(size = 0.1) +
  tm_layout(main.title = "Accessibility to eldercare: Hansen method",
            main.title.position = "center",
            main.title.size = 2,
            legend.outside = FALSE,
            legend.height = 0.45, 
            legend.width = 3.0,
            legend.format = list(digits = 6),
            legend.position = c("right", "top"),
            frame = TRUE) +
  tm_compass(type="8star", size = 2) +
  tm_scale_bar(width = 0.15) +
  tm_grid(lwd = 0.1, alpha = 0.5)
```
