index.Rmd

---
title: "Interpreting Stormwater management model (SWMM) Output"
author: "Seyed Hamed Ghodsi (seyedham@buffalo.edu)"
date: 'Deadline: November 9, 2018'
output: html_document
---

***
# Inroduction

This project provides documentation for the analysis of an anonymous city sewer collection system. The purpose of this report is to provide the interpretation of **SWMM** output files.  
EPA's Stormwater Management Model (SWMM) is used for single event or long-term simulations of water runoff quantity and quality in primarily urban areas–although there are also many applications that can be used for drainage systems in non-urban areas. It is used throughout the world for planning, analysis, and design related to stormwater runoff, combined and sanitary sewers, and other drainage systems.  
The SWMM output would be made in two files. The first one is a report file (*".rpt"*) which includes all the cumulative results from the simulation model such as runoff and outfall loading summary. The second one is a binary file (*".out"*) which include the time series analysis of different elements of the model like subcatchments, nodes and links. Usually researchers use SWMM software interface to analyze the time series result, which is stored in the binary file. However, in this report the analysis would be done using `R`. As a result, the simulation model results would be analyzed and showed in two separated parts. It should be noted that to get the output files, the input file was run manually through the SWMM model interface.

***  

# Data Analysis  

## Report Data  


The report file which would be used as a data set in this project, has been uploaded in [github](https://raw.githubusercontent.com/AdamWilsonLabEDU/geo503-2018-finalproject-HamedGhodsi90/master/data/SWMMReportFile.rpt). This file shows the SWMM output as a report file. 
```{r echo=FALSE, message=FALSE, warning=FALSE, paged.print=FALSE}

#install.packages(c("dplyr", "swmmr", "devtools", "knitr","purrr", "gridExtra","ggplot2"))

library(dplyr)
library(swmmr)
library(devtools)
library(purrr)
library(knitr)
library(gridExtra)
library(grid)
library(ggplot2)
library(readr)
```

``` {r echo = FALSE}
#.......................Needed Functions ....................
percent_Function <- function(xp, digits = 2, format = "f", ...) {
  paste0(formatC(100 * xp, format = format, digits = digits, ...), "%")
}


tt3 <- ttheme_minimal(
  core=list(bg_params = list(fill = blues9[1:2], col=NA),
            fg_params=list(fontface=3)),
  colhead=list(fg_params=list(col="navyblue", fontface=4L)),
  rowhead=list(fg_params=list(col="orange", fontface=3L)))

```

With the following lines of `R` code, the input file could be readable. The first line, would download the data set from [github](https://raw.githubusercontent.com/AdamWilsonLabEDU/geo503-2018-finalproject-HamedGhodsi90/master/data/SWMMReportFile.rpt). I made this line as a comment. Because just for the first time it is needed to be downloaded. The rest of the time, the file is in the current folder and can be read easily. Therefore, it is not needed to download it again. For the users who want to run this model for the first time, they should run this line (remove the #) 

```{r error=FALSE, message=FALSE, warning=FALSE}

download.file("https://raw.githubusercontent.com/AdamWilsonLabEDU/geo503-2018-finalproject-HamedGhodsi90/master/data/SWMMReportFile.rpt", "Report_file.rpt")

report_file <- read_rpt ("Report_file.rpt")

```

The report file has different sections which are:

```{r echo=FALSE, error=FALSE, message=FALSE, warning=FALSE}

grid.table(matrix(names(report_file[1:15]), nrow = 5), theme=tt3) 

```

The two most important sections, which have to be analyzed in the stormwater management process, are subcatchments and outfalls. 

### Subcatchment analysis:

``` {r echo=FALSE}  
# 3.1 Subcatchment Summary preparation ...................................

rpt_subcatchment_raw <- report_file $ subcatchment_summary                 # To load the subcatchment summary part

num_cor  <- which(rpt_subcatchment_raw $ Name=="*******************")                             # To find a row number befor the LID part
num_cor2 <- which(rpt_subcatchment_raw $ Rain_Gage != rpt_subcatchment_raw$Rain_Gage [1]) # To find a row number for the bug

rpt_subcatchment_pre <- rpt_subcatchment_raw %>% filter(between(row_number(), 1, num_cor[1]-1))   # To remove the LID part from the end of the list

rpt_subcatchment_pre [num_cor2[1],"Outlet"]  = rpt_subcatchment_pre [num_cor2[1],"Rain_Gage"];            # To fix a bug in the data (row number = num_cor2)
rpt_subcatchment_pre [num_cor2[1],"Rain_Gage"]  = rpt_subcatchment_pre [num_cor2[1],"Perc_Slope"];
rpt_subcatchment_pre [num_cor2[1],"Perc_Slope"]  = rpt_subcatchment_pre [num_cor2[1],"Perc_Imperv"];
rpt_subcatchment_pre [num_cor2[1],"Perc_Imperv"]  = rpt_subcatchment_pre [num_cor2[1],"Width"];
rpt_subcatchment_pre [num_cor2[1],"Width"]  = substr(rpt_subcatchment_pre $ Area[num_cor2[1]], 9,20)
rpt_subcatchment_pre [num_cor2[1],"Area"]  = substr(rpt_subcatchment_pre $ Area[num_cor2[1]], 0,8)

rpt_subcatchment_pre $ Area <- as.numeric(as.character(rpt_subcatchment_pre $ Area))                      # To change the character to number
rpt_subcatchment_pre $ Width <- as.numeric(as.character(rpt_subcatchment_pre $ Width))
rpt_subcatchment_pre $ Perc_Imperv <- as.numeric(as.character(rpt_subcatchment_pre $ Perc_Imperv))
rpt_subcatchment_pre $ Perc_Slope <- as.numeric(as.character(rpt_subcatchment_pre $ Perc_Slope))

rpt_subcatchment <- rpt_subcatchment_pre  %>% arrange(desc(Area))
# View (rpt_subcatchment_pre)
# View (rpt_subcatchment)


#3.2: Subcatchment Runoff prepareation ...................................

rpt_Runoff_subcatchment_raw <- report_file $ subcatchment_runoff_summary 

num_cor3  <- which(rpt_Runoff_subcatchment_raw $ Subcatchment == "***********************") 

rpt_Runoff_subcatchment_pre <- report_file $ subcatchment_runoff_summary %>% filter(between(row_number(), 1, 59655))

rpt_Runoff_subcatchment_pre $ Total_Precip <- as.numeric(as.character(rpt_Runoff_subcatchment_pre $ Total_Precip))
rpt_Runoff_subcatchment_pre $ Total_Runon <- as.numeric(as.character(rpt_Runoff_subcatchment_pre $ Total_Runon))
rpt_Runoff_subcatchment_pre $ Total_Evap <- as.numeric(as.character(rpt_Runoff_subcatchment_pre $ Total_Evap))
rpt_Runoff_subcatchment_pre $ Total_Infil <- as.numeric(as.character(rpt_Runoff_subcatchment_pre$ Total_Infil))
rpt_Runoff_subcatchment_pre $ Total_Runoff_Depth <- as.numeric(as.character(rpt_Runoff_subcatchment_pre $ Total_Runoff_Depth))
rpt_Runoff_subcatchment_pre $ Total_Runoff_Volume <- as.numeric(as.character(rpt_Runoff_subcatchment_pre $ Total_Runoff_Volume))
rpt_Runoff_subcatchment_pre $ Total_Peak_Runoff <- as.numeric(as.character(rpt_Runoff_subcatchment_pre $ Total_Peak_Runoff))
rpt_Runoff_subcatchment_pre $ Total_Runoff_Coeff <- as.numeric(as.character(rpt_Runoff_subcatchment_pre $ Total_Runoff_Coeff))

rpt_Runoff_subcatchment <- rpt_Runoff_subcatchment_pre %>% arrange (desc(Total_Runoff_Volume))
# View(rpt_Runoff_subcatchment_raw)
# View (rpt_Runoff_subcatchment)


#3.3: Subcatchment Summarize............................................................................
rpt_subcatchment_summary <- rpt_subcatchment %>% 
  rename (Subcatchment=Name) %>% 
    select(Subcatchment, Area, Width, Perc_Imperv, Outlet) %>%
      inner_join(rpt_Runoff_subcatchment, by=c("Subcatchment")) %>% 
        mutate(
             Runoff_Volume_perc = percent_Function(Total_Runoff_Volume/sum(Total_Runoff_Volume, na.rm=TRUE)), 
             Runoff_Depth_perc = percent_Function(Total_Runoff_Depth/sum(Total_Runoff_Depth)),
             Area_perc = percent_Function(Area/sum(Area)),
             Rank_Area = dense_rank(desc(Area)),
             Rank_Dep = dense_rank(desc(Total_Runoff_Depth)),
             Rank_Vol = dense_rank(desc(Total_Runoff_Volume))
             ) %>%
          select(
                Subcatchment, 
                Area, Area_perc, Rank_Area,
                Total_Runoff_Depth,Runoff_Depth_perc, Rank_Dep,
                Total_Runoff_Volume,Runoff_Volume_perc, Rank_Vol
                ) %>%
            arrange(desc(Total_Runoff_Volume))
# View (rpt_subcatchment_summary)


# 3.3.1: Sorting subcatchments based on total runoff volume .........................

rpt_sub_Vol_Area <- rpt_subcatchment_summary %>% 
  select(Subcatchment,Total_Runoff_Volume, Runoff_Volume_perc, Area, Rank_Area) %>%
    rename(
           "Runoff Volume (MG)"= Total_Runoff_Volume,
           "%" = Runoff_Volume_perc,
           "Area (hec)" = Area,
           "Rank Area" = Rank_Area
           )

rpt_sub_Vol_Area_top <- rpt_sub_Vol_Area %>%
  filter(between(row_number(), 1, 15))
# View(rpt_sub_Vol_Area_top)


# 3.3.2: Sorting subcatchments based on total runoff depth .........................

rpt_sub_Dep_Area <- rpt_subcatchment_summary  %>% 
  arrange(desc(Total_Runoff_Depth)) %>% 
    select(Subcatchment, Total_Runoff_Depth, Runoff_Depth_perc, Area, Rank_Area) %>%
      rename (
             "Runoff Depth (inch)"= Total_Runoff_Depth,
              "%" = Runoff_Depth_perc,
              "Area (hec)" = Area,
              "Rank Area" = Rank_Area
              ) 

rpt_sub_Dep_Area_top <- rpt_sub_Dep_Area %>% 
  filter(between(row_number(), 1, 15))
# View (rpt_sub_Dep_Area_top)


# 3.3.3: Sorting subcatchments based on Area ........... .........................

rpt_sub_Area_Vol <- rpt_subcatchment_summary  %>% 
  arrange(desc(Area)) %>%
    select(Subcatchment, Area, Area_perc, Total_Runoff_Volume, Runoff_Volume_perc) %>%
      rename(
             "Area (hec)" = Area,
             "Area %" = Area_perc,
             "Runoff Volume (MG)"= Total_Runoff_Volume,
             "Vol %" = Runoff_Volume_perc
             )

rpt_sub_Area_Vol_top <- rpt_sub_Area_Vol %>%
  filter(between(row_number(), 1, 15))
# View (rpt_sub_Area_Vol_top)

# The cumulative runoff ...................................................................
rpt_subcatchment_summarize <- summarise (rpt_subcatchment_summary, 
                                         Sum_Runoff_Volume_MG = sum(Total_Runoff_Volume), 
                                         Sum_Runoff_Depth_inch = sum(Total_Runoff_Depth)
                                         )

```

Number of subcatchments in the watershed: **`r nrow(rpt_Runoff_subcatchment)`**    
The subcatchments have been arranged based on total runoff volume, total runoff depth, and area which the top ones are shown in the following tables.  

``` {r echo = FALSE} 
# To plot the tables ..............................................
# grid.arrange(tableGrob(rpt_sub_top_Vol_Area, theme=tt3),
#              tableGrob(rpt_sub_top_Dep, theme=tt3))
grid.table(rpt_sub_Vol_Area_top, theme=tt3)  
```
``` {r echo = FALSE} 
grid.table(rpt_sub_Dep_Area_top, theme=tt3)  
```
``` {r echo = FALSE}
grid.table(rpt_sub_Area_Vol_top, theme=tt3)
```

Knowing the subcatchments with the high value of total runoff and total depth would help the managers to understand the watershed in a better way. For example in order to decrease the runoff volume, low impact developments (**LID**) could be utilized. TO implement the LIDs, these selected subcatchments could have more priority.  
``` {r echo=FALSE}

# 3.4: Upstream subcatchments ..........................................................

rpt_sub_ups <- rpt_subcatchment_summary %>% 
  select(Subcatchment, Area) %>%
    filter (Area > 2000) %>% arrange(desc(Area))
rpt_sub_ups_name <- matrix()
for (i in 1:nrow(rpt_sub_ups))
  rpt_sub_ups_name[i] <- rpt_sub_ups[[i,1]]
```

by looking at the the subcatchments's area, it could be realized that The first **`r nrow(rpt_sub_ups)`** subcatchments' area are much higher than the rest ones. These subcatchments are the watershed upsream areas, which act as boudary conditions. These subcatchments are:    
**`r rpt_sub_ups_name[[1]]`,  `r rpt_sub_ups_name[[2]]`,  `r rpt_sub_ups_name[[3]]`,  `r rpt_sub_ups_name[[4]]`, and `r rpt_sub_ups_name[[5]]`**. 

Each subcatchment was named based on a critera which are:  
``` {r echo=FALSE}
# 3.5: subcatchments name Definition and number ..........................................

rpt_sub_subname = matrix()
rpt_sub_subname <- c("Cheektowaga5", "ECSD1", "ECSD4", "WSSD13", 
                     "DetBasin", "PlantInf", "RainGardens", 
                     "Buf1", "Buf2", "GLawn", "HDC", "HLawn",
                     "HSplash", "HStreet", "Hstreet", "Lat",
                     "Main", "Parking", "RCom", "SAimp", "SAperv")

rpt_sub_run <- rpt_subcatchment_summary %>% arrange(desc(Area)) %>%
  select(Subcatchment)
# View(rpt_sub_run)

rpt_sub_subname_num = matrix()
for (i in 1: length(rpt_sub_subname)) {
  rpt_sub_subname_num [i] <- rpt_sub_run %>% 
          filter(grepl(rpt_sub_subname[i] ,Subcatchment)) %>% nrow
}     

rpt_sub_subname_num_to <- cbind(rpt_sub_subname, rpt_sub_subname_num)
```

`Cheektowaga5`: inflow from outside sewer district (OSD) (**`r rpt_sub_subname_num_to[1,2]`** subcatchment)     
`ECSD1`: inflow from outside sewer district (OSD) (**`r rpt_sub_subname_num_to[2,2]`** subcatchment)     
`ECSD4`: inflow from outside sewer district (OSD) (**`r rpt_sub_subname_num_to[3,2]`** subcatchment)       
`WSSD13`: inflow from outside sewer district (OSD) (**`r rpt_sub_subname_num_to[4,2]`** subcatchment)         
`DetBasin`: detention basin (**`r rpt_sub_subname_num_to[5,2]`** subcatchment)        
`PlantInf`: green infrastructure project (infiltration planters) (**`r rpt_sub_subname_num_to[6,2]`** subcatchments)     
`RainGardens`: green infrastructure project (rain gardens) (**`r rpt_sub_subname_num_to[7,2]`** subcatchments)    
`Buf1`: House Buffer Area (type 1, disconnected downspouts) (**`r rpt_sub_subname_num_to[8,2]`** subcatchments)    
`Buf2`: House Buffer Area (type 2, all other areas) (**`r rpt_sub_subname_num_to[9,2]`** subcatchments)  
`GLawn`: Lawn (**`r rpt_sub_subname_num_to[10,2]`** subcatchments)     
`HDC`: Roofs DC (**`r rpt_sub_subname_num_to[11,2]`** subcatchments)       
`HLawn`: Roofs to Lawn (**`r rpt_sub_subname_num_to[12,2]`** subcatchments)       
`HSplash`: Roofs Splashing (impervious component) and House Buffer Area Splashing (pervious component) (**`r rpt_sub_subname_num_to[13,2]`** subcatchments)    
`HStreet`: Roofs to Street (**`r rpt_sub_subname_num_to[14,2]`** subcatchments)     
`Hstreet`: Street (**`r rpt_sub_subname_num_to[15,2]`** subcatchments)    
`Lat`: Lateral Area  (**`r rpt_sub_subname_num_to[16,2]`** subcatchments)    
`Main`:Mains Area (**`r rpt_sub_subname_num_to[17,2]`** subcatchments)     
`Parking`: Commercial Parking Lots (**`r rpt_sub_subname_num_to[18,2]`** subcatchments)      
`RCom`: Roofs Commercial (**`r rpt_sub_subname_num_to[19,2]`** subcatchments)     
`SAimp`: impervious subarea (**`r rpt_sub_subname_num_to[20,2]`** subcatchments)   
`SAperv`: pervious subarea (**`r rpt_sub_subname_num_to[21,2]`** subcatchments)




### Outfall Analysis


























## Binary Data










***

# Conclusions




***

# References















<!-- # Materials and methods -->

<!-- Narrative and most code will go here.  Describe what you are doing and show how to do it (with code). -->

<!-- You can do bullets like this: -->

<!-- * The first most important thing -->
<!-- * The second most important thing -->
<!-- * The third most important thing -->

<!-- You can do numbers like this: -->

<!-- 1. The first most important thing -->
<!-- 2. The second most important thing -->
<!-- 3. The third most important thing -->

<!-- See [http://rmarkdown.rstudio.com/](http://rmarkdown.rstudio.com/) for all the amazing things you can do. -->


<!-- Here's my first code chunk. -->
<!-- ```{r} -->
<!-- 1+2 -->
<!-- ``` -->

<!-- Load any required packages in a code chunk (you may need to install some packages): -->

<!-- ```{r, message=F, warning=F} -->
<!-- library(dplyr) -->
<!-- library(ggplot2) -->
<!-- library(maps) -->
<!-- library(spocc) -->
<!-- knitr::opts_chunk$set(cache=TRUE)  # cache the results for quick compiling -->
<!-- ``` -->


<!-- We are going to use the `occ()` function to download _occurrence_ records for the American robin (_Turdus migratorius_) from the [Global Biodiversity Information Facility](gbif.org). -->

<!-- <img src="https://upload.wikimedia.org/wikipedia/commons/b/b8/Turdus-migratorius-002.jpg" alt="alt text" width="200"> -->
<!-- <small><small><small>Licensed under CC BY-SA 3.0 via [Wikimedia Commons](https://commons.wikimedia.org/wiki/File:Turdus-migratorius-002.jpg#/media/File:Turdus-migratorius-002.jpg)</small></small></small> -->

<!-- ```{r, warning=F}  -->
<!-- ## define which species to query -->
<!-- sp='Turdus migratorius' -->

<!-- ## run the query and convert to data.frame() -->
<!-- d = occ(query=sp, from='ebird',limit = 100) %>% occ2df() -->
<!-- ``` -->
<!-- This can take a few seconds. -->

<!-- ```{r, fig.width=6, fig.height=3, fig.cap="Map illustrating the known occurrence locations"} -->
<!-- # Load coastline -->
<!-- map=map_data("world") -->

<!-- ggplot(d,aes(x=longitude,y=latitude))+ -->
<!--   geom_polygon(aes(x=long,y=lat,group=group),data=map)+ -->
<!--   geom_point(col="red")+ -->
<!--   coord_equal() -->
<!-- ``` -->


<!-- # Results -->

<!-- Show tables, plots, etc. and describe them. -->

<!-- ```{r} -->
<!-- hist(d$latitude) -->
<!-- ``` -->

<!-- # Conclusions -->

<!-- What have you learned?  Are there any broader implications? -->

<!-- # References -->