CPLN 675 Midterm Project: Forecasting Flood Inundation
Huidi Hu
March 25, 2019
1 Motivation
The world is now exposed to more severe and prevalent flood risk than in the past decades. With climate change, the sea level is rising, and the weather is becoming more extreme. Such challenges demand an integrated flood risk management plan based on accurate and generalizable predictions, which requires a new method to predict flood risk in a timely and efficient manner. By establishing the flood inundation model, the government will be able to assign tax revenues and federal funds to the critical areas, as well as better coordinate regional flood infrastructures and future developments.
2 Data
First, load the requisite packages and create a mapTheme() function that is used to standardized the map outputs created below.
library(tidyverse)
library(sf)
library(QuantPsyc)
library(RSocrata)
library(viridis)
library(caret)
library(spatstat)
library(spdep)
library(FNN)
library(grid)
library(gridExtra)
library(knitr)
library(kableExtra)
library(plotROC)
library(pROC)
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)
)
}Download data and create base maps, including city boundary, flood map and fishnet.
limit <-
st_read("./midTermProject_Data/CALGIS_CITYBOUND_LIMIT/CALGIS_CITYBOUND_LIMIT.shp") %>%
st_transform(crs=3780)
fishnet <-
st_make_grid(limit, cellsize = 500) %>%
st_sf()
fishnet <-
fishnet[limit,] %>%
mutate(uniqueID = rownames(.)) %>%
dplyr::select(uniqueID)
river <- st_read("./calgary/river_r.shp") %>%
st_transform(crs=3780)
flood <- st_read("./calgary/flood_r.shp") %>%
st_transform(crs=3780)ggplot() +
geom_sf(data=limit) +
geom_sf(data=flood,color="red", fill="red") +
labs(title="Past Flood in Calgary") +
mapTheme()
Next, we gather, wrangle and plot the flood risk factor data.
library(raster)
landuse <- st_read("./calgary/landuse.shp") %>%
st_transform(crs=3780)
elevation <- st_read("./calgary/elevation.shp") %>%
st_transform(crs=3780)
drain <- st_read("./calgary/hydro/drainv.shp") %>%
st_transform(crs=3780)
tree <-
st_read("https://data.calgary.ca/resource/tfs4-3wwa.geojson") %>%
st_transform(crs=3780)
sloper <- raster("./calgary/sloper.tif")
slopep <-
rasterToPoints(sloper) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y"), crs = 3780)
steepp <- slopep %>%
filter(sloper > 3)
fac <- raster("./calgary/hydro/fac.tif")
facp <-
rasterToPoints(fac) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y"), crs = 3780)
band1 <- raster("./calgary/data/band1.tif")
band2 <- raster("./calgary/data/band2.tif")
band3 <- raster("./calgary/data/band3.tif")
band4 <- raster("./calgary/data/band4.tif")
band5 <- raster("./calgary/data/band5.tif")
band6 <- raster("./calgary/data/band6.tif")
band7 <- raster("./calgary/data/band7.tif")
band8 <- raster("./calgary/data/band8.tif")
band9 <- raster("./calgary/data/band9.tif")
band10 <- raster("./calgary/data/band10.tif")
band11 <- raster("./calgary/data/band11.tif")
band8 <- aggregate(band8, fact = 2)
image <- stack(band1, band2, band3, band4, band5, band6, band7,
band8, band9, band10, band11)
ndvi <- (image[[5]] - image[[4]])/(image[[5]] + image[[4]])
ndvip <-
rasterToPoints(ndvi) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y"), crs = st_crs(ndvi))%>%
st_transform(st_crs(limit))
greenp <- ndvip %>%
filter(layer > 0.35) %>%
st_transform(st_crs(limit))
streamp <- st_read("./calgary/hydro/streamp.shp")%>%
st_transform(crs=3780) We then join the flood map and the risk factors to the fishnet.
#unique(landuse$major)
tbc_net <-
landuse %>%
filter(major== "Parks, Recreation and Public Education" | major== "Future Urban Development") %>%
mutate(tbc =1) %>%
aggregate(., fishnet, mean) %>%
mutate(tbc = ifelse(is.na(tbc), 0, tbc))%>%
dplyr::select(tbc)
lowdens_net <-
landuse %>%
filter(major== "Residential - Low Density"|major== "Residential - Medium Density"|major== "Industrial"|major== "Transportation and Utility Corridor"|major== "Major Infrastructure") %>%
mutate(lowdens =1) %>%
aggregate(., fishnet, mean) %>%
mutate(lowdens = ifelse(is.na(lowdens), 0, lowdens))%>%
dplyr::select(lowdens)
highdens_net <-
landuse %>%
filter(major== "Residential - High Density"|major== "Commercial - Core"|major== "Mixed Use"|major== "Institutional"|major== "Recreational"|major== "Commercial") %>%
mutate(highdens =1) %>%
aggregate(., fishnet, mean) %>%
mutate(highdens = ifelse(is.na(highdens), 0, highdens))%>%
dplyr::select(highdens)
elevation_net <-
elevation %>%
aggregate(., fishnet, mean)%>%
mutate(GRIDCODE = ifelse(is.na(GRIDCODE), mean(na.omit(GRIDCODE)), GRIDCODE))
#plot(elevation_net)
slopep_net <-
slopep %>%
aggregate(., fishnet, mean)%>%
mutate(sloper = ifelse(is.na(sloper), mean(na.omit(sloper)), sloper))
#plot(slopep_net)
steep_net <-
steepp %>%
mutate(steep =1) %>%
aggregate(., fishnet, sum) %>%
mutate(steep = ifelse(is.na(steep), 0, steep))%>%
dplyr::select(steep)
tree_net<-
tree %>%
mutate(tree = 1) %>%
dplyr::select(tree)%>%
aggregate(., fishnet, sum) %>%
mutate(tree = ifelse(is.na(tree), 0, tree))
drain_net<-
drain %>%
mutate(drain =1) %>%
dplyr::select(drain)%>%
aggregate(., fishnet, mean) %>%
mutate(drain = ifelse(is.na(drain), 0, drain))
fac_net <-
facp %>%
aggregate(., fishnet, sum)
ndvi_net<-
ndvip %>%
aggregate(., fishnet, mean)%>%
mutate(layer = ifelse(is.na(layer), mean(na.omit(layer)), layer))
#plot(ndvi_net)
green_net<-
greenp %>%
mutate(green =1) %>%
aggregate(., fishnet, sum) %>%
mutate(green = ifelse(is.na(green), 0, green))%>%
dplyr::select(green)
vars_net <-
cbind(tbc_net,lowdens_net, highdens_net, elevation_net, tree_net, drain_net, fac_net, slopep_net, steep_net, ndvi_net, green_net)%>%
mutate(elevation=GRIDCODE,slope=sloper, ndvi=layer)%>%
dplyr::select(tbc, lowdens, highdens, elevation, tree, drain, fac, slope, steep, ndvi, green)%>%
mutate(uniqueID = rownames(.))
nn_function <- function(measureFrom,measureTo,k) {
measureFrom_Matrix <-
as.matrix(measureFrom)
measureTo_Matrix <-
as.matrix(measureTo)
nn <-
get.knnx(measureTo, measureFrom, k)$nn.dist
output <-
as.data.frame(nn) %>%
rownames_to_column(var = "thisPoint") %>%
gather(points, point_distance, V1:ncol(.)) %>%
arrange(as.numeric(thisPoint)) %>%
group_by(thisPoint) %>%
summarize(pointDistance = mean(point_distance)) %>%
arrange(as.numeric(thisPoint)) %>%
dplyr::select(-thisPoint) %>%
pull()
return(output)
}
vars_net$steep.nn =
nn_function(st_coordinates(st_centroid(vars_net)), st_coordinates(steepp), 3)
vars_net$green.nn =
nn_function(st_coordinates(st_centroid(vars_net)), st_coordinates(greenp), 3)
vars_net$stream.nn =
nn_function(st_coordinates(st_centroid(vars_net)), st_coordinates(streamp), 3)
vars_net2 <-
vars_net %>%
dplyr::select(ndvi,slope,drain,lowdens)
vars_net.long <-
vars_net2 %>%
gather(Variable, value, -geometry)
vars <- unique(vars_net.long$Variable)
mapList <- list()
for(i in vars){
mapList[[i]] <-
ggplot() +
geom_sf(data = filter(vars_net.long, Variable == i), aes(fill=value), colour=NA) +
scale_fill_viridis(name="") +
labs(title=i) +
mapTheme()}
do.call(grid.arrange,c(mapList, ncol =2, top = "Risk Factors by Fishnet"))
river_net <-
river %>%
dplyr::select() %>%
mutate(river = as.numeric(1)) %>%
aggregate(., fishnet, mean) %>%
mutate(river = ifelse(is.na(river), 0, river))
flood_net <-
flood %>%
dplyr::select() %>%
mutate(flood = as.numeric(1)) %>%
aggregate(., fishnet, mean) %>%
mutate(flood = ifelse(is.na(flood), 0, flood))
final_net <-
cbind(vars_net,flood_net)3 Model
We start modeling by building another dataset of just the variables we want to analyze, as well as creating the training and test sets.
#min(final_net$elevation)
#max(final_net$steep)
#max(final_net$green)
final <-
final_net %>%
mutate(elev = elevation-min(elevation)) %>%
mutate(stee = steep/max(steep)) %>%
mutate(gree = green/max(green)) %>%
dplyr::select(-elevation, -steep,-green,-geometry.1)
final$tbc <- as.factor(final$tbc)
final$lowdens <- as.factor(final$lowdens)
final$highdens <- as.factor(final$highdens)
set.seed(3456)
trainIndex <- createDataPartition(y = paste(final$lowdens, final$tbc, final$highdens), p = .70,
list = FALSE,
times = 1)
finalTrain <- final[ trainIndex,]
finalTest <- final[-trainIndex,]Now let’s estimate a logistic regression model.
finalModel <- glm(flood ~ .,
family="binomial"(link="logit"), data = finalTrain %>%
as.data.frame() %>%
dplyr::select(-geometry,-uniqueID))
summary(finalModel)##
## Call:
## glm(formula = flood ~ ., family = binomial(link = "logit"), data = finalTrain %>%
## as.data.frame() %>% dplyr::select(-geometry, -uniqueID))
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.25276 -0.62773 -0.30862 -0.05703 3.12503
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.10957314004 0.44415268135 7.001 0.00000000000254 ***
## tbc1 0.55912734600 0.20323694863 2.751 0.00594 **
## lowdens1 -0.97143956625 0.19577770224 -4.962 0.00000069788235 ***
## highdens1 -0.19559229040 0.13617657327 -1.436 0.15091
## tree 0.09563367812 0.05279868976 1.811 0.07010 .
## drain 0.10886677199 0.11695994208 0.931 0.35196
## fac 0.00000005070 0.00000001017 4.983 0.00000062490279 ***
## slope 0.16533788290 0.10025934101 1.649 0.09913 .
## ndvi -6.19834127423 1.53635152098 -4.034 0.00005472914683 ***
## steep.nn -0.00036079527 0.00012419961 -2.905 0.00367 **
## green.nn -0.00159370720 0.00050833067 -3.135 0.00172 **
## stream.nn -0.00098232882 0.00006885300 -14.267 < 0.0000000000000002 ***
## elev -0.01395127743 0.00144636306 -9.646 < 0.0000000000000002 ***
## stee -0.70335999732 0.65778490606 -1.069 0.28494
## gree 2.88480968562 0.54744825507 5.270 0.00000013675354 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 2782.1 on 2506 degrees of freedom
## Residual deviance: 1933.0 on 2492 degrees of freedom
## AIC: 1963
##
## Number of Fisher Scoring iterations: 6When it’s done, we output the predicted probability of a test set fishnet cell being flooded conditional on our model.
classProbs <- predict(finalModel, finalTest, type="response")
testProbs <- data.frame(obs = as.numeric(finalTest$flood),
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 flooded","Flood"),
name = "")
testProbs$predClass = ifelse(testProbs$pred > .5 ,1,0)
caret::confusionMatrix(reference = as.factor(testProbs$obs),
data = as.factor(testProbs$predClass),
positive = "1")## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 769 114
## 1 31 154
##
## Accuracy : 0.8642
## 95% CI : (0.8422, 0.8842)
## No Information Rate : 0.7491
## P-Value [Acc > NIR] : < 0.00000000000000022
##
## Kappa : 0.5974
##
## Mcnemar's Test P-Value : 0.000000000009778
##
## Sensitivity : 0.5746
## Specificity : 0.9613
## Pos Pred Value : 0.8324
## Neg Pred Value : 0.8709
## Prevalence : 0.2509
## Detection Rate : 0.1442
## Detection Prevalence : 0.1732
## Balanced Accuracy : 0.7679
##
## 'Positive' Class : 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') 
auc(testProbs$obs, testProbs$pred)## Area under the curve: 0.85We examine the goodness of fit metrics for random k-fold cross validation.
ctrl <- trainControl(method = "cv",
number = 100,
savePredictions = TRUE)
cvFit <- train(as.factor(flood) ~ ., data = finalTrain %>%
as.data.frame() %>%
dplyr::select(-geometry,-uniqueID),
method="glm", family="binomial",
trControl = ctrl)
cvFit## Generalized Linear Model
##
## 2507 samples
## 14 predictor
## 2 classes: '0', '1'
##
## No pre-processing
## Resampling: Cross-Validated (100 fold)
## Summary of sample sizes: 2481, 2482, 2482, 2482, 2481, 2482, ...
## Resampling results:
##
## Accuracy Kappa
## 0.8469077 0.5296552ggplot(as.data.frame(cvFit$resample), aes(Accuracy)) +
geom_histogram() +
scale_x_continuous(limits = c(0, 1)) +
labs(x="Accuracy",
y="Count")
4 Map predictions
Now that we have tuned out model, let’s predict for the entire dataset.
allPredictions <-
predict(cvFit, final, type="prob")[,2]
finalflood <-
cbind(final,allPredictions) %>%
mutate(allPredictions = round(allPredictions * 100))
ggplot() +
geom_sf(data=finalflood, aes(fill=factor(ntile(allPredictions,5))), colour=NA) +
scale_fill_manual(values = c("#edf8fb","#b3cde3","#8c96c6","#8856a7","#810f7c"),
labels=as.character(quantile(finalflood$allPredictions,
c(0.1,.2,.4,.6,.8),na.rm=T)),
name="Predicted\nProbabilities(%)\n(Quintile\nBreaks)") +
mapTheme() +
labs(title="")
ggplot() +
geom_sf(data=finalflood, aes(fill=factor(ntile(allPredictions,5))), colour=NA) +
scale_fill_manual(values = c("#edf8fb","#b3cde3","#8c96c6","#8856a7","#810f7c"),
labels=as.character(quantile(finalflood$allPredictions,
c(0.1,.2,.4,.6,.8),na.rm=T)),
name="Predicted\nProbabilities(%)\n(Quintile\nBreaks)") +
geom_sf(data=finalflood %>%
filter(flood == 1),
fill="red",colour=NA) +
mapTheme() +
labs(title="Predicted Probabilities and Past Flood",
subtitle="Calgary; Past Flood in red")
To further examine the spatial generalizability of our model, we plot the confusion matrix over space (for the train set).
trainPredictions <-
predict(cvFit, finalTrain, type="prob")[,2]
matrix <-
cbind(finalTrain,trainPredictions) %>%
mutate(trainPredictions = round(trainPredictions * 100)) %>%
mutate(pre = ifelse(trainPredictions>22,1,0)) %>%
mutate(compare=
ifelse(pre==0 & flood==0,"True_Negative",
ifelse(pre==1 & flood==1,"True_Positive",
ifelse(pre==0 & flood==1,"False_Negative",
ifelse(pre==1 & flood==0,"False_Positive",0)))))
ggplot() +
geom_sf(data=matrix, aes(fill=factor(compare)), colour=NA) +
scale_fill_manual(values = c("red","lightblue","grey","blue"),
labels= c("False_Negative","False_Positive","True_Negative","True_Positive"),
name="Prediction\n& Reference\n(Train Set)") +
mapTheme() 
Lastly, we present the final inundation map for Calgary (entire dataset).
finalpre <-
finalflood %>%
mutate(pre = ifelse(allPredictions<22,1,0))
ggplot() +
geom_sf(data=finalpre, aes(fill=factor(pre)), colour=NA) +
scale_fill_manual(values = c("blue","grey"),
labels=c("flood","non-flood"),
name="Final Prediction\n(entire dataset)") +
mapTheme() 
5 Comparable City: Denver
After training and validating our model on Calgary, we apply it on a comparable city: Denver.
Download data and create base maps, including city boundary, flood map and fishnet.
limit <-
st_read("./denver/limit.shp") %>%
st_transform(crs=6428)
fishnet <-
st_make_grid(limit, cellsize = 500) %>%
st_sf()
fishnet <-
fishnet[limit,] %>%
mutate(uniqueID = rownames(.)) %>%
dplyr::select(uniqueID)
#plot(fishnet)Next, we gather, wrangle and plot the flood risk factor data.
#library(raster)
#plot(landuse)
landuse <- st_read("./denver/existing_landuse_2018/existing_landuse_2018.shp") %>%
st_transform(crs=6428)
elevation <- st_read("./denver/elevation.shp") %>%
st_transform(crs=6428)
drain <- st_read("./denver/hydro/drainv.shp") %>%
st_transform(crs=6428)
tree <-
st_read("./denver/tree_inventory/tree_inventory.shp") %>%
st_transform(crs=6428)
sloper <- raster("./denver/sloper.tif")
slopep <-
rasterToPoints(sloper) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y"), crs = crs(sloper))%>%
st_transform(st_crs(limit))
steepp <- slopep %>%
filter(sloper > 3)
fac <- raster("./denver/hydro/fac.tif")
facp <-
rasterToPoints(fac) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y"), crs = st_crs(fac))%>%
st_transform(st_crs(limit))
band4 <- raster("./denver/data/band4.tif")
band5 <- raster("./denver/data/band5.tif")
image <- stack(band4, band5)
ndvi <- (image[[2]] - image[[1]])/(image[[2]] + image[[1]])
#plot(ndvi)
ndvip <-
rasterToPoints(ndvi) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y"), crs = st_crs(ndvi))%>%
st_transform(st_crs(limit))
greenp <- ndvip %>%
filter(layer > 0.35) %>%
st_transform(st_crs(limit))
streamp <- ndvip %>%
filter(layer < 0.01) %>%
st_transform(st_crs(limit))We then join the flood map and the risk factors to the fishnet.
#unique(landuse$CPD_LANDUS)
tbc_net <-
landuse %>%
filter(CPD_LANDUS== "Vacant" | CPD_LANDUS== "Agriculture"| CPD_LANDUS== "Park/Open Space") %>%
mutate(tbc =1) %>%
aggregate(., fishnet, mean) %>%
mutate(tbc = ifelse(is.na(tbc), 0, tbc))%>%
dplyr::select(tbc)
lowdens_net <-
landuse %>%
filter(CPD_LANDUS== "Single-unit Residential"|CPD_LANDUS== "Two-unit Residential"|CPD_LANDUS== "Industrial"|CPD_LANDUS== "Trans/Comm/Utilities"|CPD_LANDUS== "ROW/Road"|CPD_LANDUS== "Parking") %>%
mutate(lowdens =1) %>%
aggregate(., fishnet, mean) %>%
mutate(lowdens = ifelse(is.na(lowdens), 0, lowdens))%>%
dplyr::select(lowdens)
highdens_net <-
landuse %>%
filter(CPD_LANDUS== "Multi-unit Residential"|CPD_LANDUS== "Office"|CPD_LANDUS== "Mixed-use"|CPD_LANDUS== "Entertainment/Cultural"|CPD_LANDUS== "Commercial/Retail") %>%
mutate(highdens =1) %>%
aggregate(., fishnet, mean) %>%
mutate(highdens = ifelse(is.na(highdens), 0, highdens))%>%
dplyr::select(highdens)
elevation_net <-
elevation %>%
aggregate(., fishnet, mean)%>%
mutate(GRIDCODE = ifelse(is.na(GRIDCODE), mean(na.omit(GRIDCODE)), GRIDCODE))
#plot(elevation_net)
slopep_net <-
slopep %>%
aggregate(., fishnet, mean)%>%
mutate(sloper = ifelse(is.na(sloper), mean(na.omit(sloper)), sloper))
#plot(slopep_net)
steep_net <-
steepp %>%
mutate(steep =1) %>%
aggregate(., fishnet, sum) %>%
mutate(steep = ifelse(is.na(steep), 0, steep))%>%
dplyr::select(steep)
tree_net<-
tree %>%
mutate(tree = 1) %>%
dplyr::select(tree)%>%
aggregate(., fishnet, sum) %>%
mutate(tree = ifelse(is.na(tree), 0, tree))
drain_net<-
drain %>%
mutate(drain =1) %>%
dplyr::select(drain)%>%
aggregate(., fishnet, mean) %>%
mutate(drain = ifelse(is.na(drain), 0, drain))
fac_net <-
facp %>%
aggregate(., fishnet, sum)
ndvi_net<-
ndvip %>%
aggregate(., fishnet, mean)%>%
mutate(layer = ifelse(is.na(layer), mean(na.omit(layer)), layer))
#plot(ndvi_net)
green_net<-
greenp %>%
mutate(green =1) %>%
aggregate(., fishnet, sum) %>%
mutate(green = ifelse(is.na(green), 0, green))%>%
dplyr::select(green)
vars_net <-
cbind(tbc_net,lowdens_net, highdens_net, elevation_net, tree_net, drain_net, fac_net, slopep_net, steep_net, ndvi_net, green_net)%>%
mutate(elevation=GRIDCODE,slope=sloper, ndvi=layer)%>%
dplyr::select(tbc, lowdens, highdens, elevation, tree, drain, fac, slope, steep, ndvi, green)%>%
mutate(uniqueID = rownames(.))
nn_function <- function(measureFrom,measureTo,k) {
measureFrom_Matrix <-
as.matrix(measureFrom)
measureTo_Matrix <-
as.matrix(measureTo)
nn <-
get.knnx(measureTo, measureFrom, k)$nn.dist
output <-
as.data.frame(nn) %>%
rownames_to_column(var = "thisPoint") %>%
gather(points, point_distance, V1:ncol(.)) %>%
arrange(as.numeric(thisPoint)) %>%
group_by(thisPoint) %>%
summarize(pointDistance = mean(point_distance)) %>%
arrange(as.numeric(thisPoint)) %>%
dplyr::select(-thisPoint) %>%
pull()
return(output)
}
vars_net$steep.nn =
nn_function(st_coordinates(st_centroid(vars_net)), st_coordinates(steepp), 3)
vars_net$green.nn =
nn_function(st_coordinates(st_centroid(vars_net)), st_coordinates(greenp), 3)
vars_net$stream.nn =
nn_function(st_coordinates(st_centroid(vars_net)), st_coordinates(streamp), 3)
vars_net2 <-
vars_net %>%
dplyr::select(ndvi,slope,drain,lowdens)
vars_net.long <-
vars_net2 %>%
gather(Variable, value, -geometry)
vars <- unique(vars_net.long$Variable)
mapList <- list()
for(i in vars){
mapList[[i]] <-
ggplot() +
geom_sf(data = filter(vars_net.long, Variable == i), aes(fill=value), colour=NA) +
scale_fill_viridis(name="") +
labs(title=i) +
mapTheme()}
do.call(grid.arrange,c(mapList, ncol =2, top = "Risk Factors by Fishnet"))
Lastly, we present the final inundation map for Denver (entire dataset).
#min(final_net$elevation)
#max(final_net$steep)
#max(final_net$green)
final2 <-
vars_net %>%
mutate(elev = elevation-min(elevation)) %>%
mutate(stee = steep/max(steep)) %>%
mutate(gree = green/max(green)) %>%
dplyr::select(-elevation, -steep,-green)
final2$tbc <- as.factor(final2$tbc)
final2$lowdens <- as.factor(final2$lowdens)
final2$highdens <- as.factor(final2$highdens)
classProbs2 <- predict(finalModel, final2, type="response")
#typeof(classProbs2)
finalflood2 <-
cbind(final2,classProbs2) %>%
mutate(allPredictions2 = round(classProbs2 * 100))
ggplot() +
geom_sf(data=finalflood2, aes(fill=factor(ntile(allPredictions2,5))), colour=NA) +
scale_fill_manual(values = c("#edf8fb","#b3cde3","#8c96c6","#8856a7","#810f7c"),
labels=as.character(quantile(finalflood2$allPredictions2,
c(0.1,.2,.4,.6,.8),na.rm=T)),
name="Predicted\nProbabilities(%)\n(Quintile\nBreaks)") +
mapTheme() +
labs(title="")
finalpre2 <-
finalflood2 %>%
mutate(pre = ifelse(allPredictions2<22,1,0))
ggplot() +
geom_sf(data=finalpre2, aes(fill=factor(pre)), colour=NA) +
scale_fill_manual(values = c("blue","grey"),
labels=c("flood","non-flood"),
name="Final Prediction\n(comparable city)") +
mapTheme() 