Freight O-D.Rmd

---
title: "Freight O-D"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

Sys.setenv(PATH = paste("C:\\Python27\\ArcGIS10.7\\", Sys.getenv("PATH"), sep=";"))
```

# Tract-Tract O-D

```{r}
library(sf)
library(tidyverse)
library(odbc)

mn_tracts <- st_read("Shapefile/86096_Freight_Tract_to_Tract_10_16_19_1_origin.shp")

# Get ID of tracts, starting with state (27) designation
tracts_id <- mn_tracts %>%
  mutate(Name2 = name) %>%
  separate(Name2, into = c("Prefix", "Tract"), sep = -11)

# Get ID of metro tracts (no geometry)
gis_db <- dbConnect(odbc(), "GIS")
metro_tracts <- dbGetQuery(gis_db, "SELECT * FROM GISLibrary.dbo.CENSUS2010TIGERTRACT")
dbDisconnect(gis_db)

# Create Greater MN shapefile
metro_tracts_id_only <- metro_tracts %>%
  transmute(Tract = GEOID10)

greater_mn_tracts <- anti_join(tracts_id, metro_tracts_id_only, by = "Tract")

# Get metro tracts as shapefile
gis_db <- dbConnect(odbc(), "GIS")
metro_tracts_shp <- st_as_sf(dbGetQuery(gis_db, "SELECT *, Shape.STAsText() as wkt FROM GISLibrary.dbo.CENSUS2010TIGERTRACT"), wkt = "wkt", crs = 26915)
dbDisconnect(gis_db)

metro_dissolve <- st_union(metro_tracts_shp)

# Get centroids
greater_mn_centroids <- st_centroid(greater_mn_tracts)

metro_centroid <- st_centroid(metro_dissolve)

# od_centroids <- st_union(greater_mn_centroids, metro_centroid)
# 
# od_centroids_attr <- od_centroids %>%
#   mutate(Tract = ifelse(is.na(Tract), "Metro", Tract))

# st_write(od_centroids, "Greater_MN_Metro_Centroids.shp", "ESRI Shapefile")

st_write(greater_mn_centroids, "Greater MN Centroids.shp", "ESRI Shapefile")
st_write(metro_centroids, "Metro Centroid.shp", "ESRI Shapefile")
st_write(metro_dissolve, "Metro Outline.shp", "ESRI Shapefile")
```

# Using st_union creates shapefile that doesn't render correctly in Tableau with shape type (works fine with shapefile type) - merge the two shapefiles using Arc so CRS compatibility is fixed automatically

```{python}
import arcpy
from arcpy import env
env.overwriteOutput = True
env.workspace = "N:\MTS\Working\Planning\Highway Planning\Twin Cities Mobility Needs Analysis 2020\Freight-O-D"

# Convert shp to table view
arcpy.MakeTableView_management("Greater MN Centroids.shp", "greater_mn")
arcpy.MakeTableView_management("Metro Centroid.shp", "metro")

arcpy.Merge_management(["greater_mn", "metro"], "Greater_MN_Metro_Merge")

```

```{r}
greater_mn_metro <- st_read("Greater_MN_Metro_Merge.shp")

od_centroids <- greater_mn_metro %>%
  mutate(Tract = ifelse(is.na(Tract), "Metro", Tract))

st_write(od_centroids, "OD_Centroids.shp")
```

```{r}
library(data.table)

freight_od <- fread("86096_Freight_Tract_to_Tract_10_16_19_1_od_comm.csv")

freight_od_ids <- freight_od[, "Origin_tract" := substring(`Origin Zone Name`, 10, length(`Origin Zone Name`))][, "Destination_tract" := substring(`Destination Zone Name`, 10, length(`Destination Zone Name`))]

o_ids <- metro_tracts %>%
  transmute(O_tract = GEOID10,
            O_metro_desig = "Metro")

o_tr <- data.table(o_ids)

setkey(freight_od_ids, Origin_tract)
setkey(o_tr, O_tract)

freight_origins <- o_tr[freight_od_ids]

dest_ids <- metro_tracts %>%
  transmute(D_tract = GEOID10,
            D_metro_desig = "Metro")

d_tr <- data.table(dest_ids)

setkey(freight_origins, Destination_tract)
setkey(d_tr, D_tract)

freight_dests <- d_tr[freight_origins]

freight_od <- freight_dests[, "Dest_tract_or_metro" := ifelse(is.na(D_metro_desig), D_tract, D_metro_desig)][, "Orig_tract_or_metro" := ifelse(is.na(O_metro_desig), O_tract, O_metro_desig)]

fwrite(freight_od, "SL Results with Metro Designation.csv")

freight_od %>%
  head()
```

```{r}
# Get ID of metro tracts (no geometry)
gis_db <- dbConnect(odbc(), "GIS")
metro_tracts <- dbGetQuery(gis_db, "SELECT GEOID10 FROM GISLibrary.dbo.CENSUS2010TIGERTRACT")
dbDisconnect(gis_db)

metro_tr <- metro_tracts %>%
  transmute(Metro_desig = "Metro",
            Tract = GEOID10)

# Read in SL data
freight_od <- fread("86096_Freight_Tract_to_Tract_10_16_19_1_od_comm.csv")

freight_pared <- freight_od %>%
  filter(`Day Part` == "0: All Day (12am-12am)" & `Day Type` == "0: All Days (M-Su)") %>%
  dplyr::select(`Origin Zone Name`, `Destination Zone Name`, `Average Daily O-D Traffic (StL Index)`) %>%
  separate(`Origin Zone Name`, into = c("Prefix", "Origin_tract"), sep = -11) %>%
  select(-Prefix) %>%
  separate(`Destination Zone Name`, into = c("Prefix", "Dest_tract"), sep = -11) %>%
  select(-Prefix)

# Add column designating metro tracts as "Metro"
freight_o_desig <- left_join(freight_pared, metro_tr, by = c("Origin_tract" = "Tract"))
freight_d_desig <- left_join(freight_o_desig, metro_tr, by = c("Dest_tract" = "Tract"))

# Overwrite all metro tracts with single designation to collapse to metro area
freight_metro <- freight_d_desig %>%
  mutate(Origin_tract = ifelse(is.na(Metro_desig.x), Origin_tract, "Metro"),
         Dest_tract = ifelse(is.na(Metro_desig.y), Dest_tract, "Metro")) %>%
  rename(Avg_traffic = `Average Daily O-D Traffic (StL Index)`) %>%
  dplyr::select(Origin_tract, Dest_tract, Avg_traffic) %>%
  group_by(Origin_tract, Dest_tract) %>%
  mutate(Avg_traffic = sum(Avg_traffic)) %>%
  unique()

# Calculate total index going from origin
freight_origin <- freight_metro %>%
  group_by(Origin_tract) %>%
  mutate(Total_origin_stl = sum(Avg_traffic)) %>%
  select(Origin_tract, Total_origin_stl) %>% 
  unique()

# Calculate total index coming to destination
freight_dest <- freight_metro %>%
  group_by(Dest_tract) %>%
  mutate(Total_dest_stl = sum(Avg_traffic)) %>%
  select(Dest_tract, Total_dest_stl) %>% 
  unique()

dest_count <- freight_d_desig %>%
  mutate(Origin_tract = ifelse(is.na(Metro_desig.x), Origin_tract, "Metro"),
         Dest_tract = ifelse(is.na(Metro_desig.y), Dest_tract, "Metro")) %>%
  group_by(Dest_tract) %>%
  count() %>%
  rename(Count_dest_tract = n)

# Add totals to O-D matrix
freight_o <- left_join(freight_metro, freight_origin, by = "Origin_tract")
freight_d <- left_join(freight_o, freight_dest, by = "Dest_tract")
freight_count <- left_join(freight_d, dest_count, by = "Dest_tract")

od_perc <- freight_count %>%
  ungroup() %>%
  mutate(Perc_of_dest_traffic = Avg_traffic/Total_dest_stl*100,
         Path_id = row_number()) %>%
  gather(Origin_tract, Dest_tract, key = "Tract_category", value = "Tract")

write_excel_csv(od_perc, "Freight StL OD Percentages.csv")

  
```

# No Metro Attribution

```{r}
# Read in SL data
freight_od <- fread("86096_Freight_Tract_to_Tract_10_16_19_1_od_comm.csv")

freight_pared <- freight_od %>%
  filter(`Day Part` == "0: All Day (12am-12am)" & `Day Type` == "0: All Days (M-Su)") %>%
  dplyr::select(`Origin Zone Name`, `Destination Zone Name`, `Average Daily O-D Traffic (StL Index)`) %>%
  separate(`Origin Zone Name`, into = c("Prefix", "Origin_tract"), sep = -11) %>%
  select(-Prefix) %>%
  separate(`Destination Zone Name`, into = c("Prefix", "Dest_tract"), sep = -11) %>%
  select(-Prefix) %>%
  rename(Avg_traffic = `Average Daily O-D Traffic (StL Index)`)

# Calculate total index coming to destination
freight_dest_full <- freight_pared %>%
  group_by(Dest_tract) %>%
  mutate(Total_dest_stl = sum(Avg_traffic)) %>%
  select(Dest_tract, Total_dest_stl) %>% 
  unique()

dest_count_full <- freight_pared %>%
  group_by(Dest_tract) %>%
  count() %>%
  rename(Count_dest_tract = n)

# Add totals to O-D matrix
freight_d_full <- left_join(freight_pared, freight_dest_full, by = "Dest_tract")
freight_full_count <- left_join(freight_d_full, dest_count_full, by = "Dest_tract")

# O-D percentages
od_perc_full <- freight_full_count %>%
  mutate(Origin_tract2 = Origin_tract,
         Dest_tract2 = Dest_tract) %>%
  separate(Origin_tract2, into = c("Origin_county", "Origin_block"), sep = 5) %>%
  separate(Dest_tract2, into = c("Dest_county", "Dest_block"), sep = 5) %>%
  mutate(Origin_county = ifelse(Origin_county == 27003 |
                                  Origin_county == 27019 |
                                  Origin_county == 27037 |
                                  Origin_county == 27053 |
                                  Origin_county == 27123 |
                                  Origin_county == 27139 |
                                  Origin_county == 27163, "Metro", Origin_county),
         Dest_county = ifelse(Dest_county == 27003 |
                                  Dest_county == 27019 |
                                  Dest_county == 27037 |
                                  Dest_county == 27053 |
                                  Dest_county == 27123 |
                                  Dest_county == 27139 |
                                  Dest_county == 27163, "Metro", Dest_county),
         Metro_to_metro = ifelse(Origin_county == "Metro" & Dest_county == "Metro", 1, 0)) %>%
  filter(Metro_to_metro == 0) %>%
  ungroup() %>%
  mutate(Perc_of_dest_traffic = Avg_traffic/Total_dest_stl*100,
         Path_id = row_number()) %>%
  gather(Origin_tract, Dest_tract, key = "Tract_category", value = "Tract")

write_excel_csv(od_perc_full, "Freight StL OD Percentages No Metro-Metro.csv")


```

# County-County O-D

```{r}
library(sf)
library(tidyverse)
library(odbc)

# Get MN counties as shapefile
gis_db <- dbConnect(odbc(), "GIS")
mn_counties <- st_as_sf(dbGetQuery(gis_db, "SELECT *, Shape.STAsText() as wkt FROM GISLibrary.dbo.MNCOUNTIES"), wkt = "wkt", crs = 26915)
dbDisconnect(gis_db)

county_centroids <- st_centroid(mn_counties)

st_write(county_centroids, "MN_County_Centroids.shp")
```

```{r}
# MN Counties shapefile used for o-d analysis (contains fewer attributes, but suffices for visualization)

counties_od <- st_read("County O-D/Shapefile/87342_County_County_OD_origin.shp")

counties_od_centroids <- st_centroid(counties_od)

st_write(counties_od_centroids, "OD_County_Centroids.shp")
```

```{r}
county_ods <- read_csv("County O-D/87342_County_County_OD_od_comm.csv")

options(scipen=999)
county_ods_tidy <- county_ods %>%
  filter(`Vehicle Weight` == "Heavy") %>%
  group_by(`Destination Zone ID`) %>%
  mutate(Dest_tot_traffic = sum(`Average Daily O-D Traffic (StL Index)`)) %>%
  ungroup() %>%
  dplyr::select(`Origin Zone ID`, `Destination Zone ID`, Dest_tot_traffic, `Average Daily O-D Traffic (StL Index)`) %>%
  unique() %>%
  mutate(Perc_of_dest = `Average Daily O-D Traffic (StL Index)`/Dest_tot_traffic) %>%
  rename(Origin = `Origin Zone ID`,
         Destination = `Destination Zone ID`) %>%
  mutate(Path_ID = row_number()) %>%
  gather(Origin, Destination, key = "OD_Category", value = "Zone_ID")

county_od_count <- county_ods %>%
  filter(`Vehicle Weight` == "Heavy") %>%
  group_by(`Destination Zone ID`) %>%
  count() %>%
  rename(Count_dests = n,
         Zone_ID = `Destination Zone ID`) %>%
  mutate(OD_Category = "Destination")

county_od_w_count <- left_join(county_ods_tidy, county_od_count, by = c("OD_Category", "Zone_ID")) %>%
  mutate(Count_dests = ifelse(is.na(Count_dests), 0, Count_dests),
         Count_dests = as.numeric(Count_dests))

write_csv(county_od_w_count, "StL County O-D.csv")
```