UtahMuncipalIndustrialData-WaterPlanComment.Rmd

---
title: "Comments on 2021 draft Utah state water plan"
author: "David E. Rosenberg"
date: "11/12/2021"
output: pdf_document
---

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

## Introduction

This is an R Markdown document. This document plots Municipal and Industrial water use data downloaded from the Utah Division of Water Resources (UDWR) data portal. The plot is made in support of 2 comments on the draft 2021 Utah state water plan. Below, findings are listed and the method described. A final section shows the resultant plot.

## Findings

1. Make the state conservation goals more agressive. A 20% reducation from current use will only get 60 Utah providers to reach the 191 gallon per capita per day level that Salt Lake City currently uses. Another 267 Utah providers will still use much more than Salt Lake City.

2. Target conservation efforts to the Utah customers and Utah water providers that use the most water. For example, Garden City, Wellsville, Bringham City, Park City, and Green River all use more than double the amount of Salt Lake City on a per-capita basis. Cities in Washington County like Hurricane and St. George also use more water on a per capita basis than Salt Lake City. There is potential to mount aggressive conservation campaigns with these cities so their water use comes down to the level of the lowest 141 providers across the state.

## Analysis Steps

1. Download data for 2015 to 2018 from the data portal,  https://dwre-utahdnr.opendata.arcgis.com/pages/municipal-and-industrial-data.
2. Read yearly data into R and combine into a single data frame.
3. Sorted the 2018 data by the TotalGPCD field (Total gallons per capita per day)
4. Plot as a bar graph with providers ranked from largest to smallest GPCD.
5. Break the plot by the Utah average of 240 gallons per person per day and twice that volume (480 gpcd).
6. Annotate the bar plot with the number of providers in each block.

## Requested Citation

David E. Rosenberg (2021). "Comments on draft 2021 Utah State Water Plan."

\newpage
## Read in data (year)
```{r LoadData, echo=FALSE, warning=FALSE, message=FALSE}
# UtahMunicipalIndustiralUse_InitAnalysis.r
#
# Initial analysis of Utah's Municipal and Industrial Water Use Data
# 2018-2015
#
# Downloaded from https://dwre-utahdnr.opendata.arcgis.com/pages/municipal-and-industrial-data
#
# David E. Rosenberg
# November 28, 2020
# Updated November 12, 2021
# Utah State University
# david.rosenberg@usu.edu

rm(list = ls())  #Clear history

# Load required libraies

if (!require(tidyverse)) { 
  install.packages("tidyverse", repos="http://cran.r-project.org") 
  library(tidyverse) 
}

if (!require(readxl)) { 
  install.packages("readxl", repos="http://cran.r-project.org") 
  library(readxl) 
}

  
if (!require(RColorBrewer)) { 
  install.packages("RColorBrewer",repos="http://cran.r-project.org") 
  library(RColorBrewer) # 
}

if (!require(dplyr)) { 
  install.packages("dplyr",repos="http://cran.r-project.org") 
  library(dplyr) # 
}

if (!require(expss)) { 
  install.packages("expss",repos="http://cran.r-project.org") 
  library(expss) # 
}

if (!require(reshape)) { 
  install.packages("reshape", repos="http://cran.r-project.org") 
  library(reshape) 
}

if (!require(reshape2)) { 
  install.packages("reshape2", repos="http://cran.r-project.org") 
  library(reshape2) 
}


if (!require(pracma)) { 
  install.packages("pracma", repos="http://cran.r-project.org") 
  library(pracma) 
}

if (!require(lubridate)) { 
  install.packages("lubridate", repos="http://cran.r-project.org") 
  library(lubridate) 
}

if (!require(directlabels)) { 
  install.packages("directlabels", repo="http://cran.r-project.org")
  library(directlabels) 
}


#if (!require(plyr)) { 
#  install.packages("plyr", repo="http://cran.r-project.org")
#  library(plyr) 
#}

if (!require(ggrepel)) { 
  devtools::install_github("slowkow/ggrepel")
  library(ggrepel) 
}

#library(dygraphs)
#library(xts)          # To make the convertion data-frame / xts format
library(tidyverse)
library(lubridate)


### 0. Definitions

sFileBase <- "_Municipal_and_Industrial_Water_Use_Databases.csv"

#Define the River locations to use
dfFileYears <- data.frame(year = seq(2015,2018, by=1))
dfFileYears$filename <- paste0(dfFileYears$year,sFileBase,"")

### 1. Read IN the data files. Rbind so long list

for(i in 1:nrow(dfFileYears)){
 # i <- 1
  
  print(dfFileYears$year[i])
  
  # Read in the historical Powell data
  dfTempFile <- read.csv(file=dfFileYears$filename[i], 
                         header=TRUE, 
                         
                         stringsAsFactors=FALSE,
                         sep=",")
  
  dfTempFile$year <- dfFileYears$year[i]                     #factor(dfLocations$rivermile[i], levels = dfLocations$rivermile)
  #Rename the first column because they differ amoung files
  cColNames <- colnames(dfTempFile)
  cColNames[1] <- "FID"
  colnames(dfTempFile) <- cColNames
  
  if(dfTempFile$year[i] %in% c(2018)) { #We need to adjust some other columns
    cColNamesTemp <- cColNames
    cColNamesTemp[c(12, 25,26,27,28)] <- c("TotPotaGPCD", "TotalSup", "WellsSup","SpringsSup","SurfaceSup")
    
    colnames(dfTempFile) <- cColNamesTemp
    
    #Add blank columns for these fields
    cColsAdd <- c("ComSecGPCD", "IndSecGPCD", "InsSecGPCD", "InsSecoUse", "ResSecGPCD", "ResSecoUse", "ComSecoUse", "IndSecoUse", "IndSecoUse")
    for (i in (1:length(cColsAdd))) {
      dfTempFile[[cColsAdd[i]]] = 0
      }
    
  }
  
  if(dfTempFile$year[i] %in% c(2019)) { #We need to adjust some other columns
    cColNamesTemp <- cColNames
    cColNamesTemp[c(12, 25,26,27,28)] <- c("TotPotaGPCD", "TotalSup", "WellsSup","SpringsSup","SurfaceSup")
    
    colnames(dfTempFile) <- cColNamesTemp
    
    #Add blank columns for these fields
    cColsAdd <- c("ComSecGPCD", "IndSecGPCD", "InsSecGPCD", "InsSecoUse", "ResSecGPCD", "ResSecoUse", "ComSecoUse", "IndSecoUse", "IndSecoUse")
    for (i in (1:length(cColsAdd))) {
      dfTempFile[[cColsAdd[i]]] = 0
    }
    
  }

  if(i==1){ #first year, creat new data frame
    dfAllYears <- dfTempFile
  } else { #subsequent year, bind to existing records
  
    #bind the latest record to the existing records
    dfAllYears <- rbind(dfAllYears,dfTempFile)
    }

}



#Calculate population from per capita daily use and annual volume
dGalPerAcreFeet <- 325851
dfAllYears$Population <-  dfAllYears$TotalUse * dGalPerAcreFeet / (dfAllYears$TotalGPCD * 365)

#Calculate average percapita use statewide

dfAvgGPCD <- dfAllYears %>% group_by(year) %>% filter(is.na(Population) == FALSE, is.infinite(Population) == FALSE) %>% summarize(TotalUse = sum(TotalUse), 
                                                         Population=sum(Population),
                                                         AvgGPCD = sum(TotalUse)*dGalPerAcreFeet/(365*sum(Population)),
                                                         UsePerPersonPerYear = sum(TotalUse)/sum(Population))


#Grab color palettes
palBlues <- brewer.pal(9, "Blues")
palReds <- brewer.pal(9, "Reds")
palBlueFunc <- colorRampPalette(c(palBlues[3],palBlues[9]))
```

\newpage
## Figure 1. Utah water providers ranked by 2018 total per capita water use.
## Comparison to Salt Lake City use and state conservation targets.
```{r Figure1, echo=FALSE, warning=FALSE, message=FALSE, fig.width=9}

#Bar graph of providers ranked by Total GPCD
#Colors to split into zones. The zones are Salt Lake City's use of 191 gpcd, 2x that, and above
#Highlight listed cities as black bars

#order the data frame by the TotalGPCD field
dfAllYearsSort <- dfAllYears %>% filter(year == 2018) %>% arrange(TotalGPCD)
dfAllYearsSort$Row <- seq(1,nrow(dfAllYearsSort), by=1)

#Cities to highlight
#These cities have a range of Total GPCD, Total GPCD is often similar to Potable GPCD, and people will recognize these providers
#Read in from Excel
dfHighlightCities <- read_excel("HighlightCities.xlsx","HighlightCities","A1:B23")

# The Salt Lake City gallons per person per day reported
nAvgGPCD <- 191 #240 = state average
sAvgText <- "Salt Lake City"

# Define the state conservation goal and calculate the equivalent GPCD that
# will reach the nAvgGPCD vaue
nStateConserveGoal <- 0.2
#Calculate the GPCD that will meet the SLC target
nGPCDToMeetTarget <-  -nAvgGPCD/(nStateConserveGoal - 1)
print(paste("A utility with", sprintf("%0.f",nGPCDToMeetTarget), "gpcd can reduce by state conservation goal of", sprintf("%.0f%%", nStateConserveGoal*100), "and achieve SLC's use"))

# Break at 0, the average, the GPCD to reduce to meet the average, and max
cGPCDBreaks <-  c(0, nAvgGPCD, nGPCDToMeetTarget, max(dfAllYearsSort$TotalGPCD))
#Count water providers within breaks
dfCounts <- hist(dfAllYearsSort$TotalGPCD, breaks = cGPCDBreaks, plot=FALSE) #, breaks = c(0,100,200,300,400,500, 1000, 5000))
dfCounts$CumSum <- nrow(dfAllYearsSort) - cumsum(dfCounts$counts)

#Assign each water provider to a break
dfAllYearsSort$TotGPCDBin <- cut(dfAllYearsSort$TotalGPCD, breaks = cGPCDBreaks, labels = cGPCDBreaks[1:(length(cGPCDBreaks)-1)])
#Sum in each group
dfGPCDByVolume <- dfAllYearsSort %>% group_by(TotGPCDBin, year) %>% summarize(TotalUse = sum(TotalUse))

#Add an additional group that is the Water providers to highlight
dfAllYearsSort$TotGPCDBin <- ifelse(dfAllYearsSort$WRENAME %in% dfHighlightCities$FullName,"Highlight", dfAllYearsSort$TotGPCDBin)

#Joint the Short cities names
dfAllYearsSort <- left_join(dfAllYearsSort, dfHighlightCities, by = c("WRENAME" = "FullName"))
#Convert NAs to ""
dfAllYearsSort$ShortName <- ifelse(is.na(dfAllYearsSort$ShortName),"",dfAllYearsSort$ShortName)

#Add the volume data to the prior counts data
dfCounts$Volume <- dfGPCDByVolume$TotalUse
dfCounts$PercentVolume <- dfCounts$Volume/sum(dfCounts$Volume)
dfCounts$TextFull <- paste("These ",dfCounts$counts, "providers use\n", sprintf("%1.0f%%", dfCounts$PercentVolume*100), "of the statewide total")

dfCounts$Text <- paste(dfCounts$counts, "\nProviders")
dfCounts$TextFull <- paste(dfCounts$counts, "\nProviders\n", "(",sprintf("%1.0f%%", dfCounts$PercentVolume*100), "state volume)")


#Nullify the last text entry
dfCounts$TextFull <- dfCounts$TextFull[1:(length(dfCounts$TextFull)-1)]

dfCounts$Text <- paste(dfCounts$counts, "\nProviders")


#Build a data frame to plot labels for the groups
dfLabels <- data.frame(Rank = dfCounts$CumSum, 
                       Text = dfCounts$Text, 
                       TextFull = dfCounts$TextFull,
                       Counts = dfCounts$counts,
                       FontColor = as.factor(c(1,2,3)))
#Calculate the Mid rank as Cumulative + 1/2 current counts
dfLabels$MidRank <- dfLabels$Rank + dfLabels$Counts/2
dfLabels$YPos <- 1000

dfAllYearsSort$Blank <- ""

#Build a data frame to label the lines
dfLineLabels <- data.frame(Rank = dfLabels$Rank[1:2],
                           GPCD = cGPCDBreaks[2:3],
                           #Label = c("Salt Lake City ~ 191 gpcd", "2x Salt Lake City ~ 382 gpcd"))
                          Label = c(paste(sAvgText, "\n", nAvgGPCD, "gpcd"), paste("120%", sAvgText, "\n", sprintf("%.0f", nGPCDToMeetTarget), "gpcd")))

sBlockLabelColor <- "red"
sBlockSepColor <- "black"
sHighlightColor <- "red"

#Flip the Short Names upside down to plot labels correctly
dfAllYearsSortFlip <- dfAllYearsSort %>% map_df(rev)

ggplot(dfAllYearsSort ) +
  #Bar graph of Total gallons per person per day reversed ranked (largest at left)
  geom_bar(aes(x= reorder(WRENAME, -TotalGPCD), y=TotalGPCD, fill = TotGPCDBin, color = TotGPCDBin), stat="identity") +
  #Plot vertical line to separate blocks
  geom_vline(xintercept = dfCounts$CumSum, color=sBlockSepColor, size = 1, linetype = "dashed") +
  
  #Label the number of providers in each block
  geom_text(data = dfLabels, aes(x = MidRank, y = YPos, label = Text, color = FontColor), size = 4.5) +
  #Label the line breaks
  geom_text(data = dfLineLabels, aes(x = Rank - 1 , y = 2500, label = Label), size = 4, angle = 90, color = sBlockSepColor) +
   
  #Plot some points to see where they are
  #geom_point(data = data.frame(x = c(0,1000), y = c(0,1000)), aes(x=x,y=y), size= 10) +
  
  #geom_text(data=dfAvgGPCD %>% filter(year==2015), aes(x=AvgGPCD+30,y=25), label = "Utah average", size=6, hjust=0) +
  
  #Define scales
  scale_x_discrete(labels = dfAllYearsSortFlip$ShortName) + 
  scale_color_manual(values = c(palBlues[4], palBlues[6], palBlues[8], sHighlightColor)) +
  scale_fill_manual(values = c(palBlues[4], palBlues[6], palBlues[8], sHighlightColor)) +
  scale_y_continuous(breaks = seq(0,4500, by=500)) +
  
  #Turn the Fill guide off
  guides(fill = "none", color = "none") + 
  labs(x="Utah water provider (Ranked)", y="Gallons per person per day") +
  theme(text = element_text(size=14), legend.title=element_blank(), legend.text=element_text(size=18),
        #Remove the minor x grid lines
        panel.grid.minor.x = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        legend.key = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, color = sHighlightColor))


#Prepare a data frame to download to excel that shows:
#   1) Highlight city
#   2) Total GPCD
#   3) Potable GPCD
#   4) Percent reducation to get to Salt Lake City value

#Filter and select to get the data we want
dfReductions <- dfAllYearsSort %>%  filter(ShortName != "") %>% select(ShortName, TotalGPCD, TotPotaGPCD) %>% arrange(-TotalGPCD)
# Calculate % reduction
dfReductions$PercentReduce <- (dfReductions$TotalGPCD - nAvgGPCD)/dfReductions$TotalGPCD

dfReductions$TargetCanMeet <- ifelse(dfReductions$PercentReduce > nStateConserveGoal, "No", "Yes")
# State whether sufficient to meet conservation goal of 20%


#Write the sorted dataframe to csv
write.csv(dfReductions, "dfReductions.csv")