diff --git a/Makefile b/Makefile
index 6e05756..e69b36e 100644
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,7 @@ all: data containers cycle
 
 data: osrm-data brouter-data
 containers: osrm-container brouter-container
-cycle: cycle_brouter
+cycle: WI-schools-cycle
 
 walk: route_analysis.Rmd
 	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./route_analysis.Rmd", output_file = "./html/route_analysis.html")'
@@ -13,6 +13,12 @@ cycle_osrm: cycling_route_analysis.Rmd
 cycle_brouter: cycling_route_analysis_brouter.Rmd
 	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./cycling_route_analysis_brouter.Rmd", output_file = "./html/cycling_route_analysis.html")'
 
+route_to_school: route_to_school.Rmd
+	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./route_to_school.Rmd", output_file = "./html/route_to_school.html")'
+
+WI-schools-cycle: WI-schools-cycle.Rmd
+	R -e 'library("rmarkdown"); old_path <- Sys.getenv("PATH"); Sys.setenv(PATH = paste(old_path, "/usr/local/bin", sep = ":")); rmarkdown::render(knit_root_dir = "./", output_dir = "./html", input = "./WI-schools-cycle.Rmd", output_file = "./html/WI-schools-cycle.html")'
+
 osrm-container: ./docker/osrm/docker-compose.yml
 	cd ./docker/osrm/; docker compose up -d
 
diff --git a/R/data/.gitignore b/R/data/.gitignore
new file mode 100644
index 0000000..5e7d273
--- /dev/null
+++ b/R/data/.gitignore
@@ -0,0 +1,4 @@
+# Ignore everything in this directory
+*
+# Except this file
+!.gitignore
diff --git a/R/functions.R b/R/functions.R
index 1e9d050..b6dc7d6 100644
--- a/R/functions.R
+++ b/R/functions.R
@@ -32,3 +32,60 @@ getLTSForRoute <- function(i, route_table) {
 
   return(result)
 }
+
+
+routeChar <- function(route){
+
+    if(is.na(route$messages)){
+        return(NA)
+    }
+
+    text <- route$messages
+    text <- gsub(x = text, pattern = "\\\"", replacement = "")
+    text <- gsub(x = text, pattern = "\ ", replacement = "")
+    text <- gsub(x = text, pattern = "\\[\\[", replacement = "")
+    text <- gsub(x = text, pattern = "\\]\\]", replacement = "")
+    foobar <- strsplit(text, split = "],[", fixed = TRUE)
+    x <- lapply(foobar, function(x){strsplit(x, split = ",", fixed = TRUE)})
+    xx <- unlist(x)
+    m <- matrix(xx, ncol = 13, byrow = TRUE)
+    names.vec <- m[1,]
+
+    if(nrow(m) == 2){
+        df <- data.frame(t(m[-1,]))
+    }else{
+        df <- data.frame(m[-1,])
+    }
+
+    names(df) <- names.vec
+
+    df2 <- within(df, {
+        Time <- as.numeric(Time)
+        stageTime <- diff(c(0,Time))
+        path <- grepl("highway=path", df$WayTags)
+        residential <- grepl("highway=residential", df$WayTags)
+        footway <- grepl("highway=footway", df$WayTags)
+        primary <- grepl("highway=primary", df$WayTags)
+        service <- grepl("highway=service", df$WayTags)
+        cycleway <- grepl("highway=cycleway", df$WayTags)
+        bike <- grepl("bicycle=designated", df$WayTags)
+    })
+
+
+    foo <- function(x){
+        ifelse(x$path, "path", ifelse(x$residential, "residential", ifelse(x$footway, "footway", ifelse(x$primary, "primary", ifelse(x$service, "service", ifelse(x$cycleway, "cycleway", "other"))))))
+    }
+
+    df2 <- cbind(df2, highway = foo(df2))
+    df2 <- df2 %>% group_by(highway) %>% summarize(T = sum(stageTime))
+
+    df2 <- df2 %>% filter(!is.na(highway))
+
+
+    if(!("cycleway" %in% df2$highway)){
+        return(0)
+    }else{
+        return(df2[df2$highway == "cycleway",]$T)
+    }
+
+}
diff --git a/WI-schools-cycle.Rmd b/WI-schools-cycle.Rmd
new file mode 100644
index 0000000..3f8b5aa
--- /dev/null
+++ b/WI-schools-cycle.Rmd
@@ -0,0 +1,233 @@
+---
+title: "Wisconsin Cycling to School"
+output:
+  html_document:
+    toc: true
+    toc_depth: 5
+    toc_float:
+      collapsed: false
+      smooth_scroll: true
+editor_options:
+  chunk_output_type: console
+---
+
+# Input Data & Configuration
+
+## Libraries
+
+```{r libs, eval = TRUE, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+date()
+rm(list=ls())
+library(tidyverse)
+library(ggmap)
+library(sf)
+library(osrm)
+library(smoothr)
+library(magick)
+library(ggnewscale)
+library(rsvg)
+library(httr)
+library(jsonlite)
+library(parallel)
+fig.height <- 6
+set.seed(1)
+source("./R/functions.R")
+runLoop <- FALSE
+```
+
+## Configuration
+
+```{r config, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+options(osrm.server = "http://127.0.0.1:5001/")
+options(osrm.profile = "bike")
+brouter_url <- "http://127.0.0.1:17777/brouter"
+brouter_profile <- "safety"
+register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
+WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg"), crs = 4326)
+WI_schools <- WI_schools %>% mutate(geom = SHAPE)
+```
+
+## Subset Schools
+
+```{r subsetSchools, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+WI_schools <- subset(WI_schools, !is.na(LAT) & !is.na(LON) & GRADE_RANGE == "09-12")
+```
+
+Non-virtual, grades 9-12.
+
+## Loop through WI Schools
+
+```{r mainloop, eval = runLoop, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = TRUE}
+radius <- 3 # miles
+levels <- c(1)
+res <- 100
+threshold <- units::set_units(1, km^2)
+
+gridList <- list()
+routesList <- list()
+
+indexVec <- 1:nrow(WI_schools)
+
+jj <- 1;
+bad.school.vec <- c()
+
+for(j in indexVec){
+
+    school_location <- WI_schools[j,]
+
+    cycle_boundary_m <- radius*1609
+    cycle_boundary_poly <- osrmIsodistance( loc = school_location, breaks = cycle_boundary_m, res = res )
+    cycle_boundary_poly <- st_make_valid(cycle_boundary_poly)
+    cycle_boundary_poly <- fill_holes(cycle_boundary_poly, threshold)
+    cycle_boundary_poly <- st_transform(cycle_boundary_poly, crs = 4326)
+
+    cellsize <- 1e-2
+    grid <- st_make_grid(cycle_boundary_poly, cellsize = cellsize, what = "polygons", square = FALSE)
+    grid <- st_intersection(cycle_boundary_poly, grid)
+
+    grid <- st_make_valid(grid)
+
+    grid_pts <- st_centroid(grid)
+
+    grid_coods <- st_coordinates(grid_pts)
+    school_focus_location <- school_location %>% select(LAT, LON)
+
+    routes <- list(NULL)
+    for(i in 1:nrow(grid_coods) ) {
+        query <- paste0(
+            brouter_url,
+            "?lonlats=", grid_coods[i,1], ",",grid_coods[i,2], "|",
+            school_focus_location$LON, ",", school_focus_location$LAT,
+            "&profile=", brouter_profile,
+            "&alternativeidx=0&format=geojson"
+        )
+        response <- GET(query)
+
+        if( response$status_code == "200" ){
+            route_run <- st_read(content <- content(response, as = "text"), quiet = TRUE)
+            routes[[i]] <- route_run
+        }else{
+            routes[[i]] <- NA
+        }
+    }
+
+    bad.cell <- which(is.na(routes))
+
+    if(length(bad.cell) > 0){
+        routes <- routes[-bad.cell]
+        grid <- grid[-bad.cell,]
+    }
+
+    if(length(routes) > 0){
+        routes <- st_transform(bind_rows(routes), crs = 4326)
+        gridList[[jj]] <- grid
+        routesList[[jj]] <- routes
+        jj <- jj + 1
+    }else{
+        routes <- NA
+        bad.school.vec <- c(j, bad.school.vec)
+        cat( WI_schools$SCHOOL[j], "has zero routes to school and has been removed from analysis.\n")
+    }
+
+}
+
+if(length(bad.school.vec) > 0){
+    WI_schools <- WI_schools[-bad.school.vec,]
+}
+
+saveRDS(WI_schools, "./R/data/WI_schools.rds")
+saveRDS(gridList, "./R/data/gridList.rds")
+saveRDS(routesList, "./R/data/routesList.rds")
+```
+
+## Read List Data
+
+```{r readLists, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+routesList <- readRDS(file = "./R/data/routesList.rds")
+gridList <- readRDS(file = "./R/data/gridList.rds")
+WI_schools <- readRDS(file = "./R/data/WI_schools.rds")
+
+not.cycleway.vec <- c()
+
+for(j in 1:length(gridList)){
+
+    grid <- gridList[[j]]
+    routes <- routesList[[j]]
+
+    total.time.vec <- routes %>% pull(total.time)
+    grid <- cbind(grid, total.time = as.numeric(total.time.vec))
+
+    x.vec <- c()
+    for( i in 1:nrow(grid) ){
+        route <- routes[i,"messages"] # Grid cell i to school j
+        x <- routeChar(route)
+        x.vec <- c(x.vec, x)
+    }
+
+    grid <- cbind(grid, T.cycleway = x.vec)
+    grid <- cbind( grid, not.cycleway = (grid$total.time - grid$T.cycleway)/60)
+    gridList[[j]] <- grid
+
+    not.cycleway.vec <- c(not.cycleway.vec, median(grid$not.cycleway))
+
+}
+```
+
+## Plot List Data
+
+### Best & Worst Schools
+
+```{r plots, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+ggplot(data.frame(not.cycleway = not.cycleway.vec), aes(not.cycleway)) + geom_histogram(fill = "grey", color = "black") + theme_bw()
+
+register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
+
+zoom.level <- 15
+k <- 306 #247
+
+bbox <- st_bbox(st_buffer(gridList[[k]], dist = 500))
+bbox <- c(left = as.double(bbox[1]),
+          bottom = as.double(bbox[2]),
+          right = as.double(bbox[3]),
+          top = as.double(bbox[4]))
+basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
+
+ggmap(basemap) + geom_sf(data = gridList[[k]], aes(fill= not.cycleway), inherit.aes = FALSE) + scale_fill_gradient(low = "yellow", high = "red", limits = c(0,17), na.value = NA)
+
+k <- 247 #306
+
+bbox <- st_bbox(st_buffer(gridList[[k]], dist = 500))
+bbox <- c(left = as.double(bbox[1]),
+          bottom = as.double(bbox[2]),
+          right = as.double(bbox[3]),
+          top = as.double(bbox[4]))
+basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
+
+ggmap(basemap) + geom_sf(data = gridList[[k]], aes(fill= not.cycleway), inherit.aes = FALSE) + scale_fill_gradient(low = "yellow", high = "red", limits = c(0,17), na.value = NA)
+```
+
+### Statewide Map
+
+```{r plots2, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+
+D <- cbind(WI_schools, not.cycleway.vec)
+
+zoom.level <- 8
+
+bbox <- st_bbox(st_buffer(D, dist = 10e3))
+bbox <- c(left = as.double(bbox[1]),
+          bottom = as.double(bbox[2]),
+          right = as.double(bbox[3]),
+          top = as.double(bbox[4]))
+basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
+
+ggmap(basemap) + geom_sf(data = D, aes(size = 2, color = not.cycleway.vec), inherit.aes = FALSE) + scale_color_gradient(low = "yellow", high = "red", na.value = NA) # , limits = c(0,17)
+
+```
+
+# Archive
+
+```{r chunklast, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+date()
+sessionInfo()
+```
diff --git a/route_to_school.Rmd b/route_to_school.Rmd
new file mode 100644
index 0000000..c4db7e2
--- /dev/null
+++ b/route_to_school.Rmd
@@ -0,0 +1,192 @@
+---
+title: "East High Cycling Routes"
+output:
+  html_document:
+    toc: true
+    toc_depth: 5
+    toc_float:
+      collapsed: false
+      smooth_scroll: true
+editor_options:
+  chunk_output_type: console
+---
+
+# Input Data & Configuration
+
+## Libraries
+
+```{r libs, eval = TRUE, echo = FALSE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+date()
+rm(list=ls())
+library(tidyverse)
+library(ggmap)
+library(sf)
+library(osrm)
+library(smoothr)
+library(magick)
+library(ggnewscale)
+library(rsvg)
+library(httr)
+library(jsonlite)
+library(parallel)
+fig.height <- 6
+set.seed(1)
+source("./R/functions.R")
+```
+
+# External sources configurations
+
+## Open Source Routing Machine (OSRM)
+
+```{r osrm, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+options(osrm.server = "http://127.0.0.1:5001/")
+options(osrm.profile = "bike")
+```
+
+## Brouter options
+```{r brouter, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+brouter_url <- "http://127.0.0.1:17777/brouter"
+brouter_profile <- "safety"
+```
+
+## Stadia Maps API Key
+
+```{r stadiamaps, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+register_stadiamaps(key = substr(read_file(file = "api_keys/stadia_api_key"), 1, 36))
+```
+
+
+# Analysis
+
+## Create Bikeable Region Using OSRM
+
+```{r boundary, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+WI_schools <- st_transform(st_read(dsn = "data/Schools/Wisconsin_Public_Schools_-5986231931870160084.gpkg"), crs = 4326)
+WI_schools <- WI_schools %>% mutate(geom = SHAPE)
+
+school_focus <- data.frame(name = c("East High School"), NCES_CODE = c("550852000925"))
+#school_focus <- data.frame(name = c("IMAP"), NCES_CODE = c("550008203085"))
+school_location <- WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE)
+
+radius <- 4 # miles
+levels <- c(1)
+res <- 100
+threshold <- units::set_units(1, km^2)
+
+cycle_boundary_m <- radius*1609
+
+cycle_boundary_poly <- osrmIsodistance( loc = school_location, breaks = cycle_boundary_m, res = res )
+
+cycle_boundary_poly <- st_make_valid(cycle_boundary_poly)
+
+cycle_boundary_poly <- fill_holes(cycle_boundary_poly, threshold)
+
+cycle_boundary_poly <- st_transform(cycle_boundary_poly, crs = 4326)
+
+saveRDS(cycle_boundary_poly, "./R/data/cycle_boundary_poly.rds")
+```
+
+# Create Grid Over Bikeable Region
+
+```{r grid, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+cellsize <- 5e-3
+grid <- st_make_grid(cycle_boundary_poly, cellsize = cellsize, what = "polygons", square = FALSE)
+grid <- st_intersection(cycle_boundary_poly, grid)
+```
+
+# Compute Routes from Cell Centroids to School with brouter
+
+```{r routes, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+grid_pts <- st_centroid(grid)
+grid_coods <- st_coordinates(grid_pts)
+school_focus_location <- WI_schools %>% filter(NCES_CODE %in% school_focus$NCES_CODE) %>% select(LAT, LON)
+
+routes <- list(NULL)
+for(i in 1:nrow(grid_coods) ) {
+  query <- paste0(
+    brouter_url,
+    "?lonlats=", grid_coods[i,1], ",",grid_coods[i,2], "|",
+    school_focus_location$LON, ",", school_focus_location$LAT,
+    "&profile=", brouter_profile,
+    "&alternativeidx=0&format=geojson"
+  )
+  response <- GET(query)
+  route_run <- st_read(content <- content(response, as = "text"), quiet = TRUE)
+  route_run[["student_number"]] <- i
+  routes[[i]] <- route_run
+}
+
+routes <- st_transform(bind_rows(routes), crs = 4326)
+```
+
+Notes:
+
+- What does `st_transform(bind_rows(routes), crs = 4326)` do?
+
+
+# Generate Map for Total Time
+
+## Set boundaries and get basemap
+```{r basemap, eval = TRUE, echo = TRUE, results = "show", warning = FALSE, error = TRUE, message = FALSE}
+bbox <- st_bbox(st_buffer(cycle_boundary_poly, dist = 500))
+bbox <- c(left = as.double(bbox[1]),
+          bottom = as.double(bbox[2]),
+          right = as.double(bbox[3]),
+          top = as.double(bbox[4]))
+
+zoom.level <- 12
+basemap <- get_stadiamap(bbox = bbox, zoom = zoom.level, maptype = "stamen_toner_lite")
+```
+
+## Total Trip Time Map
+
+```{r sandbox3, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+track.length.vec <- routes %>% pull(track.length)
+grid <- cbind(grid, track.length = as.numeric(track.length.vec)/1609)
+
+total.time.vec <- routes %>% pull(total.time)
+grid <- cbind(grid, total.time = as.numeric(total.time.vec))
+
+total.energy.vec <- routes %>% pull(total.energy)
+grid <- cbind(grid, total.energy = as.numeric(total.energy.vec))
+
+gg1 <- ggmap(basemap) + geom_sf(data = subset(grid, track.length > 1), aes(fill = total.time), inherit.aes = FALSE)
+
+ggsave(gg1, filename = "./figures/route-characteristics.pdf", width = 11, height = 8, units = "in")
+
+gg1
+```
+
+## Routes Map
+
+```{r sandbox3b, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+gg2 <- ggmap(basemap) + geom_sf(data = routes, aes(color = "red"), inherit.aes = FALSE)
+ggsave(gg2, filename = "./figures/routes.pdf", width = 11, height = 8, units = "in")
+gg2
+```
+
+# Route Characteristics
+
+## Compute Percent of Trip on Cycleway
+
+```{r sandbox4, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+x.vec <- c()
+for(j in 1:nrow(routes)){
+    foobar <- routeChar(routes[j, "messages"])
+    x.vec <- c(x.vec, foobar)
+}
+
+grid <- cbind(grid, T.cycleway = x.vec)
+grid <- cbind( grid, not.cycleway = (grid$total.time - grid$T.cycleway)/60)
+
+gg3 <- ggmap(basemap) + geom_sf(data = grid, aes(fill= not.cycleway), inherit.aes = FALSE) + scale_fill_gradient(low = "yellow", high = "red", limits = c(0,17), na.value = NA)
+ggsave(gg3, filename = "./figures/cycleway.pdf", width = 11, height = 8, units = "in")
+gg3
+```
+
+# Archive
+
+```{r chunklast, eval = TRUE, echo = TRUE, results = "show", warning = TRUE, error = TRUE, message = TRUE}
+date()
+sessionInfo()
+```