---
title: "In-Class Exercise 11: Calibrating Spatial Interaction Models"
date: "27 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(sp, sf, spdep, GWmodel, tmap, carat, stplanr, reshape2, broom, tidyverse)
```

#### Reading aspatial data

```{r}
odbus <- read_csv("data/aspatial/origin_destination_bus_202210.csv")
```

**Converting numeric data into character data**

```{r}
#| eval: false
odbus$ORIGIN_PT_CODE <- 
  as.factor(odbus$ORIGIN_PT_CODE)
odbus$DESTINATION_PT_CODE <- 
  as.factor(odbus$DESTINATION_PT_CODE)
```

```{r}
#| eval: false
odbus7_9 <- odbus %>%
  filter(DAY_TYPE == "WEEKDAY") %>%
  filter(TIME_PER_HOUR >= 7 & TIME_PER_HOUR <= 9) %>%
  group_by(ORIGIN_PT_CODE, DESTINATION_PT_CODE) %>%
  summarise(TRIPS = sum(TOTAL_TRIPS))
```

**Saving output as data**

```{r}
#| eval: false
odbus7_9 <- write_rds(odbus7_9, "data/rds/odbus7_9.rds")
```

**Bringing in previously saved output**

```{r}
odbus7_9 <- read_rds("data/rds/odbus7_9.rds")
```

#### Importing geospatial data

**Bus stops**

```{r}
busstop = st_read(dsn = "data/geospatial", layer = "BusStop")
busstop <- busstop %>% st_transform(crs = 3414)
```

**Master Plan 2019**

```{r}
mpsz = st_read(dsn = "data/geospatial", layer = "MPSZ-2019")
mpsz <- mpsz %>% st_transform(crs = 3414)
```

**Combining bus stop and MPSZ data**

```{r}
busstop_mpsz <- st_intersection(busstop, mpsz) %>%
  select(BUS_STOP_N, SUBZONE_C) %>%
  st_drop_geometry()
```

## Geospatial Data Wrangling

**Appending planning subzone data onto odbus7_9 dataframe**

```{r}
#| eval: false
od_data <- left_join(odbus7_9, busstop_mpsz,
                     by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
  rename(ORIGIN_BS = ORIGIN_PT_CODE,
         ORIGIN_SZ = SUBZONE_C,
         DESTIN_BS = DESTINATION_PT_CODE)
```

**Check for duplicate records**

```{r}
#| eval: false
duplicate <- od_data %>% group_by_all() %>%
  filter(n()>1) %>% ungroup()
```

**Retaining unique records only**

```{r}
#| eval: false
od_data <- unique(od_data)
```

**Joining destination bus stop with planning subzone**

```{r}
#| eval: false
od_data <- left_join(od_data, busstop_mpsz,
                     by = c("DESTIN_BS" = "BUS_STOP_N"))
```

**Checking for duplicates again**

```{r}
#| eval: false
duplicate2 <- od_data %>% group_by_all() %>%
  filter(n()>1) %>% ungroup()
```

**Retaining unique data only**

```{r}
#| eval: false
od_data <- unique(od_data)
```

**Dropping missing values**

```{r}
#| eval: false
od_data <- od_data %>%
  rename(DESTIN_SZ = SUBZONE_C) %>%
  drop_na()
```

**Saving output**

```{r}
#| eval: false
write_rds(od_data, "data/rds/od_data.rds")
```

**Bringing in previously saved output**

```{r}
od_data <- read_rds("data/rds/od_data.rds")
```

## Visualising Geospatial data

Planning subzone

```{r}
tmap_mode("view")
tmap_options(check.and.fix = TRUE)
qtm(mpsz)
```

**Sorting by Subzone**

```{r}
mpsz <- mpsz[order(mpsz$SUBZONE_C),]
head(mpsz,10)
```

## Computing Distance Matrix

**Converting into sp version**

```{r}
mpsz_sp <- as(mpsz, "Spatial")
```

**Computing distance matrix**

```{r}
dist <- spDists(mpsz_sp)
dist #values are in metres
```

**Sorting distance matrix by planning subzone**

```{r}
sz_names <- mpsz$SUBZONE_C
```

**Attaching SUBZONE_C to row and column for distance matrix ahead**

```{r}
colnames(dist) <- paste0(sz_names)
rownames(dist) <- paste0(sz_names)
```

Pivoting distance value by SUBZONE C

```{r}
distPair <- melt(dist) %>%
  rename(dist = value)
head(distPair, 10)
```

```{r}
distPair$dist <- ifelse(distPair$dist == 0,
                        50, distPair$dist)
```

**Rename origin and destination fields**

```{r}
distPair <- distPair %>%
  rename(orig = Var1,
         dest = Var2)
```

**Completing the flow data**

```{r}
flow_data <- od_data %>%
  group_by(ORIGIN_SZ, DESTIN_SZ) %>%
  summarize(TRIPS = sum(TRIPS))
```

### Bringing in population data

```{r}
pop <- read_csv("data/aspatial/pop.csv")
```

```{r}
pop <- pop %>% left_join(mpsz,
                         by = c("PA" = "PLN_AREA_N",
                                "SZ" = "SUBZONE_N")) %>%
  select(1:6) %>%
  rename(SZ_NAME = SZ,
         SZ = SUBZONE_C)
```

```{r}
flow_data1 <- flow_data %>%
  left_join(distPair,
            by = c("ORIGIN_SZ" = "orig",
                   "DESTIN_SZ" = "dest"))
```

```{r}
flow_data1 <- flow_data1 %>%
  left_join(pop,
            by = c(ORIGIN_SZ = "SZ")) %>%
  rename(ORIGIN_AGE7_12 = AGE7_12,
         ORIGIN_AGE13_24 = AGE13_24,
         ORIGIN_AGE25_64 = AGE25_64) %>%
  select(-c(PA, SZ_NAME))
```

```{r}
flow_data1 <- flow_data1 %>%
  left_join(pop,
            by = c(DESTIN_SZ = "SZ")) %>%
  rename(DESTIN_AGE7_12 = AGE7_12,
         DESTIN_AGE13_24 = AGE13_24,
         DESTIN_AGE25_64 = AGE25_64) %>%
  select(-c(PA, SZ_NAME))
```
