midterm.Rmd

---
title: "675 Midterm"
author: "Kyle McCarthy"
date: "3/17/2021"
output: html_document
---

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

```{r setup, message=FALSE, warning=FALSE, include=TRUE, results='hide', cache=TRUE}


library(caret)
library(pscl)
library(plotROC)
library(pROC)
library(sf)
library(tidyverse)
library(knitr)
library(kableExtra)
library(FNN) 


mapTheme <- function(base_size = 12) {
  theme(
    text = element_text( color = "black"),
    plot.title = element_text(size = 14,colour = "black"),
    plot.subtitle=element_text(face="italic"),
    plot.caption=element_text(hjust=0),
    axis.ticks = element_blank(),
    panel.background = element_blank(),axis.title = element_blank(),
    axis.text = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(colour = "black", fill=NA, size=2)
  )
}

root.dir = "https://raw.githubusercontent.com/urbanSpatial/Public-Policy-Analytics-Landing/master/DATA/"
source("https://raw.githubusercontent.com/urbanSpatial/Public-Policy-Analytics-Landing/master/functions.r")

```


```{r ImportData, message=FALSE, warning=FALSE, include=TRUE, results='hide', cache=TRUE}

sr <- "+proj=utm +zone=12 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" 

Calgary <- st_read("C:/Users/Kyle McCarthy/Documents/CPLN 675/Midterm/midTermProject_Data/CALGIS_CITYBOUND_LIMIT/CALGIS_CITYBOUND_LIMIT.shp")%>% 
  st_transform(crs = sr)

inundation <- st_read("C:/Users/Kyle McCarthy/Documents/CPLN 675/Midterm/inundation.shp")%>% 
  st_transform(crs = sr)%>% 
  mutate(Label = ifelse(gridcode == 0, "No Flooding", "Flooding")) %>% 
  st_make_valid() 

imperviousness <- st_read("https://data.calgary.ca/resource/i9mt-qafb.geojson")%>%
  st_transform(crs = sr)

LandUse <- 
  st_read("C:/Users/Kyle McCarthy/Documents/CPLN 675/Midterm/LandUse.geojson")%>% 
  st_transform(crs = sr)%>% 
  filter(major == "Parks, Recreation and Public Education" | major == "Residential - Low Density")

PropVal <- st_read("C:/Users/Kyle McCarthy/Documents/CPLN 675/Midterm/midTermProject_Data/PropVal.geojson")%>% 
  st_transform(crs = sr) 
# const <- const_spatial %>%

fishnet <- 
  st_make_grid(Calgary,
               cellsize = 500) %>%
  st_sf() %>%
  mutate(uniqueID = rownames(.))%>%
  st_transform(crs = sr)



```


```{r ImportData, message=FALSE, warning=FALSE, include=TRUE, results='hide', cache=TRUE}


ggplot()+ 
  geom_sf(data = inundation, aes(fill = Label))+ 
  scale_fill_manual(values = c("#B5F8FE", "#0294A1"))+ 
  labs(title = "Calgary Flood Inundation Map")+
  theme(legend.position = "bottom")+
  mapTheme()



inundation_fishnet <- 
  inundation%>% 
  filter(gridcode == 1)%>% 
  st_intersection(fishnet, inundation)%>% 
  mutate(Area = as.numeric(st_area(.)))%>%
  st_drop_geometry()%>% 
  group_by(uniqueID)%>% 
  summarise(InundationArea = sum(Area))%>% 
  left_join(fishnet, .)%>% 
  mutate(pctInundation = InundationArea / 250000)%>% 
  mutate_all(funs(replace_na(.,0)))%>%
  mutate(Inundated = ifelse(pctInundation > .2, 1, 0))%>%
  mutate(Label = ifelse(Inundated == 1, "Inundated", "Not Inundated"))
  

```


```{r ImportData, message=FALSE, warning=FALSE, include=TRUE, results='hide', cache=TRUE}
ggplot()+ 
  geom_sf(data = inundation_fishnet, aes(fill = Label))+ 
  scale_fill_manual(values = c("#0294A1", "#B5F8FE"))+ 
  labs(title = "Calgary Flood Inundation Map")+
  theme(legend.position = "bottom")+
  mapTheme()


```


```{r ImportData, message=FALSE, warning=FALSE, include=TRUE, results='hide', cache=TRUE}

# Writing out fishnet as shapefile to perform zonal statistics as table on a euclidean distance raster from Calgary's water features 

# st_write(fishnet, "C:/Users/Kyle McCarthy/Documents/CPLN 675/Midterm/fishnet.shp")

# Reading zonal statistics Shapefile 

inundation_fishnet <- st_read("C:/Users/Kyle McCarthy/Documents/CPLN 675/Midterm/ZonalWater.dbf")%>% 
  dplyr::select(uniqueID, MEAN)%>%
  left_join(inundation_fishnet, .)%>% 
  na.omit()

ggplot()+ 
  geom_sf(data = inundation_fishnet, aes(fill = Label))+ 
  scale_fill_manual(values = c("#0294A1", "#B5F8FE"))+ 
  labs(title = "Calgary Flood Inundation Map")+
  theme(legend.position = "bottom")+
  mapTheme()

ggplot() + 
  geom_sf(data = inundation_fishnet, aes(fill = MEAN))+ 
  mapTheme()


ggplot()+ 
  geom_bar(data = inundation_fishnet, aes(x = Label, y = MEAN), stat="identity", fill = "blue")

```

``` {r Imperviousness, message=FALSE, warning=FALSE, include=TRUE, results='hide', cache=TRUE}

Impervious <- st_read("C:/Users/Kyle McCarthy/Documents/GitHub/CPLN-675-Midterm-/Data/Impervious/Imperv_cal.geojson")%>%
  st_transform(crs = sr)

# The Intersection Here Takes a While -- est. 20 min 

inundation_fishnet <- 
  Impervious%>% 
  filter(gen_surface == "Gravel" | gen_surface == "Buildings" | gen_surface == "Pavement" | gen_surface == "Bridge" | gen_surface == "Roads (Pavement)" )%>% 
  st_make_valid()%>%
  st_intersection(fishnet)%>% 
  mutate(Area = as.numeric(st_area(.)))%>% 
  st_drop_geometry()%>% 
  group_by(uniqueID)%>% 
  summarise(Area = sum(Area))%>% 
  mutate(pctImpervious = Area / 250000)%>% 
  left_join(inundation_fishnet, .)%>% 
  mutate(Imperv = ifelse(pctImpervious > 0.50, "Pervious", "Impermeable"))%>% 
  mutate_all(funs(replace_na(.,0)))

inundation_fishnet<- 
  inundation_fishnet%>% 
  mutate(Imperv = ifelse(pctImpervious > 0.50, "Impervious", "Permeable"), 
         pctImpervious = pctImpervious * 100) 

Impervioustats <- 
  inundation_fishnet%>%
  dplyr::select(pctImpervious, Label)%>% 
  mutate_all(funs(replace_na(.,0)))%>% 
  group_by(Label)%>% 
  summarise(pctImperv = mean(pctImpervious))



# Group for Graph 

ggplot()+ 
  geom_bar(data = Impervioustats, aes(x = Label, y = pctImperv), stat="identity", fill = "blue")




```


```{r Parks, message=FALSE, warning=FALSE, include=TRUE, results='hide', cache=TRUE}

st_c <- st_coordinates

inundation_Fish <- 
  PropVal %>% 
  mutate(assessed_value = as.numeric(assessed_value))%>%
  dplyr::select(assessed_value)%>%
  aggregate(., inundation_fishnet, mean)%>% 
  mutate_all(funs(replace_na(.,0)))%>% 
  mutate(uniqueID = rownames(.))%>%
  st_drop_geometry()%>% 
  left_join(inundation_fishnet, .)

PropValStats <- 
  inundation_Fish%>% 
  na.omit()%>%
  st_drop_geometry()%>% 
  group_by(Label)%>% 
  summarise(assessed_value = sum(assessed_value))

ggplot()+ 
  geom_bar(data = PropValStats, aes(x = Label, y = assessed_value), stat="identity", fill = "blue")



```

```{r LandUse, message=FALSE, warning=FALSE, include=TRUE, results='hide', cache=TRUE}

Inundation_Fishnet <- 
  LandUse%>% 
  st_make_valid()%>%
  st_intersection(inundation_Fish)%>% 
  mutate(landuArea = as.numeric(st_area(.)))%>% 
  st_drop_geometry()%>% 
  group_by(uniqueID)%>%
  summarise(landuArea = sum(landuArea))%>%
  left_join(inundation_Fish, .)%>%
  mutate(pctLandCover = landuArea / 250000)

LandStats <- 
  Inundation_Fishnet %>%
  mutate_all(funs(replace_na(.,0)))%>% 
  group_by(Label)%>% 
  summarise(Rec_low = mean(pctLandCover))


ggplot()+ 
  geom_bar(data =LandStats, aes(x = Label, y = Rec_low), stat="identity", fill = "blue")


```


```{r Building Denisty, message=FALSE, warning=FALSE, include=TRUE, results='hide', cache=TRUE}

Inundation_Fishnet<-st_read("C:/Users/Kyle McCarthy/Documents/GitHub/CPLN-675-Midterm-/Data/Buildings.geojson")%>%
  st_centroid()%>%
  st_transform(crs = sr)%>%
  mutate(BuildingCount = 1) %>%
  dplyr::select(BuildingCount)%>%
  aggregate(., inundation_fishnet, sum)%>%
  mutate_all(funs(replace_na(.,0)))%>%
  mutate(uniqueID = rownames(.))%>%
  st_drop_geometry()%>%
  left_join(Inundation_Fishnet, .)


BuidlingStats<-
  Inundation_Fishnet%>%
  st_drop_geometry()%>%
  na.omit()%>%
  group_by(Label)%>%
  summarise(BuildingSum = sum(BuildingCount),
            BuildingCountAverage = mean(BuildingCount))


ggplot()+
  geom_bar(data = BuidlingStats, aes(x = Label, y = BuildingCountAverage), stat="identity", fill = "blue")


```

```{r model, message=FALSE, warning=FALSE, include=TRUE, results='hide', cache=TRUE}


inundation_Fish <- 
  Inundation_Fishnet %>%
  dplyr::select(MEAN, pctImpervious, BuildingCount, pctLandCover, assessed_value, Inundated)%>%
  na.omit()%>% 
  st_intersection(inundationArea)%>% 
  mutate(uniqueID = rownames(.))%>%
  dplyr::select(-uniqueID)
  



inundation_Fish$Inundated <- as.factor(inundation_Fish$Inundated)

set.seed(3456)
trainIndex <- createDataPartition(inundation_Fish$Inundated, p = .70,
                                  list = FALSE,
                                  times = 1)
preserveTrain <- inundation_Fish[ -trainIndex,]
preserveTest  <- inundation_Fish[-trainIndex,]

typeof(inundation_Fish) 

preserveModel <- glm(Inundated ~ ., 
                    family="binomial"(link="logit"), data = preserveTrain %>%
                                                            as.data.frame() %>%
                                                            dplyr::select(-geometry))
summary(preserveModel)

```


```{r Model Validation, message=FALSE, warning=FALSE, include=TRUE, cache=TRUE, fig.height = 12, fig.width = 12} 
classProbs <- predict(preserveModel, preserveTest, type="response")

hist(classProbs)

testProbs <- data.frame(obs = as.numeric(preserveTest$Inundated),
                        pred = classProbs)
ggplot(testProbs, aes(x = pred, fill=as.factor(obs))) + geom_density() +
  facet_grid(obs ~ .) + xlab("Probability") + geom_vline(xintercept = .5) +
  scale_fill_manual(values = c("dodgerblue4", "darkgreen"),
                      labels = c("Not Inundated","Inundated"),
                      name = "")
```

```{r ROC message=FALSE, warning=FALSE, include=TRUE, cache=TRUE, fig.height = 12, fig.width = 12}
testProbs$predClass  = ifelse(testProbs$pred > 0.5 , 1,0)

testProbs <- 
  testProbs %>% 
  mutate(obs = ifelse(obs == 2, 1, 0))


caret::confusionMatrix(reference = as.factor(testProbs$obs), 
                       data = as.factor(testProbs$predClass), 
                       positive = "1")


ggplot(testProbs, aes(d = obs, m = pred)) + 
  geom_roc(n.cuts = 50, labels = FALSE) + 
  style_roc(theme = theme_grey) +
  geom_abline(slope = 1, intercept = 0, size = 1.5, color = 'grey')


preserveTest$uniqueID <- seq.int(nrow(preserveTest))
testProbs1 <- testProbs
testProbs$uniqueID <- seq.int(nrow(preserveTest))

try <- left_join(preserveTest, testProbs)

try <- 
  try%>% 
  mutate(result = "0")%>% 
  mutate(result = ifelse(obs == 0 & predClass == 0, "True Negative", result))%>% 
  mutate(result = ifelse(obs == 1 & predClass == 1, "True Positive", result))%>% 
  mutate(result = ifelse(obs == 0 & predClass == 1, "False Positive", result))%>% 
  mutate(result = ifelse(obs == 1 & predClass == 0, "False Negative", result))


ggplot() + 
  geom_sf(data = inundation_Fish)+ 
  geom_sf(data = try, aes(fill = result))+ 
  mapTheme() 
  

```



```{r Cross Validation, message=FALSE, warning=FALSE, include=TRUE, cache = TRUE, fig.height = 8 fig.width = 8} 

ctrl <- trainControl(method = "cv", 
                     number = 100, 
                     savePredictions = TRUE)

inundation_Fish1 <- 
  inundation_Fish %>% 
  na.omit() %>% 
  dplyr::select(MEAN, pctImpervious, BuildingCount, pctLandCover, assessed_value, Inundated)

cvFit <- train(as.factor(Inundated) ~ .,  data = inundation_Fish1 %>% 
                                                as.data.frame() %>%
                                                dplyr::select(-geometry), 
               method="glm", family="binomial",
               trControl = ctrl)

cvFit

ggplot(as.data.frame(cvFit$resample), aes(Accuracy)) + 
  geom_histogram() +
  scale_x_continuous(limits = c(0, 1)) +
  labs(x="Accuracy",
       y="Count")



```

```{r Mapping Predictions, message=FALSE, warning=FALSE, include=TRUE, cache = TRUE, fig.height = 12, fig.width = 12} 

allPredictions <- 
  predict(cvFit, inundation_Fish, type="prob")[,2]
  
preserve <- 
  cbind(inundation_Fish,allPredictions) %>%
  mutate(allPredictions = round(allPredictions * 100)) 

 ggplot() + 
    geom_sf(data=preserve, aes(fill=factor(ntile(allPredictions,5))), colour=NA) +
    scale_fill_manual(values = c("#edf8fb","#b3cde3","#8c96c6","#8856a7","#810f7c"),
                      labels=as.character(quantile(preserve$allPredictions,
                                                 c(0.1,.2,.4,.6,.8),na.rm=T)),
                      name="Predicted\nProbabilities(%)\n(Quintile\nBreaks)") +
  mapTheme() +
  labs(title="")
 
 
 
 
```