analysis_NORC.Rmd

---
title: "FRIS/IWMS Analysis: No maps or trends plots"
author: "Britta Schumacher"
date: "12/16/2020"
output: 
  html_document:
    theme: flatly
    toc: yes
    toc_float: true
    code_folding: hide
---

# Prepare data for analyses

#### Pull in data from step 11, build-panels-for-analysis.Rmd. First, we'll work with the full-panel-FRIS.RDS to draw summary statistics for each year-by-county and for all years-by-county.

Bring in packages & data built in step 11, build-panels-for-analysis.Rmd
```{r message = F, warning = F}
library(tidyverse)
library(dplyr)
library(ggplot2)
library(corrplot)
library(sp)
library(spdplyr)
library(tmap)
library(rgdal)
library(RColorBrewer)
library(gridExtra)
library(ggplot2)
library(ggpubr)
library(rgeos)
options(scipen=999) # no scientific notation

fris <- readRDS("./data/full-panel-FRIS.RDS")
```

# Build county-level summary statistics

### Find estimated average yield

Code Chunk A):

**DATA IN**: full-panel-FRIS.RDS

Find estimated average yield, water applied, and irrigation productivity across years.

As a reminder, the following data comes from *key variables*: K50, K51, K52, K53 (CORN GRAIN); K60, K61, K62, K63 (CORN SILAGE); K80, K81, K82, K83 (WHEAT); K142, K143 (ALFALFA); K150, K151, K152, K153 (OTHER HAY). These variables represent:
K50: acres harvested (I converted all acres to hectares)
K51: yield (I converted all yield to pounds, and then kilograms and Megagrams in irrigation-productivity.RDS)
K52: water use (I took all water use in acre feet/acre and converted in to inches/acre, cubic cm/hectare, and cubic meters/hectare, again, in irrigation-productivity.RDS).

#### Now, we will take 'statewij' and use it to build all county-level summary statistics as follows (thanks to Matt Fetter at NASS for these exact instructions):

1) Create df with statewij, K51 (acres harvested), K61 (reported yield/acre) [I am adding GEOID, YEAR, etc.]; I converted all reported yields to lbs/acre and kg/acre & lbs/hectare and kg/hectare;
2) create new column: multiply K50*K51 (harvested lbs OR kg);
3) create new column: expanded acres harvested (statewij*K50);
4) create new column: expanded harvested (step 2 *statewij);
5) create summary table: summarise at the county level by sum(expanded harvested) / sum(expanded acres harvested) to get a county-level yield/acre.

```{r A) WEIGHTED YIELD, warning = FALSE, message = FALSE}
### following the instructions given to me on 10/7 by Matt Fetter at NASS
# 1) Create df with statewij, K51 (acres harvested), K61 (reported yield/acre) [I am adding GEOID, YEAR, etc.]; I converted all reported yields to lbs/acre, kg/hectare, and Mg/hectare, and all water use from AF/acre to in/acre, and cm/hectare, m/hectare

weighted <- fris %>% 
  select(GEOID, STPOID, YEAR, STATE_ALPHA, STATE_ABBR, STATEFP, COUNTY_ALPHA, COUNTYFP, statewij, CROP, ACRES_HARV, HECTARES_HARV,  YIELD_LBS, YIELD_KG, YIELD_Mg, WATER_AF, WATER_IN, WATER_CM, WATER_M)  

# 2) create new column(s): multiply K50*K51 (harvested lbs OR kg)
weighted <- weighted %>% 
  mutate(HARVESTED_LBS = ACRES_HARV*YIELD_LBS,
         HARVESTED_KG = HECTARES_HARV*YIELD_KG,
         HARVESTED_Mg = HECTARES_HARV*YIELD_Mg)

# 3) create new column: expanded acres harvested (statewij*K50)
weighted <- weighted %>% 
  mutate(w_ACRES_HARV = ACRES_HARV*statewij,
         w_HECTARES_HARV = (HECTARES_HARV*statewij))

# 4) create new column: expanded harvested (step 2 *statewij)
weighted <- weighted %>% 
  mutate(w_HARVESTED_LBS = HARVESTED_LBS*statewij,
         w_HARVESTED_KG = HARVESTED_KG*statewij,
         w_HARVESTED_Mg = HARVESTED_Mg*statewij)

# 5) create summary table: summarise at the county level by sum(expanded harvested) / sum(expanded acres harvested) to get a county-level yield/acre
yield <- weighted %>% 
    group_by(GEOID, CROP, YEAR) %>% 
    summarise(sum_w_ACRES_HARV = sum(w_ACRES_HARV),
           sum_w_HECTARES_HARV = sum(w_HECTARES_HARV),
           sum_w_HARVESTED_LBS = sum(w_HARVESTED_LBS),
           sum_w_HARVESTED_KG = sum(w_HARVESTED_KG),
           sum_w_HARVESTED_Mg = sum(w_HARVESTED_Mg),
           EST_AVE_YIELD_LBS = sum_w_HARVESTED_LBS/sum_w_ACRES_HARV,
           EST_AVE_YIELD_KG = sum_w_HARVESTED_KG/sum_w_HECTARES_HARV,
           EST_AVE_YIELD_Mg = sum_w_HARVESTED_Mg/sum_w_HECTARES_HARV) 

# Just as a check on things, let's see how different these weighted stats are from the raw
yield_unweight <- weighted %>% 
    group_by(GEOID, CROP, YEAR) %>% 
    summarise(mean_YIELD_LBS = mean(YIELD_LBS),
              mean_YIELD_KG = mean(YIELD_KG),
              mean_YIELD_Mg = mean(YIELD_Mg)) 

# most estimates are very, very different but even the highest estimated alfalfa yield is within reasonable range for the irrigated west (~ 10 tons/acre)

# select just estimated average columns and SAVE
yield_sub <- yield %>% 
  select(GEOID, CROP, YEAR, EST_AVE_YIELD_LBS, EST_AVE_YIELD_KG, EST_AVE_YIELD_Mg)

```

### Find estimated average water applied

Code chunk B):

**DATA IN**: full-panel-FRIS.RDS

We are going to follow thise same instructions as listed above (Steps 2-6) to build weighted - water applied.

Additionally, we include step 7:

7) Merge average estimated yield (yield) and average estimated water use (water) into yield_water  dataframe(). These dataframes consist only of average estimated yields and water uses across different metric (e.g., cubic meters/hectare, kg/hectare) and US (e.g., acre feet/acre, lbs/acre) measurements. They include no raw data.

```{r B) WEIGHTED WATER APPLIED, warning = FALSE, message = FALSE}
### following the instructions given to me on 10/7 by Matt Fetter at NASS
# 1) Create df with statewij, K51 (acres harvested), K61 (reported yield/acre) [I am adding GEOID, YEAR, etc.]; I converted all reported yields to lbs/acre, kg/hectare, and Mg/hectare, and all water use from AF/acre to in/acre, and cm/hectare, m/hectare
## We will use the dataframe created, above in step 2), below:

# 2) create new column(s): multiply K50*K52 (water used AF-IN-CM-M)
weighted <- weighted %>% 
  mutate(WU_AF = ACRES_HARV*WATER_AF,
         WU_IN = ACRES_HARV*WATER_IN,
         WU_CM = HECTARES_HARV*WATER_CM,
         WU_M = HECTARES_HARV*WATER_M)

# 3) create new column: expanded acres harvested (statewij*K50)
## Column already created in 'weighted' df --> w_HECTARES_HARV, w_ACRES_HARV

# 4) create new column: expanded water used (step 2 *statewij)
weighted <- weighted %>% 
  mutate(w_WU_AF = WU_AF*statewij,
         w_WU_IN = WU_IN*statewij,
         w_WU_CM = WU_CM*statewij,
         w_WU_M = WU_M*statewij)


# create summary table: summarise at the county level by sum(expanded harvested) / sum(expanded acres harvested) to get a county-level yield/acre
water <- weighted %>% 
    group_by(GEOID, CROP, YEAR) %>% 
    summarise(sum_w_ACRES_HARV = sum(w_ACRES_HARV),
           sum_w_HECTARES_HARV = sum(w_HECTARES_HARV),
           sum_w_WU_AF = sum(w_WU_AF),
           sum_w_WU_IN = sum(w_WU_IN),
           sum_w_WU_CM = sum(w_WU_CM),
           sum_w_WU_M = sum(w_WU_M),
           EST_AVE_WU_AF = sum_w_WU_AF/sum_w_ACRES_HARV,
           EST_AVE_WU_IN = sum_w_WU_IN/sum_w_ACRES_HARV,
           EST_AVE_WU_CM = sum_w_WU_CM/sum_w_HECTARES_HARV,
           EST_AVE_WU_M = sum_w_WU_M/sum_w_HECTARES_HARV) 

# Just as a check on things, let's see how different these weighted stats are from the raw
water_unweight <- weighted %>% 
    group_by(GEOID, CROP, YEAR) %>% 
    summarise(mean_WU_AF = mean(WATER_AF),
              mean_WU_IN = mean(WATER_IN),
              mean_WU_CM = mean(WATER_CM),
              mean_WU_M = mean(WATER_M)) 

# most estimates are PRETTY close in AF so these estimates feel reasonable

# select just estimated average columns and SAVE
water_sub <- water %>% 
  select(GEOID, CROP, YEAR, EST_AVE_WU_AF, EST_AVE_WU_IN, EST_AVE_WU_CM, EST_AVE_WU_M)

# Merge yield-water use, then save!
yield_water <- merge(yield_sub, water_sub, by = c("GEOID", "CROP", "YEAR"))

```

### Find estimated average irrigation productivity 

Code chunk C):

**DATA IN**: full-panel-FRIS.RDS
**DATA OUT**: sumstat_indyrs.csv & sumstat_allyrs.csv

We will measure irrigation productivity as: weight harvested/unit water applied.

We'll call this step 8 and 9 of the data building process:

8) Build a function to calculate the estimated average yield per unit water applied. Here, we'll take the average estimated yield we built in code chunk A) divided by the average estimated water applied we built in code chunk B) to get an estimated average irrigation productivity. First we will complete this task for individual years (ie. county-year-crop combinations) and save these summary statistics .RDSs out to DEER as 'sumstat_indyrs.RDS'.; then,
9) We will build a function to create the same metrics, but we will summarise across all years. Thus, our function will be modified--we will still mask for all county-crop combinations with fewer than *n* >= 6 observations, but we will be summarising across *ALL* years. We will save these summary statistics .RDSs out to DEER as 'sumstat_allyrs.RDS'.

Again, these dataframes/.RDSs consist only of average estimated yields, water uses, and irrigation productivities across different metric (e.g., cubic meters/hectare, kg/hectare, kg/cubic meter) and US (e.g., acre feet/acre, lbs/acre, lbs/acre foot) measurements. They include no raw data.

```{r C) ESTIMATED AVERAGE IRRIGATION PRODUCTIVITY, warning = FALSE, message = FALSE}
# Use yield and water use at n >= 6/cty-yr, and find estimated average irrigation productivity for each crop-cty-yr combination

ss <- function(df) {
  
  use <- weighted %>% 
  group_by(GEOID, CROP, YEAR) %>% 
    summarise(unique = length(unique((STPOID))),
              max_wateruse_AF = max(WATER_AF),
              max_yield_lbs = max(YIELD_LBS),
              max_acres = max(ACRES_HARV),
              total_weighted_acres = sum(w_ACRES_HARV)) %>% 
    filter(unique >= 6)
    
  use_me <- merge(use, df, by = c("GEOID", "CROP", "YEAR"))
    
  summary_ip <- use_me %>% 
    group_by(GEOID, CROP, YEAR) %>%
    mutate(EST_IP_LBS_AF = EST_AVE_YIELD_LBS/EST_AVE_WU_AF, # estimated average pounds/acre foot
           EST_IP_LBS_IN = EST_AVE_YIELD_LBS/EST_AVE_WU_IN, # estimated average pounds/cubic inch
           EST_IP_KG_CM = EST_AVE_YIELD_KG/EST_AVE_WU_CM, # estimated average kg/cubic cm
           EST_IP_KG_M = EST_AVE_YIELD_KG/EST_AVE_WU_M) %>% # estimated average kg/cubic m
    select(-c(EST_AVE_YIELD_Mg))
}

ip_indyrs <- ss(yield_water)
#saveRDS(ip_indyrs, "./out/DEER/sumstat_indyrs.RDS")
write.csv(ip_indyrs, "./out/DEER/sumstat_indyrs.csv")

# Modify the function above to summarise the 'weighted' dataframe created above across all years. Build out all years withn = 10/cty-yr and n = 5/cty-yr

ss2 <- function(df) {
  
  use <- df %>% 
    group_by(GEOID, CROP) %>% 
    summarise(unique = length(unique((STPOID))),
              max_wateruse_AF = max(WATER_AF),
              max_yield_lbs = max(YIELD_LBS),
              max_acres = max(ACRES_HARV),
              total_weighted_acres = sum(w_ACRES_HARV)) %>% 
    filter(unique >= 6)
  
  use_me <- merge(use, df, by = c("GEOID", "CROP"))
  
  ss_panel <- use_me %>% 
    group_by(GEOID, CROP) %>% 
    summarise(sum_w_ACRES_HARV = sum(w_ACRES_HARV),
           sum_w_HECTARES_HARV = sum(w_HECTARES_HARV),
           sum_w_HARVESTED_LBS = sum(w_HARVESTED_LBS),
           sum_w_HARVESTED_KG = sum(w_HARVESTED_KG),
           sum_w_WU_AF = sum(w_WU_AF),
           sum_w_WU_IN = sum(w_WU_IN),
           sum_w_WU_CM = sum(w_WU_CM),
           sum_w_WU_M = sum(w_WU_M),
           EST_AVE_YIELD_LBS = sum_w_HARVESTED_LBS/sum_w_ACRES_HARV,
           EST_AVE_YIELD_KG = sum_w_HARVESTED_KG/sum_w_HECTARES_HARV,
           EST_AVE_WU_AF = sum_w_WU_AF/sum_w_ACRES_HARV,
           EST_AVE_WU_IN = sum_w_WU_IN/sum_w_ACRES_HARV,
           EST_AVE_WU_CM = sum_w_WU_CM/sum_w_HECTARES_HARV,
           EST_AVE_WU_M = sum_w_WU_M/sum_w_HECTARES_HARV,
           EST_IP_LBS_AF = EST_AVE_YIELD_LBS/EST_AVE_WU_AF, # estimated average pounds/acre foot
           EST_IP_LBS_IN = EST_AVE_YIELD_LBS/EST_AVE_WU_IN, # estimated average pounds/cubic inch
           EST_IP_KG_CM = EST_AVE_YIELD_KG/EST_AVE_WU_CM, # estimated average kg/cubic cm
           EST_IP_KG_M = EST_AVE_YIELD_KG/EST_AVE_WU_M) %>%  # estimated average kg/cubic m
    select(GEOID, CROP, EST_AVE_YIELD_LBS, EST_AVE_YIELD_KG, EST_AVE_WU_AF, EST_AVE_WU_IN, EST_AVE_WU_CM, EST_AVE_WU_M, EST_IP_LBS_AF, EST_IP_LBS_IN, EST_IP_KG_CM, EST_IP_KG_M)
  
  ss_panel <- merge(use, ss_panel, by = c("GEOID", "CROP"))
  
  
}

ip_allyears <- ss2(weighted)
#saveRDS(ip_allyears, "./out/DEER/sumstat_allyrs.RDS")
write.csv(ip_allyears, "./out/DEER/sumstat_allyrs.csv")
```

# Spatial Visualizations: Maps

Code chunk D): *Will rebuild with approved datasets outside the NORC!*

**DATA IN**: sumstat_allyrs.csv and western_ctys.RDS

With the county-level summary data I built in code chunks A), B), and C) we are first going to create county-level maps of estimated average yield, estimated average water use, and estimated irrigation productivity. I follow these steps and have annotated code chunk D) accordingly:

1) Bring in necesary packages, the spatial data for the western US (literally just counties with a "GEOID" column to merge with), built in step 1, and the sumstat_allyrs.csv datasets;
2) Merge the ctys dataframe with fris_sp_allyrs (summary statistics created above with *n* >= 6 unique county-crop observations);
3) Create two functions: one that I used to initially visualize the data with quantile style binning; and a second to create unique bin cut-offs for each  that evenly break (on a zero), where the upper bound of the greatest bin is greater than the most extreme county value, and where the high end of each bin is at least 111% of the bins midpoint value. (It feels inconsistent and untruthful and wrong to leave the greatest bin unbounded, as Matt Fetter requested back in March. Feel free to reach out to discuss this.). These functions take on the following calls:

function_name(dataframe, "CROP", "VARIABLE OF INTEREST", c(BREAKS), "COLOR PALETTE", "TITLE", "./file_path/file_name.png")

All visualizations are directed to a specific folder and have unique filenames. The naming convention is as follows:"

file_name == crop_variable_units_map.png (e.g., alfalfa_yield_lbs_map.png, corngrain_wu_AF_map.png)

4) Build maps for estimated average yield in pounds/acre and kg/hectare;
5) Build maps for estimated average water applied/used in AF/acre, and cubic meters;
6) Build maps for estimated average irrigation productivity in AF/acre and kg/cubic meter.

```{r D) SPATIAL VISUALIZATIONS-IP-YIELD-WU, message = F, warning = F, eval=FALSE}
# Read in packages & data
library(sp)
library(spdplyr)
library(tmap)
library(rgdal)
library(RColorBrewer)
library(gridExtra)
library(ggplot2)
library(ggpubr)
library(rgeos)

# move ss_allyrs.csv and ss_indyrs.csv FROM out/DEER to the results folder

ctys <- readRDS("./data/western_ctys.RDS")
fris_sp_allyrs <- read.csv("./out/DEER/sumstat_allyrs.csv")

# Build spatial across crops using the panel dataset (summarizes across all years)

sp_allyrs_fris <- sp::merge(ctys, fris_sp_allyrs, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)

# use to initially visualize the data -- determine bin cut-offs with a quantile style (each bin will have equal numbers of counties) 
build_sp_allyrs <- function(df, crop, variable, colors, legend, place) {
  
  x <- df %>% 
    filter(CROP == crop)

  crop_yield_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(x) + 
    tm_polygons(col= variable, style = "quantile", n = 6, palette = colors, border.col = "white",
              legend.hist = TRUE, title = legend) +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

  tmap_save(crop_yield_map, place)
  
  return(crop_yield_map)

}

# use to finalize the visualization of the data -- determine bin cut-offs with a quantile-esque style (each bin will have equal-ish numbers of counties). I say quantile-esque because I attempted to satisfy Matt Fetter's request that the top of each bin be 111% of the midpoint. I also pushed the top of the bin beyond the highest value. Matt asked that the top bin be unbounded.... but this feels entirely inconsistent and a bit of a data "lie".
build_sp_allyrs2 <- function(df, crop, variable, breaks, colors, legend, place) {
  
  x <- df %>% 
    filter(CROP == crop)

  crop_yield_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(x) + 
    tm_polygons(col= variable, breaks = breaks, palette = colors, border.col = "white",
              legend.hist = TRUE, title = legend) +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

  tmap_save(crop_yield_map, place)
  
  return(crop_yield_map)

}

###################################################################################################################
### ESTIMATED AVERAGE YIELD
###################################################################################################################
### Estimated average yield across panel - US (lbs)
build_sp_allyrs2(sp_allyrs_fris, "ALFALFA", "EST_AVE_YIELD_LBS", c(4500,6250,7750,9500,11750,14500,26000), "Greens", "Average Estimated\nAlfalfa Yield\n(lbs/acre)", "./viz/mean_yield_panel/alfalfa_yield_lbs_map.png")
build_sp_allyrs2(sp_allyrs_fris, "OTHER_HAY", "EST_AVE_YIELD_LBS", c(1750,4000,5000,6500,8000,10000,20500), "Reds", "Average Estimated\nHay Yield\n(lbs/acre)", "./viz/mean_yield_panel/hay_yield_lbs_map.png")
build_sp_allyrs2(sp_allyrs_fris, "WHEAT", "EST_AVE_YIELD_LBS", c(2500,3250,4250,5250,6500,8000), "Oranges", "Average Estimated\nHay Yield\n(lbs/acre)", "./viz/mean_yield_panel/wheat_yield_lbs_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_GRAIN", "EST_AVE_YIELD_LBS", c(7000,8750,10750,14500), "Blues", "Average Estimated\nCorn Grain Yield\n(lbs/acre)", "./viz/mean_yield_panel/corngrain_yield_lbs_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_SILAGE", "EST_AVE_YIELD_LBS", c(11250,14500,17500,23000), "Purples", "Average Estimated\nCorn Silage Yield\n(lbs/acre)", "./viz/mean_yield_panel/cornsilage_yield_lbs_map.png")

###################################################################################################################
### Estimated average yield across panel - Metric (kg)
build_sp_allyrs2(sp_allyrs_fris, "ALFALFA", "EST_AVE_YIELD_KG", c(5000,7000,8500,10000,12000,15000,29000), "Greens", "Average Estimated\nAlfalfa Yield\n(kg/hectare)", "./viz/mean_yield_panel/alfalfa_yield5_kg_map.png")
build_sp_allyrs2(sp_allyrs_fris, "OTHER_HAY", "EST_AVE_YIELD_KG", c(2000,4000,5500,7000,9000,12000,23000), "Reds", "Average Estimated\nHay Yield\n(kg/hectare)", "./viz/mean_yield_panel/hay_yield5_kg_map.png")
build_sp_allyrs2(sp_allyrs_fris, "WHEAT", "EST_AVE_YIELD_KG", c(3000,4500,5500,6500,7500,8500), "Oranges", "Average Estimated\nWheat Yield\n(kg/hectare)", "./viz/mean_yield_panel/wheat_yield5_kg_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_GRAIN", "EST_AVE_YIELD_KG", c(8250,10000,11000,12000,13000,14000,16250), "Blues", "Average Estimated\nCorn Grain Yield\n(kg/hectare)", "./viz/mean_yield_panel/corngrain_yield5_kg_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_SILAGE", "EST_AVE_YIELD_KG", c(12500,18000,19000,20000,21500,23000, 26000), "Purples", "Average Estimated\nCorn Silage Yield\n(kg/hectare)", "./viz/mean_yield_panel/cornsilage_yield5_kg_map.png")

###################################################################################################################
### ESTIMATED AVERAGE WATER APPLICATION
###################################################################################################################
### Mean Water Application across panel - US
build_sp_allyrs2(sp_allyrs_fris, "ALFALFA", "EST_AVE_WU_AF", c(0.25,1.25,1.75,2.25,2.75,3.5,6.75), "Greens", "Estimated Average Water\nApplication in Alfalfa\n(acre feet/acre)", "./viz/mean_waterapplication_panel/alfalfa_wu_AF_map.png")
build_sp_allyrs2(sp_allyrs_fris, "OTHER_HAY", "EST_AVE_WU_AF", c(0.5,1.,1.5,2,2.5,3,5.25), "Reds", "Estimated Average Water\nApplication in Hay\n(acre feet/acre)", "./viz/mean_waterapplication_panel/hay_wu_AF_map.png")
build_sp_allyrs2(sp_allyrs_fris, "WHEAT", "EST_AVE_WU_AF", c(0.25,1.25,1.5,1.75,2,2.5,4.25), "Oranges", "Estimated Average Water\nApplication in Wheat\n(acre feet/acre)", "./viz/mean_waterapplication_panel/wheat_wu_AF_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_GRAIN", "EST_AVE_WU_AF", c(1,1.25,1.75,2.25,2.75,3.5,5), "Blues", "Estimated Average Water\nApplication in Corn Grain\n(acre feet/acre)", "./viz/mean_waterapplication_panel/corngrain_wu_AF_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_SILAGE", "EST_AVE_WU_AF", c(0.75,1.25,1.75,2.25,2.75,3.5,4.75), "Purples", "Estimated Average Water\nApplication in Corn Silage\n(acre feet/acre)", "./viz/mean_waterapplication_panel/cornsilage_wu_AF_map.png")

###################################################################################################################
### Mean Water Application across panel - Metric (m)
build_sp_allyrs2(sp_allyrs_fris, "ALFALFA", "EST_AVE_WU_M", c(2000,4000,5000,6500,8000,10000,20000), "Greens", "Estimated Average Water\nApplication in Alfalfa\n(cubic meters/hectare)", "./viz/mean_waterapplication_panel/alfalfa_wu5_M_map.png")
build_sp_allyrs2(sp_allyrs_fris, "OTHER_HAY", "EST_AVE_WU_M", c(1500,3500,4500,6000,7500,9000,15750), "Reds", "Estimated Average Water\nApplication in Hay\n(cubic meters/hectare)", "./viz/mean_waterapplication_panel/hay_wu5_M_map.png")
build_sp_allyrs2(sp_allyrs_fris, "WHEAT", "EST_AVE_WU_M", c(1250,3000,4000,5000,6000,8000,12250), "Oranges", "Estimated Average Water\nApplication in Wheat\n(cubic meters/hectare)", "./viz/mean_waterapplication_panel/wheat_wu5_M_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_GRAIN", "EST_AVE_WU_M", c(3500,4500,5500,6500,7500,9000,15250), "Blues", "Estimated Average Water\nApplication in Corn Grain\n(cubic meters/hectare)", "./viz/mean_waterapplication_panel/corngrain_wu5_M_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_SILAGE", "EST_AVE_WU_M", c(2500,4500,5500,6500,7500,9000,15000), "Purples", "Estimated Average Water\nApplication in Corn Silage\n(cubic meters/hectare)", "./viz/mean_waterapplication_panel/cornsilage_wu5_M_map.png")

###################################################################################################################
### IRRIGATION PRODUCTIVITY
###################################################################################################################
### Mean Irrigation Productivity across panel - US (pounds/acre foot)
build_sp_allyrs2(sp_allyrs_fris, "ALFALFA", "EST_IP_LBS_AF", c(2000,3000,4000,5000,6250,8000,23000), "Greens", "Estimated Average Irrigation\nProductivity in Alfalfa\n(pounds/acre foot)", "./viz/mean_ip_panel/alfalfa_ip_LBSperAF_map.png") 
build_sp_allyrs2(sp_allyrs_fris, "OTHER_HAY", "EST_IP_LBS_AF", c(1000,2500,3250,4000,5000,6250,25500), "Reds", "Estimated Average Irrigation\nProductivity in Hay\n(pounds/acre foot)", "./viz/mean_ip_panel/hay_ip_LBSperAF_map.png")
build_sp_allyrs2(sp_allyrs_fris, "WHEAT", "EST_IP_LBS_AF", c(1000,2000,3000,3750,4000,5000,11000), "Oranges", "Estimated Average Irrigation\nProductivity in Wheat\n(pounds/acre foot)", "./viz/mean_ip_panel/wheat_ip_LBSperAF_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_GRAIN", "EST_IP_LBS_AF", c(1750,3000,4000,5000,6250,8000), "Blues", "Estimated Average Irrigation\nProductivity in Corn Grain\n(pounds/acre foot)", "./viz/mean_ip_panel/corngrain_ip_LBSperAF_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_SILAGE", "EST_IP_LBS_AF", c(4000,5000,6250,7500,9250,12000,20000), "Purples", "Estimated Average Irrigation\nProductivity in Corn Silage\n(pounds/acre foot)", "./viz/mean_ip_panel/cornsilage_ip_LBSperAF_map.png")

###################################################################################################################
### Mean Irrigation Productivity across panel - Metric (kg/m)
build_sp_allyrs2(sp_allyrs_fris, "ALFALFA", "EST_IP_KG_M", c(0.75,1,1.25,1.5,1.75,2.25,8.5), "Greens", "Estimated Average Irrigation\nProductivity in Alfalfa\n(kilograms/cubic meter)", "./viz/mean_ip_panel/alfalfa_ip5_KGperM_map.png")
build_sp_allyrs2(sp_allyrs_fris, "OTHER_HAY", "EST_IP_KG_M", c(0.25,0.75,1,1.25,1.5,2.25,9.5), "Reds", "Estimated Average Irrigation\nProductivity in Hay\n(kilograms/cubic meter)", "./viz/mean_ip_panel/hay_ip5_KGperM_map.png")
build_sp_allyrs2(sp_allyrs_fris, "WHEAT", "EST_IP_KG_M", c(0.25,0.75,1,1.25,1.5,1.75,4), "Oranges", "Estimated Average Irrigation\nProductivity in Wheat\n(kilograms/cubic meter)", "./viz/mean_ip_panel/wheat_ip5_KGperM_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_GRAIN", "EST_IP_KG_M", c(0.5,1.25,1.5,1.75,2,2.5,3), "Blues", "Estimated Average Irrigation\nProductivity in Corn Grain\n(kilograms/cubic meter)", "./viz/mean_ip_panel/corngrain_ip5_KGperM_map.png")
build_sp_allyrs2(sp_allyrs_fris, "CORN_SILAGE", "EST_IP_KG_M", c(1.5,2.25,2.75,3.25,3.75,4.5,7.25), "Purples", "Estimated Average Irrigation\nProductivity in Corn Silage\n(kilograms/cubic meter)", "./viz/mean_ip_panel/cornsilage_ip5_KGperM_map.png")
```

Code chunk E): *Will rebuild with approved datasets outside the NORC!*

**DATA**: sumstat_indyrs.csv and western_ctys.RDS

1) I build two functions (nearly identical to the functions built in step 3), code chunk D) for visualizing data specifically in 2018 using the sumstat_indyrs.csv datasets. These functions take on the following calls:

function_name(dataframe, "CROP", YEAR, "VARIABLE OF INTEREST", c(BREAKS), "COLOR PALETTE", "TITLE", "./file_path/file_name.png") and have the same naming convention as above.

2) Build maps for estimated average irrigation productivity in AF/acre and kg/cubic meter.

```{r E) SPATIAL VISUALIZATIONS-IP-2018, message = F, warning = F, eval=FALSE}
####################################################################################################################
### ESTIMATED AVERAGE IRRIGATION PRODUCTIVITY IN 2018
####################################################################################################################

# Build spatial across crops using the _indyrs datasets (summarizes across individual years)
fris_sp_indyrs <- read.csv("./out/DEER/sumstat_indyrs.csv")

sp_indyrs_fris <- sp::merge(ctys, fris_sp_indyrs, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)

build_sp_indyrs <- function(df, crop, year, variable, colors, legend, place) {
  
  x <- df %>% 
    filter(CROP == crop) %>% 
    filter(YEAR == year)

  crop_ip_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(x) + 
    tm_polygons(col=variable, style = "quantile", n = 6, palette = colors, border.col = "white",
              legend.hist = TRUE, title = legend) +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

  tmap_save(crop_ip_map, place)
  
  return(crop_ip_map)

}

build_sp_indyrs2 <- function(df, crop, year, variable, breaks, colors, legend, place) {
  
  x <- df %>% 
    filter(CROP == crop) %>% 
    filter(YEAR == year)

  crop_ip_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(x) + 
    tm_polygons(col=variable, breaks = breaks, palette = colors, border.col = "white",
              legend.hist = TRUE, title = legend) +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

  tmap_save(crop_ip_map, place)
  
  return(crop_ip_map)

}


###################################################################################################################
### Irrigation Productivity in 2018 - US (pounds/AF)
build_sp_indyrs2(sp_indyrs_fris, "ALFALFA", 2018, "EST_IP_LBS_AF", c(2000,3000,3750,4750,5750,7250,16250), "Greens", "Estimated Average Irrigation\nProductivity in Alfalfa\n(pounds/acre foot)", "./viz/mean_ip_2018/alfalfa_ip_LBSperAF_map.png")
build_sp_indyrs2(sp_indyrs_fris, "OTHER_HAY", 2018, "EST_IP_LBS_AF", c(1250,2250,2750,3500,4500,5500,13000), "Reds", "Estimated Average Irrigation\nProductivity in Hay\n(pounds/acre foot)", "./viz/mean_ip_2018/hay_ip_LBSperAF_map.png")
build_sp_indyrs2(sp_indyrs_fris, "WHEAT", 2018, "EST_IP_LBS_AF", c(1500,2000,2750,3500,4500,5500,8500), "Oranges", "Estimated Average Irrigation\nProductivity in Wheat\n(pounds/acre foot)", "./viz/mean_ip_2018/wheat_ip_LBSperAF_map.png")
build_sp_indyrs2(sp_indyrs_fris, "CORN_GRAIN", 2018, "EST_IP_LBS_AF", c(3000,3750,4750,6000,7500,9250), "Blues", "Estimated Average Irrigation\nProductivity in Corn Grain\n(pounds/acre foot)", "./viz/mean_ip_2018/corngrain_ip_LBSperAF_map.png")
build_sp_indyrs2(sp_indyrs_fris, "CORN_SILAGE", 2018, "EST_IP_LBS_AF", c(4000,6000,7500,9250,12000,17500), "Purples", "Estimated Average Irrigation\nProductivity in Corn Silage\n(pounds/acre foot)", "./viz/mean_ip_2018/cornsilage_ip_LBSperAF_map.png")

###################################################################################################################
### Irrigation Productivity in 2018 - Metric (cubic meters/hectare)
build_sp_indyrs2(sp_indyrs_fris, "ALFALFA", 2018, "EST_IP_KG_M", c(0.75,1.25,1.5,1.75,2.25,2.75,6.0), "Greens", "Estimated Average Irrigation\nProductivity in Alfalfa\n(kilograms/cubic meter)", "./viz/mean_ip_2018/alfalfa_ip_KGperM_map.png")
build_sp_indyrs2(sp_indyrs_fris, "OTHER_HAY", 2018, "EST_IP_KG_M", c(0.25,0.75,1,1.25,1.5,2,4.75), "Reds", "Estimated Average Irrigation\nProductivity in Hay\n(kilograms/cubic meter)", "./viz/mean_ip_2018/hay_ip_KGperM_map.png")
build_sp_indyrs2(sp_indyrs_fris, "WHEAT", 2018, "EST_IP_KG_M", c(0.5,0.75,1,1.25,1.5,1.75,3.25), "Oranges", "Estimated Average Irrigation\nProductivity in Wheat\n(kilograms/cubic meter)", "./viz/mean_ip_2018/wheat_ip_KGperM_map.png")
build_sp_indyrs2(sp_indyrs_fris, "CORN_GRAIN", 2018, "EST_IP_KG_M", c(1,1.25,1.5,1.75,2,2.5,3.25), "Blues", "Estimated Average Irrigation\nProductivity in Corn Grain\n(kilograms/cubic meter)", "./viz/mean_ip_2018/corngrain_ip_KGperM_map.png")
build_sp_indyrs2(sp_indyrs_fris, "CORN_SILAGE", 2018, "EST_IP_KG_M", c(1.5,2.25,2.75,3.25,3.75,4.5,6.5), "Purples", "Estimated Average Irrigation\nProductivity in Corn Silage\n(kilograms/cubic meter)", "./viz/mean_ip_2018/cornsilage_ip_KGperM_map.png")

```

Code chunk F): *Will rebuild with approved datasets outside the NORC!*

**DATA**: sumstat_indyrs.csv and western_ctys.RDS

Below, I am interested in finding the differences in estimated average yield, estimated average water application/use, and estimated average irrigation productivity between 2018 (normal year in rainfall) and 2013 (drought year across the west). We want to see how responsive irrigation is to differences in rainfall. So the code works as follows:

1) I build a dataset that merges filtered data from ONLY 2018 and 2013 from the fris_sp_indyrs and fris5_Sp_indyrs dataframes. I then bind/merge these data with our ctys spatial layer so that we can map differences between 2013 and 2018;
2) I build two functions, again, where the first is used to visualize the data in strict quantiles and the second is used to build unique bins using Matt Fetter's criteria. This function mutates new columns into the spatial dataframe that are: 2013_variable - 2018_variable. For example, if the resulting value in the mutated water use column is positive, less water was put on it 2018 than in 2013. I feed the function a "RdBu" palette to visualize differences.;
3) Build maps for estimated average water applied/used in AF/acre and cubic meters;
4) Build maps for estimated average yield in pounds/acre and kg/hectare;
5) Build maps for estimated average irrigation productivity in AF/acre and kg/cubic meter.

```{r F) SPATIAL VISUALIZATIONS-IP-YIELD-WU_DIFFERENCE, message = F, warning = F, eval=FALSE}
######################################################################################################################## DIFFERENCe IN MEAN IP, ESTIMATED AVERAGE YIELD, & ESTIMATED AVERAGE WATER USE BETWEEN 2013 & 2018!
####################################################################################################################

# Build spatial across crops using the _indyrs dataset (summarizes across individual years)
# build-out difference in mean IP/WA/Y between 2013 and 2018

# n >= 10 observations/cty-yr
 # to find differences, use code within function, below

# build spatial
sp_diff <- sp::merge(ctys, y, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)

# build spatial differences maps
build_sp_diff <- function(df, crop, variable, legend, place) {
  
  x <- df %>% 
    filter(CROP == crop) %>% 
    mutate(yield_diff_LBS = (EST_AVE_YIELD_LBS.y - EST_AVE_YIELD_LBS.x), # yield (pounds/acre)
           yield_diff_KG = (EST_AVE_YIELD_KG.y - EST_AVE_YIELD_KG.x), # yield (kg/hectare)
           wa_diff_AF = (EST_AVE_WU_AF.y - EST_AVE_WU_AF.x),    # water applied (acre feet/acre)
           wa_diff_M = (EST_AVE_WU_M.y - EST_AVE_WU_M.x), # water applied (cubic meter/hectare)
           ip_diff_LBSAF = (EST_IP_LBS_AF.y - EST_IP_LBS_AF.x), # irrigation productivity (pounds / AF)
           ip_diff_KGM = (EST_IP_KG_M.y - EST_IP_KG_M.x)) # irrigation productivity (kg / cubic meter)

  crop_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(x) + 
    tm_polygons(col= variable, style = "quantile", n = 6, palette = "RdBu", border.col = "white",
              legend.hist = TRUE, title = legend) +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

  tmap_save(crop_map, place)
  
  return(crop_map)

}

# build spatial differences maps ~ custom breaks
build_sp_diff2 <- function(df, crop, variable, breaks, palette, legend, place) {
  
  x <- df %>% 
    filter(CROP == crop) %>% 
    mutate(yield_diff_LBS = (EST_AVE_YIELD_LBS.y - EST_AVE_YIELD_LBS.x), # yield (pounds/acre)
           yield_diff_KG = (EST_AVE_YIELD_KG.y - EST_AVE_YIELD_KG.x), # yield (kg/hectare)
           wa_diff_AF = (EST_AVE_WU_AF.y - EST_AVE_WU_AF.x),    # water applied (acre feet/acre)
           wa_diff_M = (EST_AVE_WU_M.y - EST_AVE_WU_M.x), # water applied (cubic meter/hectare)
           ip_diff_LBSAF = (EST_IP_LBS_AF.y - EST_IP_LBS_AF.x), # irrigation productivity (pounds / AF)
           ip_diff_KGM = (EST_IP_KG_M.y - EST_IP_KG_M.x)) # irrigation productivity (kg / cubic meter)

  crop_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(x) + 
    tm_polygons(col= variable, breaks = breaks, palette = palette, border.col = "white",
              legend.hist = TRUE, title = legend) +
    tm_legend(outside = TRUE, hist.width = 1.5) +
    tm_layout(frame= FALSE, legend.hist.size = 0.5)

  tmap_save(crop_map, place)
  
  return(crop_map)

}

####################################################################################################################
######  ESTIMATED AVERAGE WATER APPLIED
####################################################################################################################
### Differences in water applied between 2018 (normal year) and 2013 (drought year)

# Average Estimated Water Applied - US AF/acre
build_sp_diff2(sp_diff, "ALFALFA", "wa_diff_AF", c(-2.5,-0.75,-0.5,-0.25,0,0.25,0.5,1,2.5), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC"), "Difference in Water\nApplied in Alfalfa\n(acre feet/acre)", "./viz/mean_diff_waterapplication/alfalfa_diffAF_map.png")
build_sp_diff2(sp_diff, "OTHER_HAY", "wa_diff_AF", c(-1,-0.5,-0.25,0,0.25,0.5,1,2.5), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC"), "Difference in Water\n Applied in Hay\n(acre feet/acre)", "./viz/mean_diff_waterapplication/hay_diffAF_map.png")
build_sp_diff2(sp_diff, "WHEAT", "wa_diff_AF", c(-1.5,-0.5,-0.25,0,0.25,0.5,0.75), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#67A9CF", "#2166AC"), "Difference in Water\nApplied in Wheat\n(acre feet/acre)", "./viz/mean_diff_waterapplication/wheat_diffAF_map.png")
build_sp_diff2(sp_diff, "CORN_GRAIN", "wa_diff_AF", c(-1.5,-0.5,-0.25,0,0.25,0.5,0.75), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#67A9CF", "#2166AC"), "Difference in Water\nApplied in Corn Grain\n(acre feet/acre)", "./viz/mean_diff_waterapplication/corngrain_diffAF_map.png")
build_sp_diff2(sp_diff, "CORN_SILAGE", "wa_diff_AF", c(-2.25,-0.5,-0.25,0,0.25,0.5,0.75,1.25), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC"), "Difference in Water\nApplied in Corn Silage\n(acre feet/acre)", "./viz/mean_diff_waterapplication/cornsilage_diffAF_map.png")

####################################################################################################################
# Average Estimated Water Applied - Metric cubic m/hectare

build_sp_diff2(sp_diff, "ALFALFA", "wa_diff_M", c(-7000,-2500,-1500,-500,0,1000,2000,4000, 7250), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC"), "Difference in Water\nApplied in Alfalfa\n(cubic meters/hectare)", "./viz/mean_diff_waterapplication/alfalfa_diffM_map.png")
build_sp_diff2(sp_diff, "OTHER_HAY", "wa_diff_M", c(-2750,-2000,-1000,-500,0,1000,2000,3000,7000), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC"), "Difference in Water\nApplied in Hay\n(cubic meters/hectare)", "./viz/mean_diff_waterapplication/hay_diffM_map.png")
build_sp_diff2(sp_diff, "WHEAT", "wa_diff_M",  c(-4250,-1500,-500,0,1000,1500, 2250), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#67A9CF", "#2166AC"), "Difference in Water\nApplied in Wheat\n(cubic meters/hectare)", "./viz/mean_diff_waterapplication/wheat_diffM_map.png")
build_sp_diff2(sp_diff, "CORN_GRAIN", "wa_diff_M",  c(-4250,-2500,-1000,0,500,1000,1750), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#67A9CF", "#2166AC"), "Difference in Water\nApplied in Corn Grain\n(cubic meters/hectare)", "./viz/mean_diff_waterapplication/corngrain_diffM_map.png")
build_sp_diff2(sp_diff, "CORN_SILAGE", "wa_diff_M", c(-6500,-2000,-1000,0,1000,2000,3000,3750), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC"), "Difference in Water\nApplied in Corn Silage\n(cubic meters/hectare)", "./viz/mean_diff_waterapplication/cornsilage_diffM_map.png")

####################################################################################################################
######  ESTIMATED AVERAGE YIELD
####################################################################################################################
# Average Estimated Yield - US lbs/acre 

build_sp_diff2(sp_diff, "ALFALFA", "yield_diff_LBS", c(-6000,-2000,-1000,0,1000,2000,4000), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#67A9CF", "#2166AC"), "Difference in Yield \nin Alfalfa (pounds/acre)", "./viz/mean_diff_yield/alfalfa_diffLBS_map.png")
build_sp_diff2(sp_diff, "OTHER_HAY", "yield_diff_LBS", c(-8500,-2000,-1000,-500,0,500,1000,2000,3250), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC"), "Difference in Yield \nin Hay (pounds/acre)", "./viz/mean_diff_yield/hay_diffLBS_map.png")
build_sp_diff2(sp_diff, "WHEAT", "yield_diff_LBS", c(-3750,-2000,-1000, -500,0,500,750,1750), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#2166AC"), "Difference in Yield \nin Wheat (pounds/acre)", "./viz/mean_diff_yield/wheat_diffLBS_map.png")
build_sp_diff2(sp_diff, "CORN_GRAIN", "yield_diff_LBS", c(-4000,-2500,-2000,-1500,-500,0,500,1250), c("#67001F", "#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#2166AC"), "Difference in Yield\nin Corn Grain\n(pounds/acre)", "./viz/mean_diff_yield/corngrain_diffLBS_map.png")
build_sp_diff2(sp_diff, "CORN_SILAGE", "yield_diff_LBS", c(-4500,-2500,-500,0,500,2000,4000), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#67A9CF", "#2166AC"), "Difference in Yield\nin Corn Silage\n(pounds/acre)", "./viz/mean_diff_yield/cornsilage_diffLBS_map.png")

####################################################################################################################
# Average Estimated Yield - Metric kg/hectare

build_sp_diff2(sp_diff, "ALFALFA", "yield_diff_KG", c(-6000,-2000,-1000,-500,0,500,1000,2000,4500), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC"), "Difference in Yield\nin Alfalfa (kilograms/\nhectare)", "./viz/mean_diff_yield/alfalfa_diffKG_map.png")
build_sp_diff2(sp_diff, "OTHER_HAY", "yield_diff_KG", c(-9500,-2000, -1000,-500,0,500,1500,3500), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#2166AC"), "Difference in Yield\nin Hay (kilograms/\nhectare)", "./viz/mean_diff_yield/hay_diffKG_map.png")
build_sp_diff2(sp_diff, "WHEAT", "yield_diff_KG", c(-3000,-1500,-1000,-350,0,650,1000,1750), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#2166AC"), "Difference in Yield\nin Wheat (kilograms/\nhectare)", "./viz/mean_diff_yield/wheat_diffKG_map.png")
build_sp_diff2(sp_diff, "CORN_GRAIN", "yield_diff_KG", c(-4250,-2500,-2250,-1750,-750,0,750,1250), c("#67001F", "#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#2166AC"),"Difference in Yield\nin Corn Grain\n(kilograms/hectare)", "./viz/mean_diff_yield/corngrain_diffKG_map.png")
build_sp_diff2(sp_diff, "CORN_SILAGE", "yield_diff_KG", c(-5000,-3000,-2500,-500,0,1000,3000,4250), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#2166AC"), "Difference in Yield\nin Corn Silage\n(kilograms/hectare)", "./viz/mean_diff_yield/cornsilage_diffKG_map.png")

####################################################################################################################
######  IRRIGATION PRODUCTIVITY
####################################################################################################################
# Average Estimated Irrigation Productivity - US lbs/acre foot 

build_sp_diff2(sp_diff, "ALFALFA", "ip_diff_LBSAF", c(-7750,-3000,-2000,-500,0,500,1000,2000,8500), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC"), "Difference in Irrigation\nProductivity in Alfalfa\n(pounds/acre foot)", "./viz/mean_diff_ip/alfalfa_diffLBSAF_map.png")
build_sp_diff2(sp_diff, "OTHER_HAY", "ip_diff_LBSAF", c(-7500,-3000,-1000,-500,0,500,1000,2750), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#2166AC"), "Difference in Irrigation\nProductivity in Hay\n(pounds/acre foot)", "./viz/mean_diff_ip/hay_diffLBSAF_map.png")
build_sp_diff2(sp_diff, "WHEAT", "ip_diff_LBSAF", c(-3500,-1500,-500,0,250,750,1750), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#67A9CF", "#2166AC"), "Difference in Irrigation\nProductivity in Wheat\n(pounds/acre foot)", "./viz/mean_diff_ip/wheat_diffLBSAF_map.png")
build_sp_diff2(sp_diff, "CORN_GRAIN", "ip_diff_LBSAF", c(-2750,-1250,-750,-250,0,500,2250), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#2166AC"), "Difference in Irrigation\nProductivity in Corn Grain\n(pounds/acre foot)", "./viz/mean_diff_ip/corngrain_diffLBSAF_map.png")
build_sp_diff2(sp_diff, "CORN_SILAGE", "ip_diff_LBSAF", c(-8500,-4000,-2000,-1000,0,1000,2000,2250), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#2166AC"), "Difference in Irrigation\nProductivity in Corn Silage\n(pounds/acre foot)", "./viz/mean_diff_ip/cornsilage_diffLBSAF_map.png")

####################################################################################################################
# Average Estimated Irrigation Productivity - Metric kg/M

build_sp_diff2(sp_diff, "ALFALFA", "ip_diff_KGM", c(-3.0,-1.5,-0.75,-0.25,0,0.25,0.5,1,3.25), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3", "#2166AC"), "Difference in Irrigation\nProductivity in Alfalfa\n(kilograms/cubic meter)", "./viz/mean_diff_ip/alfalfa_diffKGM_map.png")
build_sp_diff2(sp_diff, "OTHER_HAY", "ip_diff_KGM", c(-2.75,-1.5,-0.5,-0.25,0,0.25,0.75,1), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#2166AC"), "Difference in Irrigation\nProductivity in Hay\n(kilograms/cubic meter)", "./viz/mean_diff_ip/hay_diffKGM_map.png")
build_sp_diff2(sp_diff, "WHEAT", "ip_diff_KGM", c(-1.25,-0.5,-0.25,0,0.25,0.5,0.75), c("#B2182B", "#EF8A62", "#FDDBC7", "#D1E5F0", "#67A9CF", "#2166AC"), "Difference in Irrigation\nProductivity in Wheat\n(kilograms/cubic meter)", "./viz/mean_diff_ip/wheat_diffKGM_map.png")
build_sp_diff2(sp_diff, "CORN_GRAIN", "ip_diff_KGM", c(-1,-0.75,-0.5,-0.25,0,0.25,1), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#2166AC"), "Difference in Irrigation\nProductivity in Corn Grain\n(kilograms/cubic meter)", "./viz/mean_diff_ip/corngrain_diffKGM_map.png")
build_sp_diff2(sp_diff, "CORN_SILAGE", "ip_diff_KGM", c(-3.25,-2,-1,-0.5,0,0.5,1), c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#2166AC"), "Difference in Irrigation\nProductivity in Corn Silage\n(kilograms/cubic meter)", "./viz/mean_diff_ip/cornsilage_diffKGM_map.png")
```

# Trends 

Code chunk G): *Will rebuild with approved datasets outside the NORC!*

**DATA**: sumstat_indyrs.csv and western_ctys.RDS

Let's show trends in Estimated Average Yield, Water Use, and Irrigation Productivity over time. These trends are built from the sum_stat_indyrs_X.RDS datasets and are built as line graphs, facetted by state. I am STILL using the same *key variables* (yield, water use).

1) Pull in datasets (again!) selecting only those county-year-crop combinations that show up across all four FRIS years. We will then merge these datasets (have only GEOID as a county identifier) with a distinct() version of the FRIS (has COUNTY_ALPHA, STATE_ABBR). This ensures that we have a "County, State Abbreviation" column to use in our visualizations, below.;
2) I then build two functions for irrigation productivity. I had to do this because my functions were failing to accept new parameters. The first function builds out a line graph  by crop-county-year by state for each metric; the third function builds out by crop-county-year but not by state. I use the third and fourth function to build line graphs for crops that have too few observations to warrant visualization by state.;

These functions take the arguments:

function(data_frame, "CROP", "STATE AABREVIATION", "TITLE--always blank", "Y Axis Label", "./file_path/file_name.png")

Visualizations are sent to crop-specific folders within variable specific trend folders and are built for both metric and US (same graph, different scale) measurements.

Below is a whole lot of code, but its all doing the same exact action on different datasets and different variables and spitting out the results into nested crop-specific and variable-specific folders! Everything below works on estimated average, summarised data--NO raw data is visualized. We are trying to get a sense here of how these different variables change and vary spatially and temporally.

```{r G) TEMPORAL TRENDS VISUALIZATIONS, warning = FALSE, message = FALSE, eval=FALSE}
# Build line graph over time using the _indyrs datasets (to show trends over time)
library(RColorBrewer)
fris_sp_indyrs <- read.csv("./out/DEER/sumstat_indyrs.csv")
fris_sp_indyrs$GEOID <- sprintf("%05d", fris_sp_indyrs$GEOID) 

# find counties in both datasets that have data across all four years of the fris for different crops
trends_fris <- fris_sp_indyrs %>% 
  mutate(count = 1) %>% 
  group_by(GEOID, CROP) %>% 
  summarise(count_yrs = sum(count)) %>% 
  filter(count_yrs == 4) 

# We will use the above as a key to pull data from fris_sp_indyrs and fris5_sp_indyrs, keeping only those county-crop combinations where there are observations in all four years
trends_fris <- merge(trends_fris, fris_sp_indyrs, by = c("GEOID", "CROP"))

# the dataframes above have only county-crop combinations where there is an observation in each year.
# let's merge them with our full-panel-FRIS.RDS to build in a county name identifier
fris <- readRDS("./data/full-panel-FRIS.RDS")
fris_names <- fris %>% 
  select(GEOID, COUNTY_ALPHA, STATE_ABBR) %>% 
  distinct() 

# Convert County name to title case
library(stringr)
fris_names$COUNTY_ALPHA = str_to_title(fris_names$COUNTY_ALPHA)

fris_names <- transform(fris_names, County = paste(COUNTY_ALPHA, STATE_ABBR, sep= ", "))
  
# give these df the pretty names from above!

trends_fris <- merge(trends_fris, fris_names, by = "GEOID")

####################################################################################################################
####### Build functions
####################################################################################################################
# build line graph by crop-county BY STATE
build_trends <- function(df, crop, abbr, var, ylab, place) {
  
  x <- df %>% 
    filter(CROP == crop) %>% 
    filter(STATE_ABBR == abbr)

  var <- enquo(var)
  
  crop_trends <- ggplot(x, aes(x = YEAR, y = !! var)) +
    geom_point(size = 3, aes(color = County)) +
    geom_line(aes(color = County)) +
    scale_color_brewer(palette = "Paired") +
    #guides(col=guide_legend(ncol=2)) +
    xlab("Year") +
    ylab(ylab) +
    scale_x_continuous(breaks = seq(2003, 2018, by = 5)) +
    theme_classic()

  ggsave(place, width=6, height=4,dpi=300, units="in", device='png')
  
  return(crop_trends)

}

# build line graph by crop-county BY STATE (NO PALLETE + 2 column legend)
build_trendsa <- function(df, crop, abbr, var, ylab, place) {
  
  x <- df %>% 
    filter(CROP == crop) %>% 
    filter(STATE_ABBR == abbr)

  var <- enquo(var)
  
  crop_trends <- ggplot(x, aes(x = YEAR, y = !! var)) +
    geom_point(size = 3, aes(color = County)) +
    geom_line(aes(color = County)) +
    #scale_color_brewer(palette = "Paired") +
    guides(col=guide_legend(ncol=2)) +
    xlab("Year") +
    ylab(ylab) +
    scale_x_continuous(breaks = seq(2003, 2018, by = 5)) +
    theme_classic()

  ggsave(place, width=6, height=4,dpi=300, units="in", device='png')
  
  return(crop_trends)

}

####################################################################################################################
####### Irrigation Productivity (Metric)
####################################################################################################################
######## ALFALFA
build_trends(trends_fris, "ALFALFA", "AZ", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/AZ_ip.png")
build_trends(trends_fris, "ALFALFA", "CA", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/CA_ip.png")
build_trends(trends_fris, "ALFALFA", "CO",EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/CO_ip.png")
build_trendsa(trends_fris, "ALFALFA", "ID",EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/ID_ip.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "MT",EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/MT_ip.png") 
build_trends(trends_fris, "ALFALFA", "NV",EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/NV_ip.png")
build_trends(trends_fris, "ALFALFA", "NM",EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/NM_ip.png")
build_trends(trends_fris, "ALFALFA", "OR",EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/OR_ip.png")
build_trendsa(trends_fris, "ALFALFA", "UT",EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/UT_ip.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "WA",EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/WA_ip.png")
build_trends(trends_fris, "ALFALFA", "WY",EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/alfalfa/WY_ip.png")

######## HAY
build_trends(trends_fris, "OTHER_HAY", "AZ", EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/hay/AZ_ip.png")
build_trends(trends_fris, "OTHER_HAY", "CA", EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/hay/CA_ip.png")
build_trends(trends_fris, "OTHER_HAY", "CO", EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/hay/CO_ip.png")
build_trends(trends_fris, "OTHER_HAY", "ID", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/hay/ID_ip.png")
build_trends(trends_fris, "OTHER_HAY", "MT", EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/hay/MT_ip.png")
build_trends(trends_fris, "OTHER_HAY", "NV", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/hay/NV_ip.png")
build_trends(trends_fris, "OTHER_HAY", "OR", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/hay/OR_ip.png")
build_trends(trends_fris, "OTHER_HAY", "UT", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/hay/UT_ip.png")
build_trends(trends_fris, "OTHER_HAY", "WA", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/hay/WA_ip.png")
build_trends(trends_fris, "OTHER_HAY", "WY", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/hay/WY_ip.png")

##### WHEAT
build_trends(trends_fris, "WHEAT", "AZ",EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/wheat/AZ_ip.png")
build_trends(trends_fris, "WHEAT", "CA",EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/wheat/CA_ip.png")
build_trends(trends_fris, "WHEAT", "CO",EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/wheat/CO_ip.png")
build_trends(trends_fris, "WHEAT", "ID",EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/wheat/ID_ip.png") 
build_trends(trends_fris, "WHEAT", "OR",EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/wheat/OR_ip.png")
build_trends(trends_fris, "WHEAT", "UT",EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/wheat/UT_ip.png") 
build_trends(trends_fris, "WHEAT", "WA",EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/wheat/WA_ip.png")


##### Corn Grain
build_trends(trends_fris, "CORN_GRAIN", "CA", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/corngrain/CA_ip.png")
build_trends(trends_fris, "CORN_GRAIN", "CO", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/corngrain/CO_ip.png")
build_trends(trends_fris, "CORN_GRAIN", "ID", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/corngrain/ID_ip.png") 
build_trends(trends_fris, "CORN_GRAIN", "OR", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/corngrain/OR_ip.png")
build_trends(trends_fris, "CORN_GRAIN", "WA", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/corngrain/WA_ip.png")

##### Corn Silage
build_trends(trends_fris, "CORN_SILAGE", "AZ", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/cornsilage/AZ_ip.png")
build_trends(trends_fris, "CORN_SILAGE", "CA", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/cornsilage/CA_ip.png")
build_trends(trends_fris, "CORN_SILAGE", "CO", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/cornsilage/CO_ip.png")
build_trends(trends_fris, "CORN_SILAGE", "ID", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/cornsilage/ID_ip.png") 
build_trends(trends_fris, "CORN_SILAGE", "NM", EST_IP_KG_M,"Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/cornsilage/NM_ip.png")
build_trends(trends_fris, "CORN_SILAGE", "OR", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/cornsilage/OR_ip.png")
build_trends(trends_fris, "CORN_SILAGE", "UT", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/cornsilage/UT_ip.png") 
build_trends(trends_fris, "CORN_SILAGE", "WA", EST_IP_KG_M, "Irrigation Productivity (kilograms/cubic meter)", "./viz/trends_ip/cornsilage/WA_ip.png")

####################################################################################################################
####### Irrigation Productivity (US)
####################################################################################################################
######## ALFALFA
build_trends(trends_fris, "ALFALFA", "AZ", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/AZ_ip_US.png")
build_trends(trends_fris, "ALFALFA", "CA", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/CA_ip_US.png")
build_trends(trends_fris, "ALFALFA", "CO", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/CO_ip_US.png")
build_trendsa(trends_fris, "ALFALFA", "ID", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/ID_ip_US.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "MT", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/MT_ip_US.png") 
build_trends(trends_fris, "ALFALFA", "NV", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/NV_ip_US.png")
build_trends(trends_fris, "ALFALFA", "NM", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/NM_ip_US.png")
build_trends(trends_fris, "ALFALFA", "OR", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/OR_ip_US.png")
build_trendsa(trends_fris, "ALFALFA", "UT", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/UT_ip_US.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "WA", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/WA_ip_US.png")
build_trends(trends_fris, "ALFALFA", "WY", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/alfalfa/WY_ip_US.png")

######## HAY
build_trends(trends_fris, "OTHER_HAY", "AZ", EST_IP_LBS_AF,"Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/hay/AZ_ip_US.png")
build_trends(trends_fris, "OTHER_HAY", "CA", EST_IP_LBS_AF,"Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/hay/CA_ip_US.png")
build_trends(trends_fris, "OTHER_HAY", "CO", EST_IP_LBS_AF,"Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/hay/CO_ip_US.png")
build_trends(trends_fris, "OTHER_HAY", "ID", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/hay/ID_ip_US.png")
build_trends(trends_fris, "OTHER_HAY", "MT", EST_IP_LBS_AF,"Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/hay/MT_ip_US.png")
build_trends(trends_fris, "OTHER_HAY", "NV", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/hay/NV_ip_US.png")
build_trends(trends_fris, "OTHER_HAY", "OR", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/hay/OR_ip_US.png")
build_trends(trends_fris, "OTHER_HAY", "UT", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/hay/UT_ip_US.png")
build_trends(trends_fris, "OTHER_HAY", "WA", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/hay/WA_ip_US.png")
build_trends(trends_fris, "OTHER_HAY", "WY", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/hay/WY_ip_US.png")

##### WHEAT
build_trends(trends_fris, "WHEAT", "AZ",EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/wheat/AZ_ip_US.png")
build_trends(trends_fris, "WHEAT", "CA",EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/wheat/CA_ip_US.png")
build_trends(trends_fris, "WHEAT", "CO",EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/wheat/CO_ip_US.png")
build_trends(trends_fris, "WHEAT", "ID",EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/wheat/ID_ip_US.png") 
build_trends(trends_fris, "WHEAT", "OR",EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/wheat/OR_ip_US.png")
build_trends(trends_fris, "WHEAT", "UT",EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/wheat/UT_ip_US.png") 
build_trends(trends_fris, "WHEAT", "WA",EST_IP_LBS_AF,"Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/wheat/WA_ip_US.png")


##### Corn Grain
build_trends(trends_fris, "CORN_GRAIN", "CA", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/corngrain/CA_ip_US.png")
build_trends(trends_fris, "CORN_GRAIN", "CO", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/corngrain/CO_ip_US.png")
build_trends(trends_fris, "CORN_GRAIN", "ID", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/corngrain/ID_ip_US.png") 
build_trends(trends_fris, "CORN_GRAIN", "OR", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/corngrain/OR_ip_US.png")
build_trends(trends_fris, "CORN_GRAIN", "WA", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/corngrain/WA_ip_US.png")

##### Corn Silage
build_trends(trends_fris, "CORN_SILAGE", "AZ", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/cornsilage/AZ_ip_US.png")
build_trends(trends_fris, "CORN_SILAGE", "CA", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/cornsilage/CA_ip_US.png")
build_trends(trends_fris, "CORN_SILAGE", "CO", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/cornsilage/CO_ip_US.png")
build_trends(trends_fris, "CORN_SILAGE", "ID", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/cornsilage/ID_ip_US.png") 
build_trends(trends_fris, "CORN_SILAGE", "NM", EST_IP_LBS_AF,"Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/cornsilage/NM_ip_US.png")
build_trends(trends_fris, "CORN_SILAGE", "OR", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/cornsilage/OR_ip_US.png")
build_trends(trends_fris, "CORN_SILAGE", "UT", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/cornsilage/UT_ip_US.png") 
build_trends(trends_fris, "CORN_SILAGE", "WA", EST_IP_LBS_AF, "Irrigation Productivity (pounds/acre foot)", "./viz/trends_ip/cornsilage/WA_ip_US.png")

####################################################################################################################
####### ESTIMATED AVERAGE YIELD (Metric)
####################################################################################################################
######## ALFALFA
build_trends(trends_fris, "ALFALFA", "AZ", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/AZ_yield.png")
build_trends(trends_fris, "ALFALFA", "CA", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/CA_yield.png")
build_trends(trends_fris, "ALFALFA", "CO", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/CO_yield.png")
build_trendsa(trends_fris, "ALFALFA", "ID", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/ID_yield.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "MT", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/MT_yield.png") 
build_trends(trends_fris, "ALFALFA", "NV", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/NV_yield.png")
build_trends(trends_fris, "ALFALFA", "NM", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/NM_yield.png")
build_trends(trends_fris, "ALFALFA", "OR", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/OR_yield.png")
build_trendsa(trends_fris, "ALFALFA", "UT", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/UT_yield.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "WA", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/WA_yield.png")
build_trends(trends_fris, "ALFALFA", "WY", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/alfalfa/WY_yield.png")

######## HAY
build_trends(trends_fris, "OTHER_HAY", "AZ", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/hay/AZ_yield.png")
build_trends(trends_fris, "OTHER_HAY", "CA", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/hay/CA_yield.png")
build_trends(trends_fris, "OTHER_HAY", "CO", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/hay/CO_yield.png")
build_trends(trends_fris, "OTHER_HAY", "ID", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/hay/ID_yield.png")
build_trends(trends_fris, "OTHER_HAY", "MT", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/hay/MT_yield.png")
build_trends(trends_fris, "OTHER_HAY", "NV", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/hay/NV_yield.png")
build_trends(trends_fris, "OTHER_HAY", "OR", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/hay/OR_yield.png")
build_trends(trends_fris, "OTHER_HAY", "UT", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/hay/UT_yield.png")
build_trends(trends_fris, "OTHER_HAY", "WA", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/hay/WA_yield.png")
build_trends(trends_fris, "OTHER_HAY", "WY", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/hay/WY_yield.png")

##### WHEAT
build_trends(trends_fris, "WHEAT", "AZ", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/wheat/AZ_yield.png")
build_trends(trends_fris, "WHEAT", "CA", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/wheat/CA_yield.png")
build_trends(trends_fris, "WHEAT", "CO", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/wheat/CO_yield.png")
build_trends(trends_fris, "WHEAT", "ID", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/wheat/ID_yield.png")
build_trends(trends_fris, "WHEAT", "OR", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/wheat/OR_yield.png")
build_trends(trends_fris, "WHEAT", "UT", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/wheat/UT_yield.png")
build_trends(trends_fris, "WHEAT", "WA", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/wheat/WA_yield.png")


##### Corn Grain
build_trends(trends_fris, "CORN_GRAIN", "CA", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/corngrain/CA_yield.png")
build_trends(trends_fris, "CORN_GRAIN", "CO", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/corngrain/CO_yield.png")
build_trends(trends_fris, "CORN_GRAIN", "ID", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/corngrain/ID_yield.png") 
build_trends(trends_fris, "CORN_GRAIN", "OR", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/corngrain/OR_yield.png")
build_trends(trends_fris, "CORN_GRAIN", "WA", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/corngrain/WA_yield.png")

##### Corn Silage
build_trends(trends_fris, "CORN_SILAGE", "AZ", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/cornsilage/AZ_yield.png")
build_trends(trends_fris, "CORN_SILAGE", "CA", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/cornsilage/CA_yield.png")
build_trends(trends_fris, "CORN_SILAGE", "CO", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/cornsilage/CO_yield.png")
build_trends(trends_fris, "CORN_SILAGE", "ID", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/cornsilage/ID_yield.png") 
build_trends(trends_fris, "CORN_SILAGE", "NM", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/cornsilage/NM_yield.png")
build_trends(trends_fris, "CORN_SILAGE", "OR", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/cornsilage/OR_yield.png")
build_trends(trends_fris, "CORN_SILAGE", "UT", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/cornsilage/UT_yield.png") 
build_trends(trends_fris, "CORN_SILAGE", "WA", EST_AVE_YIELD_KG, "Yield (kilograms/hectare)", "./viz/trends_yield/cornsilage/WA_yield.png")


####################################################################################################################
####### Estimated Average Yield (US)
####################################################################################################################
######## ALFALFA
build_trends(trends_fris, "ALFALFA", "AZ", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/AZ_yield_US.png")
build_trends(trends_fris, "ALFALFA", "CA", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/CA_yield_US.png")
build_trends(trends_fris, "ALFALFA", "CO", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/CO_yield_US.png")
build_trendsa(trends_fris, "ALFALFA", "ID", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/ID_yield_US.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "MT", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/MT_yield_US.png") 
build_trends(trends_fris, "ALFALFA", "NV", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/NV_yield_US.png")
build_trends(trends_fris, "ALFALFA", "NM", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/NM_yield_US.png")
build_trends(trends_fris, "ALFALFA", "OR", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/OR_yield_US.png")
build_trendsa(trends_fris, "ALFALFA", "UT", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/UT_yield_US.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "WA", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/WA_yield_US.png")
build_trends(trends_fris, "ALFALFA", "WY", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/alfalfa/WY_yield_US.png")

######## HAY
build_trends(trends_fris, "OTHER_HAY", "AZ", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/hay/AZ_yield_US.png")
build_trends(trends_fris, "OTHER_HAY", "CA", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/hay/CA_yield_US.png")
build_trends(trends_fris, "OTHER_HAY", "CO", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/hay/CO_yield_US.png")
build_trends(trends_fris, "OTHER_HAY", "ID", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/hay/ID_yield_US.png")
build_trends(trends_fris, "OTHER_HAY", "MT", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/hay/MT_yield_US.png")
build_trends(trends_fris, "OTHER_HAY", "NV", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/hay/NV_yield_US.png")
build_trends(trends_fris, "OTHER_HAY", "OR", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/hay/OR_yield_US.png")
build_trends(trends_fris, "OTHER_HAY", "UT", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/hay/UT_yield_US.png")
build_trends(trends_fris, "OTHER_HAY", "WA", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/hay/WA_yield_US.png")
build_trends(trends_fris, "OTHER_HAY", "WY", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/hay/WY_yield_US.png")

##### WHEAT
build_trends(trends_fris, "WHEAT", "AZ", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/wheat/AZ_yield_US.png")
build_trends(trends_fris, "WHEAT", "CA", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/wheat/CA_yield_US.png")
build_trends(trends_fris, "WHEAT", "CO", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/wheat/CO_yield_US.png")
build_trends(trends_fris, "WHEAT", "ID", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/wheat/ID_yield_US.png") 
build_trends(trends_fris, "WHEAT", "OR", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/wheat/OR_yield_US.png")
build_trends(trends_fris, "WHEAT", "UT", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/wheat/UT_yield_US.png") 
build_trends(trends_fris, "WHEAT", "WA", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/wheat/WA_yield_US.png")


##### Corn Grain
build_trends(trends_fris, "CORN_GRAIN", "CA", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/corngrain/CA_yield_US.png")
build_trends(trends_fris, "CORN_GRAIN", "CO", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/corngrain/CO_yield_US.png")
build_trends(trends_fris, "CORN_GRAIN", "ID", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/corngrain/ID_yield_US.png") 
build_trends(trends_fris, "CORN_GRAIN", "OR", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/corngrain/OR_yield_US.png")
build_trends(trends_fris, "CORN_GRAIN", "WA", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/corngrain/WA_yield_US.png")

##### Corn Silage
build_trends(trends_fris, "CORN_SILAGE", "AZ", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/cornsilage/AZ_yield_US.png")
build_trends(trends_fris, "CORN_SILAGE", "CA", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/cornsilage/CA_yield_US.png")
build_trends(trends_fris, "CORN_SILAGE", "CO", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/cornsilage/CO_yield_US.png")
build_trends(trends_fris, "CORN_SILAGE", "ID", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/cornsilage/ID_yield_US.png") 
build_trends(trends_fris, "CORN_SILAGE", "NM", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/cornsilage/NM_yield_US.png")
build_trends(trends_fris, "CORN_SILAGE", "OR", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/cornsilage/OR_yield_US.png")
build_trends(trends_fris, "CORN_SILAGE", "UT", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/cornsilage/UT_yield_US.png") 
build_trends(trends_fris, "CORN_SILAGE", "WA", EST_AVE_YIELD_LBS, "Yield (pounds/acre)", "./viz/trends_yield/cornsilage/WA_yield_US.png")

####################################################################################################################
####### ESTIMATED AVERAGE WATER USE (Metric)
####################################################################################################################
######## ALFALFA
build_trends(trends_fris, "ALFALFA", "AZ", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/AZ_wateruse.png")
build_trends(trends_fris, "ALFALFA", "CA", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/CA_wateruse.png")
build_trends(trends_fris, "ALFALFA", "CO", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/CO_wateruse.png")
build_trendsa(trends_fris, "ALFALFA", "ID", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/ID_wateruse.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "MT", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/MT_wateruse.png") 
build_trends(trends_fris, "ALFALFA", "NV", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/NV_wateruse.png")
build_trends(trends_fris, "ALFALFA", "NM", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/NM_wateruse.png")
build_trends(trends_fris, "ALFALFA", "OR", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/OR_wateruse.png")
build_trendsa(trends_fris, "ALFALFA", "UT", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/UT_wateruse.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "WA", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/WA_wateruse.png")
build_trends(trends_fris, "ALFALFA", "WY", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/alfalfa/WY_wateruse.png")

######## HAY
build_trends(trends_fris, "OTHER_HAY", "AZ", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/AZ_wateruse.png")
build_trends(trends_fris, "OTHER_HAY", "CA", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/CA_wateruse.png")
build_trends(trends_fris, "OTHER_HAY", "CO", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/CO_wateruse.png")
build_trends(trends_fris, "OTHER_HAY", "ID", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/ID_wateruse.png")
build_trends(trends_fris, "OTHER_HAY", "MT", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/MT_wateruse.png")
build_trends(trends_fris, "OTHER_HAY", "NV", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/NV_wateruse.png")
build_trends(trends_fris, "OTHER_HAY", "OR", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/OR_wateruse.png")
build_trends(trends_fris, "OTHER_HAY", "UT", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/UT_wateruse.png")
build_trends(trends_fris, "OTHER_HAY", "WA", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/WA_wateruse.png")
build_trends(trends_fris, "OTHER_HAY", "WY", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/WY_wateruse.png")

##### WHEAT
build_trends(trends_fris, "WHEAT", "AZ", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/wheat/AZ_wateruse.png")
build_trends(trends_fris, "WHEAT", "CA", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/wheat/CA_wateruse.png")
build_trends(trends_fris, "WHEAT", "CO", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/wheat/CO_wateruse.png")
build_trends(trends_fris, "WHEAT", "ID", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/wheat/ID_wateruse.png") 
build_trends(trends_fris, "WHEAT", "OR", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/wheat/OR_wateruse.png")
build_trends(trends_fris, "WHEAT", "UT", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/wheat/UT_wateruse.png") 
build_trends(trends_fris, "WHEAT", "WA", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/wheat/WA_wateruse.png")


##### Corn Grain
build_trends(trends_fris, "CORN_GRAIN", "CA", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/corngrain/CA_wateruse.png")
build_trends(trends_fris, "CORN_GRAIN", "CO", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/corngrain/CO_wateruse.png")
build_trends(trends_fris, "CORN_GRAIN", "ID", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/corngrain/ID_wateruse.png") 
build_trends(trends_fris, "CORN_GRAIN", "OR", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/corngrain/OR_wateruse.png")
build_trends(trends_fris, "CORN_GRAIN", "WA", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/corngrain/WA_wateruse.png")

##### Corn Silage
build_trends(trends_fris, "CORN_SILAGE", "AZ", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/cornsilage/AZ_wateruse.png")
build_trends(trends_fris, "CORN_SILAGE", "CA", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/cornsilage/CA_wateruse.png")
build_trends(trends_fris, "CORN_SILAGE", "CO", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/cornsilage/CO_wateruse.png")
build_trends(trends_fris, "CORN_SILAGE", "ID", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/cornsilage/ID_wateruse.png") 
build_trends(trends_fris, "CORN_SILAGE", "NM", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/cornsilage/NM_wateruse.png")
build_trends(trends_fris, "CORN_SILAGE", "OR", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/cornsilage/OR_wateruse.png")
build_trends(trends_fris, "CORN_SILAGE", "UT", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/cornsilage/UT_wateruse.png") 
build_trends(trends_fris, "CORN_SILAGE", "WA", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/cornsilage/WA_wateruse.png")

####################################################################################################################
####### Estimated Average Water Use (AF)
####################################################################################################################
######## ALFALFA
build_trends(trends_fris, "ALFALFA", "AZ", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/AZ_wateruse_US.png")
build_trends(trends_fris, "ALFALFA", "CA", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/CA_wateruse_US.png")
build_trends(trends_fris, "ALFALFA", "CO", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/CO_wateruse_US.png")
build_trendsa(trends_fris, "ALFALFA", "ID", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/ID_wateruse_US.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "MT", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/MT_wateruse_US.png") 
build_trends(trends_fris, "ALFALFA", "NV", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/NV_wateruse_US.png")
build_trends(trends_fris, "ALFALFA", "NM", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/NM_wateruse_US.png")
build_trends(trends_fris, "ALFALFA", "OR", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/OR_wateruse_US.png")
build_trendsa(trends_fris, "ALFALFA", "UT", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/UT_wateruse_US.png") ## too many counties for palette, use alternate function 
build_trends(trends_fris, "ALFALFA", "WA", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/WA_wateruse_US.png")
build_trends(trends_fris, "ALFALFA", "WY", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/alfalfa/WY_wateruse_US.png")

######## HAY
build_trends(trends_fris, "OTHER_HAY", "AZ", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/hay/AZ_wateruse_US.png")
build_trends(trends_fris, "OTHER_HAY", "CA", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/hay/CA_wateruse_US.png")
build_trends(trends_fris, "OTHER_HAY", "CO", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/hay/CO_wateruse_US.png")
build_trends(trends_fris, "OTHER_HAY", "ID", EST_AVE_WU_M, "Water Application (cubic meters/hectare)", "./viz/trends_wateruse/hay/ID_wateruse_US.png")
build_trends(trends_fris, "OTHER_HAY", "MT", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/hay/MT_wateruse_US.png")
build_trends(trends_fris, "OTHER_HAY", "NV", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/hay/NV_wateruse_US.png")
build_trends(trends_fris, "OTHER_HAY", "NM", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/hay/NM_wateruse_US.png")
build_trends(trends_fris, "OTHER_HAY", "OR", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/hay/OR_wateruse_US.png")
build_trends(trends_fris, "OTHER_HAY", "UT", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/hay/UT_wateruse_US.png")
build_trends(trends_fris, "OTHER_HAY", "WA", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/hay/WA_wateruse_US.png")
build_trends(trends_fris, "OTHER_HAY", "WY", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/hay/WY_wateruse_US.png")

##### WHEAT
build_trends(trends_fris, "WHEAT", "AZ", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/wheat/AZ_wateruse_US.png")
build_trends(trends_fris, "WHEAT", "CA", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/wheat/CA_wateruse_US.png")
build_trends(trends_fris, "WHEAT", "CO", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/wheat/CO_wateruse_US.png")
build_trends(trends_fris, "WHEAT", "ID", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/wheat/ID_wateruse_US.png") 
build_trends(trends_fris, "WHEAT", "OR", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/wheat/OR_wateruse_US.png")
build_trends(trends_fris, "WHEAT", "UT", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/wheat/UT_wateruse_US.png") 
build_trends(trends_fris, "WHEAT", "WA", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/wheat/WA_wateruse_US.png")


##### Corn Grain
build_trends(trends_fris, "CORN_GRAIN", "CA", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/corngrain/CA_wateruse_US.png")
build_trends(trends_fris, "CORN_GRAIN", "CO", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/corngrain/CO_wateruse_US.png")
build_trends(trends_fris, "CORN_GRAIN", "ID", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/corngrain/ID_wateruse_US.png") 
build_trends(trends_fris, "CORN_GRAIN", "OR", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/corngrain/OR_wateruse_US.png")
build_trends(trends_fris, "CORN_GRAIN", "WA", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/corngrain/WA_wateruse_US.png")

##### Corn Silage
build_trends(trends_fris, "CORN_SILAGE", "AZ", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/cornsilage/AZ_wateruse_US.png")
build_trends(trends_fris, "CORN_SILAGE", "CA", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/cornsilage/CA_wateruse_US.png")
build_trends(trends_fris, "CORN_SILAGE", "CO", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/cornsilage/CO_wateruse_US.png")
build_trends(trends_fris, "CORN_SILAGE", "ID", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/cornsilage/ID_wateruse_US.png") 
build_trends(trends_fris, "CORN_SILAGE", "NM", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/cornsilage/NM_wateruse_US.png")
build_trends(trends_fris, "CORN_SILAGE", "OR", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/cornsilage/OR_wateruse_US.png")
build_trends(trends_fris, "CORN_SILAGE", "UT", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/cornsilage/UT_wateruse_US.png") 
build_trends(trends_fris, "CORN_SILAGE", "WA", EST_AVE_WU_AF, "Water Application (acre feet/acre)", "./viz/trends_wateruse/cornsilage/WA_wateruse_US.png")
```

# Surpising counties

**DATA IN**: sumstat_indyrs.csv
**DATA OUT**: alfalfa_SDcat.csv

Let's find counties that lie +1/-1 SD from the mean of Estimated Average Irrigation Productivity in Alfalfa (we'll use these county distinctions for cty-level attribute analysis happening outside of the NORC).

Code Chunk H):

1) Pull in sumstat_indyrs_X.RDS (again!). Filter for only alfalfa data in 2018 and reduce the dataframes to contain ONLY variables of interest!;
2) Find the mean and standard deviations across all observations.;
3) Use these calculated means and SDs to build new columns "above" = +1SD, "below" = -1SD and "within" (well, everything else) from the mean for water use and irrigation productivity;
4) Save these datasets to request out of the NORC.

Here, we are interested in exploring what demographics characterize these counties--we'll do this analysis later, outside of the NORC using the counties identified in the .RDS files we requested!

```{r H) SURPRISING CTYS, warning = FALSE, message = FALSE, eval=FALSE}
fris_sp_indyrs <- read.csv("./out/DEER/sumstat_indyrs.csv")

# select only for alfalfa in 2018
alfalfa_ss <- fris_sp_indyrs %>% filter(CROP == "ALFALFA") %>% filter(YEAR == 2018) 
  
# now let's reduce the dataframe (only interested in irrigation productivity AND water use)
alfalfa_ss <- alfalfa_ss %>% select(-c("EST_AVE_WU_IN", "EST_AVE_WU_CM", "EST_IP_LBS_IN", "EST_IP_KG_CM"))

# now let's find the mean and sd of water use and ip across counties in 2018
# water use (US)
mean_WU_AF <- mean(alfalfa_ss$EST_AVE_WU_AF)
mean_WU_AF
sd_WU_AF <- sd(alfalfa_ss$EST_AVE_WU_AF)
sd_WU_AF

# yield (US)
mean_Y_LBS <- mean(alfalfa_ss$EST_AVE_YIELD_LBS)
mean_Y_LBS
sd_Y_LBS <- sd(alfalfa_ss$EST_AVE_YIELD_LBS)
sd_Y_LBS

# irrigation productivity (US)
mean_ip_US <- mean(alfalfa_ss$EST_IP_LBS_AF)
mean_ip_US
sd_ip_US <- sd(alfalfa_ss$EST_IP_LBS_AF)
sd_ip_US

# Now let's find counties above/below 1SD from the mean in each category
alfalfa_ss$SD_cat_WU <- ifelse(alfalfa_ss$EST_AVE_WU_AF > (mean_WU_AF + sd_WU_AF), "above",
                            ifelse(alfalfa_ss$EST_AVE_WU_AF < (mean_WU_AF - sd_WU_AF), "below",
                            "within"))

alfalfa_ss$SD_cat_Y <- ifelse(alfalfa_ss$EST_AVE_YIELD_LBS > (mean_Y_LBS + sd_Y_LBS), "above",
                            ifelse(alfalfa_ss$EST_AVE_YIELD_LBS < (mean_Y_LBS - sd_Y_LBS), "below",
                            "within"))

alfalfa_ss$SD_cat_IP <- ifelse(alfalfa_ss$EST_IP_LBS_AF > (mean_ip_US + sd_ip_US), "above",
                            ifelse(alfalfa_ss$EST_IP_LBS_AF < (mean_ip_US - sd_ip_US), "below",
                            "within"))

# save full .RDSs and reduced .RDSs with only county identifier and +/-SD category
write.csv(alfalfa_ss,"./out/DEER/alfalfa_SDcat.csv")
```

# Barriers to Implementing Improvements

**DATA IN**: barriers-to-improvements.RDS & western_ctys.RDS
**DATA OUT**: barrier-prop-allyrs.csv & barrier-prop-indyrs.csv

Let's determine which Barriers are greatest by county, facetted by State!

Code chunk I):

1) Ensure that all *key variables*: K1070-K1078 (2013 FRIS and 2018 IWMS) and K695 (2003 & 2008 FRIS) are weighted properly by their farm-level weight, statewij (responses from K695 were separated into K1070-K1078 to merge the four survey years in a different .Rmd).;
2) I summarise across all years and individual years, and all data is masked if *n* >= 6 unique STPOIDs was not satisfied. 
3) I build a proportional bar chart facetted by state where the proportional bar chart shows *only* the proportion of responses in each category across the four census years to give stakeholders a sense of how these vary across western counties. Variables are both categorical (barriers cited) and continuous (proportions; raw numbers are NOT reported).;
4) I take the barrier with the highest proportion or responses and map those responses spatially across all years and in 2018.

We will use the tables of proportions requested to DEER to discuss differences visualized in the proportional bar chart. 
```{r I) BARRIERS-PROP, warning = FALSE, message = FALSE, eval=FALSE}
barriers_key <- readRDS("./data/barriers-to-improvements.RDS")

# create key column -- if 0, respondent did not answer to barriers questions
barriers_key <- barriers_key %>% 
  mutate(empty = K1070+K1071+K1072+K1073+K1074+K1075+K1076+K1077+K1078) %>% 
  filter(!COUNTY_ALPHA == 0) #dunno why these are here, but they are not real GEOIDs

# create new df where if key above is 0, the row is deleted
barriers <- barriers_key %>% filter(empty != 0)

# LEt's weight everything by statewij
barriers <- barriers %>% 
  mutate(K1070_w = K1070*statewij, 
            K1071_w = K1071*statewij, 
            K1072_w = K1072*statewij, 
            K1073_w = K1073*statewij, 
            K1074_w = K1074*statewij, 
            K1075_w = K1075*statewij, 
            K1076_w = K1076*statewij, 
            K1077_w = K1077*statewij,
            K1078_w = K1078*statewij) 

# summarise 

sum_K1070 <- barriers %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1070_w)) %>% 
  summarise(K1070 = sum(K1070_w))
sum_K1071 <- barriers %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1071_w)) %>% 
  summarise(K1071 = sum(K1071_w))
sum <- merge(sum_K1070, sum_K1071, by = c("GEOID", "YEAR"), all=T)
sum_K1072 <- barriers %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1072_w)) %>% 
  summarise(K1072 = sum(K1072_w))
sum <- merge(sum, sum_K1072, by = c("GEOID", "YEAR"), all=T)
sum_K1073 <- barriers %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1073_w)) %>% 
  summarise(K1073 = sum(K1073_w))
sum <- merge(sum, sum_K1073, by = c("GEOID", "YEAR"), all=T)
sum_K1074 <- barriers %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1074_w)) %>% 
  summarise(K1074 = sum(K1074_w))
sum <- merge(sum, sum_K1074, by = c("GEOID", "YEAR"), all=T)
sum_K1075 <- barriers %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1075_w)) %>% 
  summarise(K1075 = sum(K1075_w))
sum <- merge(sum, sum_K1075, by = c("GEOID", "YEAR"), all=T)
sum_K1076 <- barriers %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1076_w)) %>% 
  summarise(K1076 = sum(K1076_w))
sum <- merge(sum, sum_K1076, by = c("GEOID", "YEAR"), all=T)
sum_K1077 <- barriers %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1077_w)) %>% 
  summarise(K1077 = sum(K1077_w))
sum <- merge(sum, sum_K1077, by = c("GEOID", "YEAR"), all=T)
sum_K1078 <- barriers %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1078_w)) %>% 
  summarise(K1078 = sum(K1078_w))
sum <- merge(sum, sum_K1078, by = c("GEOID", "YEAR"), all=T)

# Number of contributing records for each county
records <- barriers %>% 
  group_by(GEOID, YEAR) %>%
  summarise(records = length(STPOID))
sum <- merge(sum, records, by = c("GEOID", "YEAR"), all=T)

sum$n <- rowSums(sum[, 3:11], na.rm=TRUE)

# aggregate to the county-level (across FRIS years)
l <- sum
l[is.na(l)] <- 0
length <- l %>% 
  group_by(GEOID) %>% 
  summarise(K1070 = sum(K1070), 
            K1071 = sum(K1071), 
            K1072 = sum(K1072), 
            K1073 = sum(K1073), 
            K1074 = sum(K1074), 
            K1075 = sum(K1075), 
            K1076 = sum(K1076), 
            K1077 = sum(K1077),
            K1078 = sum(K1078),
            records = sum(records)) 

length$n <- rowSums(length[, 2:10], na.rm=TRUE)

# aggregate to the county-level (across individual FRIS years)
length_i <- sum

# mask all columns calculated using <6 observations across all years
length_sub <- length %>% 
  mutate(K1070 = ifelse(records < 6, NA, K1070),
         K1071 = ifelse(records < 6, NA, K1071),
         K1072 = ifelse(records < 6, NA, K1072),
         K1073 = ifelse(records < 6, NA, K1073),
         K1074 = ifelse(records < 6, NA, K1074),
         K1075 = ifelse(records < 6, NA, K1075),
         K1076 = ifelse(records < 6, NA, K1076),
         K1077 = ifelse(records < 6, NA, K1077),
         K1078 = ifelse(records < 6, NA, K1078))

# mask all columns calculated using <6 observations across individual years
length_i_sub <- length_i %>% 
  mutate(K1070 = ifelse(records < 6, NA, K1070),
         K1071 = ifelse(records < 6, NA, K1071),
         K1072 = ifelse(records < 6, NA, K1072),
         K1073 = ifelse(records < 6, NA, K1073),
         K1074 = ifelse(records < 6, NA, K1074),
         K1075 = ifelse(records < 6, NA, K1075),
         K1076 = ifelse(records < 6, NA, K1076),
         K1077 = ifelse(records < 6, NA, K1077),
         K1078 = ifelse(records < 6, NA, K1078))

# find the proportion distributed across FRIS years
proportion <- length_sub %>% 
  group_by(GEOID) %>% 
  summarise(K1070 = K1070/n*100, 
            K1071 = K1071/n*100, 
            K1072 = K1072/n*100, 
            K1073 = K1073/n*100, 
            K1074 = K1074/n*100, 
            K1075 = K1075/n*100, 
            K1076 = K1076/n*100, 
            K1077 = K1077/n*100, 
            K1078 = K1078/n*100,
            records = records)

# find the proportion distributed across individual FRIS years
proportion_i <- length_i_sub %>% 
  group_by(GEOID, YEAR) %>% 
  summarise(K1070 = K1070/n*100, 
            K1071 = K1071/n*100, 
            K1072 = K1072/n*100, 
            K1073 = K1073/n*100, 
            K1074 = K1074/n*100, 
            K1075 = K1075/n*100, 
            K1076 = K1076/n*100, 
            K1077 = K1077/n*100, 
            K1078 = K1078/n*100,
            records = records)

# All NAs in 'proportion' result from having fewer than 6 unique POIDs in a county over the 4 FRIS years
# All NAs in 'proportion_i' results from having fewer than 6 unique POIDs in a county in each FRIS year

# Do all proportions add to 100%? YES!
proportion$total <- rowSums(proportion[, 2:10], na.rm=TRUE)
proportion_i$total <- rowSums(proportion_i[, 3:11], na.rm=TRUE)

# found proportion to two significant digits
proportion[,-1] <- round(proportion[,-1],2)
proportion_i[,-1] <- round(proportion_i[,-1],2)

# How many counties are represented each year if we do this thing yearly? ie. why are we aggregating based on the four years of data we have for visualizations?
barr_2003 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2003) # 72/408
barr_2008 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2008) # 73/408
barr_2013 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2013) # 123/408
barr_2018 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2018) # 160/408

# Let's build in the State so we can facet by State for the proportional bar chart!
state <- barriers_key %>% select(GEOID, COUNTY_ALPHA, STATE_ALPHA) %>% distinct()

# let's merge so we have all county-states represented
# we'll build all years AND 2018!
prop_merge <- merge(state, proportion, by = "GEOID", all =T)
prop_merge_save <- prop_merge %>% select(!records)

### save this proportion table for export!
#saveRDS(prop_merge_save, "./out/DEER/barrier-prop-allyrs.RDS")
write.csv(prop_merge_save, "./out/DEER/barrier-prop-allyrs.csv")
# build individual years and row bind!
prop_merge_18 <- merge(state, barr_2018, by="GEOID", all=T)
prop_merge_13 <- merge(state, barr_2013, by="GEOID", all=T)
prop_merge_08 <- merge(state, barr_2008, by="GEOID", all=T)
prop_merge_03 <- merge(state, barr_2003, by="GEOID", all=T)

prop_merge_i <- rbind(prop_merge_03, prop_merge_08, prop_merge_13, prop_merge_18)
prop_merge_i_save <- prop_merge_i %>% select(!records)

### save this proportion table for export!
#saveRDS(prop_merge_i_save, "./out/DEER/barrier-prop-indyrs.RDS")
write.csv(prop_merge_i_save, "./out/DEER/barrier-prop-indyrs.csv")
################################################################################################################################################# BUILD FACETTED PROPORTIONAL BAR CHARTS
####################################################################################################################
# build proportion table above into suitable data format for proportion barchart (ie. long format)
v <- prop_merge_save %>% select(!c("total")) %>% gather(category, proportion, 4:12)
w <- prop_merge_18 %>% select(!YEAR) %>% select(!c("records", "total")) %>% gather(category, proportion, 4:12)


barriers_plots <- function(df, state, title, place) {
  
  x <- df %>% 
    filter(STATE_ALPHA == state) 

  barriers_plot <- ggplot(x, aes(fill=category, y=proportion,x=COUNTY_ALPHA)) +
    geom_bar(position = "fill", stat = "identity") +
    scale_fill_manual(values = c("#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3", "#A6D854", "#FFD92F", "#E5C494", "#696969", "#2F2F2F"), name = "Category") +
    #scale_fill_brewer(type = "qual", palette = 3, name = "Category") +
    ggtitle(title) +
    theme_classic() +
    xlab("County") +
    ylab("Proportion") +
    coord_flip() 


  ggsave(place, height = 7, width = 7, dpi=300, units="in", device='png')
  
  return(barriers_plot)

}

# All years
barriers_plots(v, "ARIZONA", "Arizona", "./viz/barriers/barriers_AZ.png") 
barriers_plots(v, "CALIFORNIA", "California", "./viz/barriers/barriers_CA.png") 
barriers_plots(v, "COLORADO", "Colorado", "./viz/barriers/barriers_CO.png")
barriers_plots(v, "IDAHO", "Idaho", "./viz/barriers/barriers_ID.png")
barriers_plots(v, "MONTANA", "Montana", "./viz/barriers/barriers_MT.png")
barriers_plots(v, "NEW MEXICO", "New Mexico", "./viz/barriers/barriers_NM.png")
barriers_plots(v, "NEVADA", "Nevada", "./viz/barriers/barriers_NV.png") 
barriers_plots(v, "OREGON", "Oregon", "./viz/barriers/barriers_OR.png")
barriers_plots(v, "UTAH", "Utah", "./viz/barriers/barriers_UT.png") 
barriers_plots(v, "WASHINGTON", "Washington", "./viz/barriers/barriers_WA.png")
barriers_plots(v, "WYOMING", "Wyoming", "./viz/barriers/barriers_WY.png")

# 2018
barriers_plots(w, "ARIZONA", "Arizona", "./viz/barriers/barriers_AZ18.png") 
barriers_plots(w, "CALIFORNIA", "California", "./viz/barriers/barriers_CA18.png") 
barriers_plots(w, "COLORADO", "Colorado", "./viz/barriers/barriers_CO18.png")
barriers_plots(w, "IDAHO", "Idaho", "./viz/barriers/barriers_ID18.png")
barriers_plots(w, "MONTANA", "Montana", "./viz/barriers/barriers_MT18.png")
barriers_plots(w, "NEW MEXICO", "New Mexico", "./viz/barriers/barriers_NM18.png")
barriers_plots(w, "NEVADA", "Nevada", "./viz/barriers/barriers_NV18.png") 
barriers_plots(w, "OREGON", "Oregon", "./viz/barriers/barriers_OR18.png")
barriers_plots(w, "UTAH", "Utah", "./viz/barriers/barriers_UT18.png") 
barriers_plots(w, "WASHINGTON", "Washington", "./viz/barriers/barriers_WA18.png")
barriers_plots(w, "WYOMING", "Wyoming", "./viz/barriers/barriers_WY18.png")
################################################################################################################################################# BUILD MAP OF MOST IMPORTANT BARRIER
####################################################################################################################
# Select most important source by proportion for each county
find <- v %>% group_by(GEOID) %>% filter(proportion == max(proportion))
find_18 <- w %>% group_by(GEOID) %>% filter(proportion == max(proportion))

# Number of counties in each category
find_summary <- find %>% 
  group_by(category) %>% 
  summarise(n = n())
find_summary

find18_summary <- find_18 %>% 
  group_by(category) %>% 
  summarise(n = n())
find18_summary

# Change colnames for visuzaliation
colnames(find)[4] <- "Barrier to Implementation"
colnames(find_18)[4] <- "Barrier to Implementation"
ctys <- readRDS("./data/western_ctys.RDS")

barrier_sp <- sp::merge(ctys, find, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)
barrier18_sp <- sp::merge(ctys, find_18, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)

# Make a map for summarised data across all years
barrier_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(barrier_sp) + 
    tm_polygons("Barrier to Implementation", palette = c("#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3", "#A6D854", "#FFD92F", "#E5C494", "#696969", "#2F2F2F"), border.col = "white")  +
    tm_legend(outside = TRUE) +
    tm_layout(frame= FALSE)

tmap_save(barrier_map, "./viz/barriers/barrier_map.png")

# Make a map for just 2018
barrier18_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(barrier18_sp) + 
    tm_polygons("Barrier to Implementation", palette = c("#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3", "#A6D854", "#FFD92F", "#E5C494", "#696969", "#2F2F2F"), border.col = "white")  +
    tm_legend(outside = TRUE) +
    tm_layout(frame= FALSE)

tmap_save(barrier18_map, "./viz/barriers/barrier18_map.png")

#display.brewer.all(n=NULL, type="all",select=NULL, colorblindFriendly=T)
```

# Scheduling Methods

**DATA IN**: scheduling-practices.RDS & western_ctys.RDS
**DATA OUT**: sched-prop-allyrs.csv & sched-prop-indyrs.csv

Let's determine which Scheduling Methods are most important by county, facetted by State!

Code chunk J):

1) Ensure that all *key variables*: K1020-K1029 (2013 FRIS & 2018 IWMS) and K527 (2003 & 2008 FRIS) are weighted properly by their farm-level weight, statewij (responses from K527 were separated into K1020-K1029 to merge the four survey years in a different .Rmd).;
2) I summarise across all years and individual years, and all data is masked if *n* >= 6 was not satisfied.;
3) I build a proportional bar chart facetted by state where the proportional bar chart shows *only* the proportion of responses in each category across the four census years to give stakeholders a sense of how these vary across western counties. Variables are both categorical (scheduling method cited) and continuous (proportions; raw numbers are NOT reported).;
4) I take the scheduling method with the highest proportion of responses and map those responses spatially across all years and in only 2018.

We will use the tables of proportions (again, NO RAW DATA) requested to DEER to discuss differences visualized in the proportional bar chart. 
```{r J) SCHEDULING METHODS-PROP, warning = FALSE, message = FALSE, eval=FALSE}

sched_key <- readRDS("./data/scheduling-practices.RDS")

sched_key[is.na(sched_key)] <- 0
# create key column -- if 0, respondent did not answer to scheduling methods questions
sched_key <- sched_key %>% 
  mutate(empty = K1020+K1021+K1022+K1023+K1024+K1025+K1026+K1027+K1028+K1029) %>% 
  filter(!COUNTY_ALPHA == 0) #dunno why these are here, but they are not real GEOIDs

# create new df where if key above is 0, the row is deleted
sched <- sched_key %>% filter(empty != 0)

sched <- sched %>% 
  mutate(K1020_w = K1020*statewij, 
            K1021_w = K1021*statewij, 
            K1022_w = K1022*statewij, 
            K1023_w = K1023*statewij, 
            K1024_w = K1024*statewij, 
            K1025_w = K1025*statewij, 
            K1026_w = K1026*statewij, 
            K1027_w = K1027*statewij,
            K1028_w = K1028*statewij,
            K1029_w = K1029*statewij) 

# summarise (ensure that plyr is detached to summarise across groups)
sum_K1020 <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1020_w)) %>% 
  summarise(K1020 = sum(K1020_w))
sum_K1021 <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1021_w)) %>% 
  summarise(K1021 = sum(K1021_w))
sum <- merge(sum_K1020, sum_K1021, by = c("GEOID", "YEAR"), all=T)
sum_K1022 <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1022_w)) %>% 
  summarise(K1022 = sum(K1022_w))
sum <- merge(sum, sum_K1022, by = c("GEOID", "YEAR"), all=T)
sum_K1023 <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1023_w)) %>% 
  summarise(K1023 = sum(K1023_w))
sum <- merge(sum, sum_K1023, by = c("GEOID", "YEAR"), all=T)
sum_K1024 <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1024_w)) %>% 
  summarise(K1024 = sum(K1024_w))
sum <- merge(sum, sum_K1024, by = c("GEOID", "YEAR"), all=T)
sum_K1025 <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1025_w)) %>% 
  summarise(K1025 = sum(K1025_w))
sum <- merge(sum, sum_K1025, by = c("GEOID", "YEAR"), all=T)
sum_K1026 <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1026_w)) %>% 
  summarise(K1026 = sum(K1026_w))
sum <- merge(sum, sum_K1026, by = c("GEOID", "YEAR"), all=T)
sum_K1027 <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1027_w)) %>% 
  summarise(K1027 = sum(K1027_w))
sum <- merge(sum, sum_K1027, by = c("GEOID", "YEAR"), all=T)
sum_K1028 <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1028_w)) %>% 
  summarise(K1028 = sum(K1028_w))
sum <- merge(sum, sum_K1028, by = c("GEOID", "YEAR"), all=T)
sum_K1029 <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1029_w)) %>% 
  summarise(K1029 = sum(K1029_w))
sum <- merge(sum, sum_K1029, by = c("GEOID", "YEAR"), all=T)

# Number of contributing records for each county
records <- sched %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(GEOID)) %>% 
  summarise(records = length(STPOID))
sum <- merge(sum, records, by = c("GEOID", "YEAR"), all=T)

sum$n <- rowSums(sum[, 3:12], na.rm=TRUE)

# aggregate to the county-level (across FRIS years)
l <- sum
l[is.na(l)] <- 0
length <- l %>% 
  dplyr::group_by(GEOID) %>% 
  summarise(K1020 = sum(K1020), 
            K1021 = sum(K1021), 
            K1022 = sum(K1022), 
            K1023 = sum(K1023), 
            K1024 = sum(K1024), 
            K1025 = sum(K1025), 
            K1026 = sum(K1026), 
            K1027 = sum(K1027), 
            K1028 = sum(K1028),
            K1029 = sum(K1029),
            records = sum(records)) 

length$n <- rowSums(length[, 2:11], na.rm=TRUE)

# aggregate to the county-level (across individual FRIS years)
length_i <- sum

# mask all columns calculated using <10 observations across all years
length_sub <- length %>% 
  mutate(K1020 = ifelse(records < 6, NA, K1020),
         K1021 = ifelse(records < 6, NA, K1021),
         K1022 = ifelse(records < 6, NA, K1022),
         K1023 = ifelse(records < 6, NA, K1023),
         K1024 = ifelse(records < 6, NA, K1024),
         K1025 = ifelse(records < 6, NA, K1025),
         K1026 = ifelse(records < 6, NA, K1026),
         K1027 = ifelse(records < 6, NA, K1027),
         K1028 = ifelse(records < 6, NA, K1028),
         K1029 = ifelse(records < 6, NA, K1029))

# mask all columns calculated using <10 observations across individual years
length_i_sub <- length_i %>% 
  mutate(K1020 = ifelse(records < 6, NA, K1020),
         K1021 = ifelse(records < 6, NA, K1021),
         K1022 = ifelse(records < 6, NA, K1022),
         K1023 = ifelse(records < 6, NA, K1023),
         K1024 = ifelse(records < 6, NA, K1024),
         K1025 = ifelse(records < 6, NA, K1025),
         K1026 = ifelse(records < 6, NA, K1026),
         K1027 = ifelse(records < 6, NA, K1027),
         K1028 = ifelse(records < 6, NA, K1028),
         K1029 = ifelse(records < 6, NA, K1029))

# All NAs result from having fewer than 6 unique STPOIDs in a county over the 4 FRIS years

# find the proportion distributed across FRIS years
proportion <- length_sub %>% 
  group_by(GEOID) %>% 
  summarise(K1020 = K1020/n*100, 
            K1021 = K1021/n*100, 
            K1022 = K1022/n*100, 
            K1023 = K1023/n*100, 
            K1024 = K1024/n*100, 
            K1025 = K1025/n*100, 
            K1026 = K1026/n*100, 
            K1027 = K1027/n*100, 
            K1028 = K1028/n*100, 
            K1029 = K1029/n*100)

# find the proportion distributed across individual FRIS years
proportion_i <- length_i_sub %>% 
  group_by(GEOID, YEAR) %>% 
  summarise(K1020 = K1020/n*100, 
            K1021 = K1021/n*100, 
            K1022 = K1022/n*100, 
            K1023 = K1023/n*100, 
            K1024 = K1024/n*100, 
            K1025 = K1025/n*100, 
            K1026 = K1026/n*100, 
            K1027 = K1027/n*100, 
            K1028 = K1028/n*100, 
            K1029 = K1029/n*100)

# Do all proportions add to 100%? YES!
proportion$total <- rowSums(proportion[, 2:11], na.rm=TRUE)
proportion_i$total <- rowSums(proportion_i[, 3:12], na.rm=TRUE)

# found proportion to two significant digits
proportion[,-1] <- round(proportion[,-1],2)
proportion_i[,-1] <- round(proportion_i[,-1],2)

# How many counties are represented each year if we do this thing yearly? ie. why are we aggregating based on the four years of data we have for visualizations?
sched_2003 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2003) # 212/408
sched_2008 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2008) # 198/408
sched_2013 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2013) # 258/408
sched_2018 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2018) # 290/408

# Let's build in the State so we can facet by State for the proportional bar chart!
state <- barriers_key %>% select(GEOID, COUNTY_ALPHA, STATE_ALPHA) %>% distinct()

# let's merge so we have all county-states represented
# we'll build all years AND 2018!
prop_merge_save <- merge(state, proportion, by = "GEOID", all = T)

### save this proportion table for export!
#saveRDS(prop_merge_save, "./out/DEER/sched-prop-allyrs.RDS")
write.csv(prop_merge_save, "./out/DEER/sched-prop-allyrs.csv")
# build individual years and row bind!
prop_merge_18 <- merge(state, sched_2018, by="GEOID", all=T)
prop_merge_13 <- merge(state, sched_2013, by="GEOID", all=T)
prop_merge_08 <- merge(state, sched_2008, by="GEOID", all=T)
prop_merge_03 <- merge(state, sched_2003, by="GEOID", all=T)

prop_merge_i_save <- rbind(prop_merge_03, prop_merge_08, prop_merge_13, prop_merge_18)

### save this proportion table for export!
#saveRDS(prop_merge_i_save, "./out/DEER/sched-prop-indyrs.RDS")
write.csv(prop_merge_i_save, "./out/DEER/sched-prop-indyrs.csv")
################################################################################################################################################# BUILD FACETTED PROPORTIONAL BAR CHARTS
####################################################################################################################
# build proportion table above into suitable data format for proportion barchart (ie. long format)
v <- prop_merge_save %>% select(!c("total")) %>% gather(category, proportion, 4:13)
w <- prop_merge_18 %>% select(!c("total")) %>% gather(category, proportion, 5:14)


scheduling_plots <- function(df, state, title, place) {
  
  x <- df %>% 
    filter(STATE_ALPHA == state) 

  scheduling_plot <- ggplot(x, aes(fill=category, y=proportion,x=COUNTY_ALPHA)) +
    geom_bar(position = "fill", stat = "identity") +
    scale_fill_brewer(palette = "Paired", name = "Category") +
    ggtitle(title) +
    theme_classic() +
    xlab("County") +
    ylab("Proportion") +
    coord_flip() 


  ggsave(place, height = 7, width = 7, dpi=300, units="in", device='png')
  
  return(scheduling_plot)

}

scheduling_plots(v, "ARIZONA", "Arizona", "./viz/scheduling/scheduling_AZ.png") 
scheduling_plots(v,"CALIFORNIA", "California", "./viz/scheduling/scheduling_CA.png") 
scheduling_plots(v, "COLORADO", "Colorado", "./viz/scheduling/scheduling_CO.png")
scheduling_plots(v, "IDAHO", "Idaho", "./viz/scheduling/scheduling_ID.png")
scheduling_plots(v, "MONTANA", "Montana", "./viz/scheduling/scheduling_MT.png")
scheduling_plots(v, "NEW MEXICO", "New Mexico", "./viz/scheduling/scheduling_NM.png")
scheduling_plots(v, "NEVADA", "Nevada", "./viz/scheduling/scheduling_NV.png") 
scheduling_plots(v, "OREGON", "Oregon", "./viz/scheduling/scheduling_OR.png") 
scheduling_plots(v, "UTAH", "Utah", "./viz/scheduling/scheduling_UT.png") 
scheduling_plots(v, "WASHINGTON", "Washington", "./viz/scheduling/scheduling_WA.png")
scheduling_plots(v, "WYOMING", "Wyoming", "./viz/scheduling/scheduling_WY.png")

scheduling_plots(w, "ARIZONA", "Arizona", "./viz/scheduling/scheduling_AZ18.png") 
scheduling_plots(w,"CALIFORNIA", "California", "./viz/scheduling/scheduling_CA18.png") 
scheduling_plots(w, "COLORADO", "Colorado", "./viz/scheduling/scheduling_CO18.png")
scheduling_plots(w, "IDAHO", "Idaho", "./viz/scheduling/scheduling_ID18.png")
scheduling_plots(w, "MONTANA", "Montana", "./viz/scheduling/scheduling_MT18.png")
scheduling_plots(w, "NEW MEXICO", "New Mexico", "./viz/scheduling/scheduling_NM18.png")
scheduling_plots(w, "NEVADA", "Nevada", "./viz/scheduling/scheduling_NV18.png") 
scheduling_plots(w, "OREGON", "Oregon", "./viz/scheduling/scheduling_OR18.png") 
scheduling_plots(w, "UTAH", "Utah", "./viz/scheduling/scheduling_UT18.png") 
scheduling_plots(w, "WASHINGTON", "Washington", "./viz/scheduling/scheduling_WA18.png")
scheduling_plots(w, "WYOMING", "Wyoming", "./viz/scheduling/scheduling_WY18.png")

################################################################################################################################################# BUILD MAP OF MOST IMPORTANT SCHEDULING METHOD
####################################################################################################################
# Select most important source by proportion for each county
find <- v %>% group_by(GEOID) %>% filter(proportion == max(proportion))
find_18 <- w %>% group_by(GEOID) %>% filter(proportion == max(proportion))

# how many counties are in each category on the maps?
find_summary <- find %>% 
  group_by(category) %>% 
  summarise(n = n())
find_summary

find_summary18 <- find_18 %>% 
  group_by(category) %>% 
  summarise(n = n())
find_summary18

# rename columns for visualizations
colnames(find)[4] <- "Scheduling Method"
colnames(find_18)[5] <- "Scheduling Method"

# build spatial using the western_ctys RDS
ctys <- readRDS("./data/western_ctys.RDS")

sched_sp <- sp::merge(ctys, find, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)
sched18_sp <- sp::merge(ctys, find_18, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)

# build maps
sched_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(sched_sp) + 
    tm_polygons("Scheduling Method", palette = c("#A6CEE3", "#1F78B4", "#B2DF8A", "#FB9A99", "#FDBF6F", "#FF7F00"), border.col = "white")  +
    tm_legend(outside = TRUE) +
    tm_layout(frame= FALSE)

tmap_save(sched_map, "./viz/scheduling/sched_map.png")

sched18_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(sched18_sp) + 
    tm_polygons("Scheduling Method", palette = c("#A6CEE3", "#1F78B4", "#B2DF8A", "#FB9A99", "#FDBF6F", "#FF7F00", "#6A3D9A"), border.col = "white")  +
    tm_legend(outside = TRUE) +
    tm_layout(frame= FALSE)

tmap_save(sched18_map, "./viz/scheduling/sched18_map.png")

```

# Sources of Information

**DATA IN**: sources-of-information.RDS & western_ctys.RDS
**DATA OUT**: info-prop-allyrs.csv & info-prop-indyrs.csv

Let's determine which Sources of Information are most important by county, facetted by State!

Code chunk K):

1) Ensure that all *key variables*: K1040-K1047 (2013 FRIS & 2018 IWMS) and K696 (2003 & 2008 FRIS) are weighted properly by their farm-level weight, statewij (responses from K695 were separated into K1040-K1047 to merge the four survey years in a different .Rmd).;
2) I summarise across all years and individual years, and all data is masked if *n* >= 6 unique STPOIDs was not satisfied.;
3) I build a proportional bar chart facetted by state where the proportional bar chart shows *only* the proportion of responses in each category across the four census years to give stakeholders a sense of how these vary across western counties. Variables are both categorical (information source cited) and continuous (proportions; raw numbers are NOT reported).;
4) I take the source of information with the highest proportion of responses and map those responses spatially across all years and in only 2018.

We will use the tables of proportions (again, NO RAW DATA) requested to DEER to discuss differences visualized in the proportional bar chart. 
```{r K) SOURCES OF INFORMATION-PROP, warning = FALSE, message = FALSE, eval=FALSE}
info_key <- readRDS("./data/sources-of-information.RDS")
library(viridis)

info_key[is.na(info_key)] <- 0
# create key column -- if 0, respondent did not answer to scheduling methods questions
info_key <- info_key %>% 
  mutate(empty = K1040+K1041+K1042+K1043+K1044+K1045+K1046+K1047) %>% 
  filter(!COUNTY_ALPHA == 0) #dunno why these are here, but they are not real GEOIDs

# create new df where if key above is 0, the row is deleted
info <- info_key %>% filter(empty != 0)

info <- info %>% 
  mutate(K1040_w = K1040*statewij, 
            K1041_w = K1041*statewij, 
            K1042_w = K1042*statewij, 
            K1043_w = K1043*statewij, 
            K1044_w = K1044*statewij, 
            K1045_w = K1045*statewij, 
            K1046_w = K1046*statewij, 
            K1047_w = K1047*statewij) 

# summarise (ensure that plyr is detached to summarise across groups)
sum_K1040 <- info %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1040_w)) %>% 
  summarise(K1040 = sum(K1040_w))
sum_K1041 <- info %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1041_w)) %>% 
  summarise(K1041 = sum(K1041_w))
sum <- merge(sum_K1040, sum_K1041, by = c("GEOID", "YEAR"), all=T)
sum_K1042 <- info %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1042_w)) %>% 
  summarise(K1042 = sum(K1042_w))
sum <- merge(sum, sum_K1042, by = c("GEOID", "YEAR"), all=T)
sum_K1043 <- info %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1043_w)) %>% 
  summarise(K1043 = sum(K1043_w))
sum <- merge(sum, sum_K1043, by = c("GEOID", "YEAR"), all=T)
sum_K1044 <- info %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1044_w)) %>% 
  summarise(K1044 = sum(K1044_w))
sum <- merge(sum, sum_K1044, by = c("GEOID", "YEAR"), all=T)
sum_K1045 <- info %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1045_w)) %>% 
  summarise(K1045 = sum(K1045_w))
sum <- merge(sum, sum_K1045, by = c("GEOID", "YEAR"), all=T)
sum_K1046 <- info %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1046_w)) %>% 
  summarise(K1046 = sum(K1046_w))
sum <- merge(sum, sum_K1046, by = c("GEOID", "YEAR"), all=T)
sum_K1047 <- info %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K1047_w)) %>% 
  summarise(K1047 = sum(K1047_w))
sum <- merge(sum, sum_K1047, by = c("GEOID", "YEAR"), all=T)


# Number of contributing records for each county
records <- info %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(GEOID)) %>% 
  summarise(records = length(STPOID))
sum <- merge(sum, records, by = c("GEOID", "YEAR"), all=T)

sum$n <- rowSums(sum[, 3:10], na.rm=TRUE)

# aggregate to the county-level (across FRIS years)
l <- sum
l[is.na(l)] <- 0
length <- l %>% 
  group_by(GEOID) %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047), 
            records = sum(records)) 

length$n <- rowSums(length[, 2:9], na.rm=TRUE)

# aggregate to the county-level (across individual FRIS years)
length_i <- sum

# mask all columns calculated using <10 observations across all years
length_sub <- length %>% 
  dplyr::mutate(K1040 = ifelse(records < 6, NA, K1040),
         K1041 = ifelse(records < 6, NA, K1041),
         K1042 = ifelse(records < 6, NA, K1042),
         K1043 = ifelse(records < 6, NA, K1043),
         K1044 = ifelse(records < 6, NA, K1044),
         K1045 = ifelse(records < 6, NA, K1045),
         K1046 = ifelse(records < 6, NA, K1046),
         K1047 = ifelse(records < 6, NA, K1047))

# mask all columns calculated using <10 observations across individual years
length_i_sub <- length_i %>% 
  dplyr::mutate(K1040 = ifelse(records < 6, NA, K1040),
         K1041 = ifelse(records < 6, NA, K1041),
         K1042 = ifelse(records < 6, NA, K1042),
         K1043 = ifelse(records < 6, NA, K1043),
         K1044 = ifelse(records < 6, NA, K1044),
         K1045 = ifelse(records < 6, NA, K1045),
         K1046 = ifelse(records < 6, NA, K1046),
         K1047 = ifelse(records < 6, NA, K1047))

# All NAs result from having fewer than 10 records in a county over the 4 FRIS years

proportion <- length_sub %>% 
  group_by(GEOID) %>% 
  summarise(K1040 = K1040/n*100, 
            K1041 = K1041/n*100, 
            K1042 = K1042/n*100, 
            K1043 = K1043/n*100, 
            K1044 = K1044/n*100, 
            K1045 = K1045/n*100, 
            K1046 = K1046/n*100, 
            K1047 = K1047/n*100)

proportion_i <- length_i_sub %>% 
  group_by(GEOID, YEAR) %>% 
  summarise(K1040 = K1040/n*100, 
            K1041 = K1041/n*100, 
            K1042 = K1042/n*100, 
            K1043 = K1043/n*100, 
            K1044 = K1044/n*100, 
            K1045 = K1045/n*100, 
            K1046 = K1046/n*100, 
            K1047 = K1047/n*100)

# Do all proportions add to 100%? YES!
proportion$total <- rowSums(proportion[, 2:9], na.rm=TRUE)
proportion_i$total <- rowSums(proportion_i[, 3:10], na.rm=TRUE)

# found proportion to two significant digits
proportion[,-1] <- round(proportion[,-1],2)
proportion_i[,-1] <- round(proportion_i[,-1],2)

# How many counties are represented each year if we do this thing yearly? ie. why are we aggregating based on the four years of data we have for visualizations?
info_2003 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2003) # 158/408
info_2008 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2008) # 129/408
info_2013 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2013) # 142/408
info_2018 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2018) # 149/408

# Let's build in the State so we can facet by State for the proportional bar chart!
state <- barriers_key %>% select(GEOID, COUNTY_ALPHA, STATE_ALPHA) %>% distinct()

# let's merge so we have all county-states represented
# we'll build all years AND 2018!
prop_merge_save <- merge(state, proportion, by = "GEOID", all = T)

### save this proportion table for export!
#saveRDS(prop_merge_save, "./out/DEER/info-prop-allyrs.RDS")
write.csv(prop_merge_save, "./out/DEER/info-prop-allyrs.csv")
# build individual years and row bind!
prop_merge_18 <- merge(state, info_2018, by="GEOID", all=T)
prop_merge_13 <- merge(state, info_2013, by="GEOID", all=T)
prop_merge_08 <- merge(state, info_2008, by="GEOID", all=T)
prop_merge_03 <- merge(state, info_2003, by="GEOID", all=T)

prop_merge_i_save <- rbind(prop_merge_03, prop_merge_08, prop_merge_13, prop_merge_18)

### save this proportion table for export!
#saveRDS(prop_merge_i_save, "./out/DEER/info-prop-indyrs.RDS")
write.csv(prop_merge_i_save, "./out/DEER/info-prop-indyrs.csv")
################################################################################################################################################# BUILD FACETTED PROPORTIONAL BAR CHARTS
####################################################################################################################
# build proportion table above into suitable data format for proportion barchart (ie. long format)
v <- prop_merge_save %>% select(!c("total")) %>% gather(category, proportion, 4:11)
w <- prop_merge_18 %>% select(!c("total")) %>% gather(category, proportion, 5:12)


information_plots <- function(df, state, title, place) {
  
  x <- df %>% 
    filter(STATE_ALPHA == state) 

  information_plot <- ggplot(x, aes(fill=category, y=proportion,x=COUNTY_ALPHA)) +
    geom_bar(position = "fill", stat = "identity") +
    scale_fill_brewer(palette = "Dark2", name = "Category") +
    ggtitle(title) +
    theme_classic() +
    xlab("County") +
    ylab("Proportion") +
    coord_flip() 


  ggsave(place, height = 7, width = 7, dpi=300, units="in", device='png')
  
  return(information_plot)

}

information_plots(v, "ARIZONA", "Arizona", "./viz/information/information_AZ.png") 
information_plots(v, "CALIFORNIA", "California", "./viz/information/information_CA.png") 
information_plots(v, "COLORADO", "Colorado", "./viz/information/information_CO.png")
information_plots(v, "IDAHO", "Idaho", "./viz/information/information_ID.png")
information_plots(v, "MONTANA", "Montana", "./viz/information/information_MT.png")
information_plots(v, "NEW MEXICO", "New Mexico", "./viz/information/information_NM.png")
information_plots(v, "NEVADA", "Nevada", "./viz/information/information_NV.png") 
information_plots(v, "OREGON", "Oregon", "./viz/information/information_OR.png")
information_plots(v, "UTAH", "Utah", "./viz/information/information_UT.png") 
information_plots(v, "WASHINGTON", "Washington", "./viz/information/information_WA.png")
information_plots(v, "WYOMING", "Wyoming", "./viz/information/information_WY.png")

information_plots(w, "ARIZONA", "Arizona", "./viz/information/information_AZ18.png") 
information_plots(w, "CALIFORNIA", "California", "./viz/information/information_CA18.png") 
information_plots(w, "COLORADO", "Colorado", "./viz/information/information_CO18.png")
information_plots(w, "IDAHO", "Idaho", "./viz/information/information_ID18.png")
information_plots(w, "MONTANA", "Montana", "./viz/information/information_MT18.png")
information_plots(w, "NEW MEXICO", "New Mexico", "./viz/information/information_NM18.png")
information_plots(w, "NEVADA", "Nevada", "./viz/information/information_NV18.png") 
information_plots(w, "OREGON", "Oregon", "./viz/information/information_OR18.png")
information_plots(w, "UTAH", "Utah", "./viz/information/information_UT18.png") 
information_plots(w, "WASHINGTON", "Washington", "./viz/information/information_WA18.png")
information_plots(w, "WYOMING", "Wyoming", "./viz/information/information_WY18.png")

################################################################################################################################################# BUILD MAP OF MOST IMPORTANT INFORMATION SOURCE
####################################################################################################################
# Select most important source by proportion for each county
find <- v %>% group_by(GEOID) %>% filter(proportion == max(proportion))
find_18 <- w %>% group_by(GEOID) %>% filter(proportion == max(proportion))

find_summary <- find %>% 
  group_by(category) %>% 
  summarise(n = n())
find_summary

find_summary18 <- find_18 %>% 
  group_by(category) %>% 
  summarise(n = n())
find_summary18

colnames(find)[4] <- "Information Source"
colnames(find_18)[5] <- "Information Source"
ctys <- readRDS("./data/western_ctys.RDS")

info_sp <- sp::merge(ctys, find, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)
info_sp18 <- sp::merge(ctys, find_18, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)

library(tmap)
info_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(info_sp) + 
    tm_polygons("Information Source", palette = "Dark2", border.col = "white")  +
    tm_legend(outside = TRUE) +
    tm_layout(frame= FALSE)

tmap_save(info_map, "./viz/information/info_map.png")

info_map18 <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(info_sp18) + 
    tm_polygons("Information Source", palette = "Dark2", border.col = "white")  +
    tm_legend(outside = TRUE) +
    tm_layout(frame= FALSE)

tmap_save(info_map18, "./viz/information/info_map18.png")

```

# Sources of Assistance

**DATA IN**: gov-fin-tech-assistance-wideformat.RDS & western_ctys.RDS
**DATA OUT**: assist-prop-allyrs.csv & assist-prop-indyrs.csv

Let's determine which Sources of Assistance are most important by county, facetted by State!

Code chunk L):

1) Ensure that all *key variables*: K215-K219 & K235-K239 (2013 FRIS & 2018 IWMS) are weighted properly by their farm-level weight, statewij (responses from 2003 and 2008 were not included in this analysis).;
2) I summarise across all years and individual years, and all data is masked if *n* >= 6 unique STPOIDs was not satisfied.;
3) I build a proportional bar chart facetted by state where the proportional bar chart shows *only* the proportion of responses in each category across the census years to give stakeholders a sense of how these vary across western counties. Variables are both categorical (source of assistance cited) and continuous (proportions; raw numbers are NOT reported).;
4) I take the source of assistance with the highest proportion of responses and map those responses spatially across all years.

We will use the tables of proportions (again, NO RAW DATA) requested to DEER to discuss differences visualized in the proportional bar chart. 
```{r L) SOURCES OF ASSISTANCe-PROP, warning=FALSE, message=FALSE, eval=FALSE}
assist_start <- readRDS("./data/gov-fin-tech-assistance-wideformat.RDS")
assist_key <- assist_start

assist_key[is.na(assist_key)] <- 0
# create key column -- if 0, respondent did not answer to scheduling methods questions
assist_key <- assist_key %>% 
  mutate(empty = K215+K216+K217+K218+K219+K235+K236+K237+K238+K239) %>% 
  filter(!COUNTY_ALPHA == 0) #dunno why these are here, but they are not real GEOIDs

# create new df where if key above is 0, the row is deleted
assist <- assist_key %>% filter(empty != 0)

assist <- assist %>% 
  mutate(K215_w = K215*statewij, 
            K216_w = K216*statewij, 
            K217_w = K217*statewij, 
            K218_w = K218*statewij, 
            K219_w = K219*statewij, 
            K235_w = K235*statewij, 
            K236_w = K236*statewij, 
            K237_w = K237*statewij,
            K238_w = K237*statewij,
            K239_w = K237*statewij) 

# summarise (ensure that plyr is detached to summarise across groups)
sum_K215 <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K215_w)) %>% 
  summarise(K215 = sum(K215_w))
sum_K216 <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K216_w)) %>% 
  summarise(K216 = sum(K216_w))
sum <- merge(sum_K215, sum_K216, by = c("GEOID", "YEAR"), all=T)
sum_K217 <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K217_w)) %>% 
  summarise(K217 = sum(K217_w))
sum <- merge(sum, sum_K217, by = c("GEOID", "YEAR"), all=T)
sum_K218 <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K218_w)) %>% 
  summarise(K218 = sum(K218_w))
sum <- merge(sum, sum_K218, by = c("GEOID", "YEAR"), all=T)
sum_K219 <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K219_w)) %>% 
  summarise(K219 = sum(K219_w))
sum <- merge(sum, sum_K219, by = c("GEOID", "YEAR"), all=T)
sum_K235 <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K235_w)) %>% 
  summarise(K235 = sum(K235_w))
sum <- merge(sum, sum_K235, by = c("GEOID", "YEAR"), all=T)
sum_K236 <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K236_w)) %>% 
  summarise(K236 = sum(K236_w))
sum <- merge(sum, sum_K236, by = c("GEOID", "YEAR"), all=T)
sum_K237 <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K237_w)) %>% 
  summarise(K237 = sum(K237_w))
sum <- merge(sum, sum_K237, by = c("GEOID", "YEAR"), all=T)
sum_K238 <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K238_w)) %>% 
  summarise(K238 = sum(K238_w))
sum <- merge(sum, sum_K238, by = c("GEOID", "YEAR"), all=T)
sum_K239 <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(K239_w)) %>% 
  summarise(K239 = sum(K239_w))
sum <- merge(sum, sum_K239, by = c("GEOID", "YEAR"), all=T)

# Number of contributing records for each county
records <- assist %>% 
  group_by(GEOID, YEAR) %>%
  filter(!is.na(GEOID)) %>% 
  summarise(records = length(STPOID))
sum <- merge(sum, records, by = c("GEOID", "YEAR"), all=T)

sum$n <- rowSums(sum[, 3:12], na.rm=TRUE)

# aggregate to the county-level (across FRIS years)
l <- sum
l[is.na(l)] <- 0
length <- l %>% 
  group_by(GEOID) %>% 
  summarise(K215 = sum(K215), 
            K216 = sum(K216), 
            K217 = sum(K217), 
            K218 = sum(K218), 
            K219 = sum(K219), 
            K235 = sum(K235), 
            K236 = sum(K236), 
            K237 = sum(K237),
            K238 = sum(K238),
            K239 = sum(K239),
            records = sum(records)) 

length$n <- rowSums(length[, 2:11], na.rm=TRUE)

# aggregate to the county-level (across individual FRIS years)
length_i <- sum

# mask all columns calculated using < 6 unique POIDs across all years
length_sub <- length %>% 
  dplyr::mutate(K215 = ifelse(records < 6, NA, K215),
         K216 = ifelse(records < 6, NA, K216),
         K217 = ifelse(records < 6, NA, K217),
         K218 = ifelse(records < 6, NA, K218),
         K219 = ifelse(records < 6, NA, K219),
         K235 = ifelse(records < 6, NA, K235),
         K236 = ifelse(records < 6, NA, K236),
         K237 = ifelse(records < 6, NA, K237),
         K238 = ifelse(records < 6, NA, K238),
         K239 = ifelse(records < 6, NA, K239))


# mask all columns calculated using <6 unique POIDs  across individual years
length_i_sub <- length_i %>% 
  dplyr::mutate(K215 = ifelse(records < 6, NA, K215),
         K216 = ifelse(records < 6, NA, K216),
         K217 = ifelse(records < 6, NA, K217),
         K218 = ifelse(records < 6, NA, K218),
         K219 = ifelse(records < 6, NA, K219),
         K235 = ifelse(records < 6, NA, K235),
         K236 = ifelse(records < 6, NA, K236),
         K237 = ifelse(records < 6, NA, K237),
         K238 = ifelse(records < 6, NA, K238),
         K239 = ifelse(records < 6, NA, K239))

# All NAs result from having fewer than 10 records in a county over the 4 FRIS years

proportion <- length_sub %>% 
  group_by(GEOID) %>% 
  summarise(K215 = K215/n*100, 
            K216 = K216/n*100, 
            K217 = K217/n*100, 
            K218 = K218/n*100, 
            K219 = K219/n*100, 
            K235 = K235/n*100, 
            K236 = K236/n*100, 
            K237 = K237/n*100,
            K238 = K238/n*100,
            K239 = K239/n*100)

proportion_i <- length_i_sub %>% 
  group_by(GEOID, YEAR) %>% 
  summarise(K215 = K215/n*100, 
            K216 = K216/n*100, 
            K217 = K217/n*100, 
            K218 = K218/n*100, 
            K219 = K219/n*100, 
            K235 = K235/n*100, 
            K236 = K236/n*100, 
            K237 = K237/n*100,
            K238 = K238/n*100,
            K239 = K239/n*100)

# Do all proportions add to 100%? YES!
proportion$total <- rowSums(proportion[, 2:11], na.rm=TRUE)
proportion_i$total <- rowSums(proportion_i[, 3:12], na.rm=TRUE)

# found proportion to two significant digits
proportion[,-1] <- round(proportion[,-1],2)
proportion_i[,-1] <- round(proportion_i[,-1],2)

# How many counties are represented each year if we do this thing yearly? ie. why are we aggregating based on the four years of data we have for visualizations?
assist_2013 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2013) # 21/408
assist_2018 <- proportion_i %>% filter(total != 0) %>% filter(YEAR == 2018) # 17/408

# Let's build in the State so we can facet by State for the proportional bar chart!
state <- barriers_key %>% select(GEOID, COUNTY_ALPHA, STATE_ALPHA) %>% distinct()

# let's merge so we have all county-states represented
# we'll build all years AND 2018!
prop_merge_save <- merge(state, proportion, by = "GEOID", all = T)

### save this proportion table for export!
#saveRDS(prop_merge_save, "./out/DEER/assist-prop-allyrs.RDS")
write.csv(prop_merge_save, "./out/DEER/assist-prop-allyrs.csv")
# build individual years and row bind!
prop_merge_18 <- merge(state, assist_2018, by="GEOID", all=T)
prop_merge_13 <- merge(state, assist_2013, by="GEOID", all=T)

prop_merge_i_save <- rbind(prop_merge_13, prop_merge_18)

### save this proportion table for export!
#saveRDS(prop_merge_i_save, "./out/DEER/assist-prop-indyrs.RDS")
write.csv(prop_merge_i_save, "./out/DEER/assist-prop-indyrs.csv")
################################################################################################################################################# BUILD FACETTED PROPORTIONAL BAR CHARTS
####################################################################################################################
# build proportion table above into suitable data format for proportion barchart (ie. long format)
v <- prop_merge_save %>% select(!c("total")) %>% gather(category, proportion, 4:13)

library(viridis)
assistance_plots <- function(df, state, title, place) {
  
  x <- df %>% 
    filter(STATE_ALPHA == state) 

  information_plot <- ggplot(x, aes(fill=category, y=proportion,x=COUNTY_ALPHA)) +
    geom_bar(position = "fill", stat = "identity") +
    scale_fill_viridis(discrete = TRUE, name = "Category") +
    ggtitle(title) +
    theme_classic() +
    xlab("County") +
    ylab("Proportion") +
    coord_flip() 


  ggsave(place, height = 7, width = 7, dpi=300, units="in", device='png')
  
  return(information_plot)

}

assistance_plots(v, "ARIZONA", "Arizona", "./viz/assistance/assistance_AZ.png") 
assistance_plots(v, "CALIFORNIA", "California", "./viz/assistance/assistance_CA.png") 
assistance_plots(v, "COLORADO", "Colorado", "./viz/assistance/assistance_CO.png")
assistance_plots(v, "IDAHO", "Idaho", "./viz/assistance/assistance_ID.png")
assistance_plots(v, "MONTANA", "Montana", "./viz/assistance/assistance_MT.png")
assistance_plots(v, "NEW MEXICO", "New Mexico", "./viz/assistance/assistance_NM.png")
assistance_plots(v, "NEVADA", "Nevada", "./viz/assistance/assistance_NV.png") 
assistance_plots(v, "OREGON", "Oregon", "./viz/assistance/assistance_OR.png")
assistance_plots(v, "UTAH", "Utah", "./viz/assistance/assistance_UT.png") 
assistance_plots(v, "WASHINGTON", "Washington", "./viz/assistance/assistance_WA.png")
assistance_plots(v, "WYOMING", "Wyoming", "./viz/assistance/assistance_WY.png")

################################################################################################################################################# BUILD MAP OF MOST IMPORTANT assistance SOURCE
####################################################################################################################
# Select most important source by proportion for each county
find <- v %>% group_by(GEOID) %>% filter(proportion == max(proportion))

find_summary <- find %>% 
  group_by(category) %>% 
  summarise(n = n())
find_summary


colnames(find)[4] <- "Type of Assistance"
ctys <- readRDS("./data/western_ctys.RDS")

assist_sp <- sp::merge(ctys, find, by = "GEOID", duplicateGeoms = TRUE, all = TRUE)

library(tmap)
assist_map <- tm_shape(ctys) + tm_polygons(col="grey", border.col = "white") +
  tm_shape(assist_sp) + 
    tm_polygons("Type of Assistance", palette = "viridis", border.col = "white")  +
    tm_legend(outside = TRUE) +
    tm_layout(frame= FALSE)

tmap_save(assist_map, "./viz/assistance/assist_map.png")

```

# Chi-Square

**DATA IN**: barriers-to-improvements.RDS, sources-of-information.RDS, gov-fin-tech-assistance-wideformat.RDS

Let's see how/if barriers, information, and assistance as associated--are they independent? OR, are the counts so different from what we would expect that we MUST conclude that there is an association?

*Key variables Code Chunk M)*: K1040-K1047 (2013 FRIS & 2018 IWMS) and K696 (2003 & 2008 FRIS), K1070-K1078 (2013 FRIS and 2018 IWMS) and K695 (2003 & 2008 FRIS), and K211 (2013 FRIS & 2018 IWMS, Y:N binary).

Code chunk M):

1) Pull in data built in other .Rmd's for barriers to improvements, sources of information, and government assistance.;
2) Create contingency tables of barriers (K1070-K1078) v. assistance (Y v. N), information source (K1040-K1047) v. assistance (Y v. N), and barriers (K1070-K1078) v. information source (K1040-K1047) in 2018.;
3) Build out row-wise and column-wise proportions of each and complete chi() (chi square) analysis.

```{r M) CHI SQUARE, warning = FALSE, message = FALSE}
b <- readRDS("./data/barriers-to-improvements.RDS")
i <- readRDS("./data/sources-of-information.RDS")
a <- readRDS("./data/gov-fin-tech-assistance-wideformat.RDS")
library(tidyverse)

###################################################################################################################
# barriers + assistance
###################################################################################################################
ba <- merge(b,a,by=c("GEOID", "YEAR", "STPOID", "STATE_ALPHA", "STATE_ABBR", "STATEFP", "COUNTY_ALPHA", "COUNTYFP", "statewij"))

# only interested in comparing when irrigator implemented improvements in 2018
ba <- ba %>% filter(IMPROVEMENTS == "YES") %>% filter(YEAR == 2018)
ba[is.na(ba)] <- 0

YES <- ba %>% 
  filter(ASSISTANCE == "YES") %>% 
  summarise(K1070 = sum(K1070), 
            K1071 = sum(K1071), 
            K1072 = sum(K1072), 
            K1073 = sum(K1073), 
            K1074 = sum(K1074), 
            K1075 = sum(K1075), 
            K1076 = sum(K1076), 
            K1077 = sum(K1077),
            K1078 = sum(K1078))

YES <- as.numeric(YES)

NO <- ba %>% 
  filter(ASSISTANCE == "NO") %>% 
  summarise(K1070 = sum(K1070), 
            K1071 = sum(K1071), 
            K1072 = sum(K1072), 
            K1073 = sum(K1073), 
            K1074 = sum(K1074), 
            K1075 = sum(K1075), 
            K1076 = sum(K1076), 
            K1077 = sum(K1077),
            K1078 = sum(K1078))

NO <- as.numeric(NO)

ba_contingency <- rbind(YES, NO)
colnames(ba_contingency) <- c("K1070", "K1071", "K1072", "K1073", "K1074", "K1075", "K1076", "K1077", "K1078")

# row-wise proportions
prop <- prop.table(ba_contingency, margin = 1)
prop

# column-wise proportions
prop2 <- prop.table(ba_contingency, margin = 2)
prop2

# chi square -- YES!--There IS an association between the variables!
ba_chi <- chisq.test(ba_contingency)
ba_chi

###################################################################################################################
# information + assistance
###################################################################################################################
ia <- merge(i,a,by=c("GEOID", "YEAR", "STPOID", "STATE_ALPHA", "STATE_ABBR", "STATEFP", "COUNTY_ALPHA", "COUNTYFP", "statewij"))

# only interested in comparing when irrigator implemented improvements in 2018
ia <- ia %>% filter(IMPROVEMENTS == "YES") %>% filter(YEAR == 2018)
ia[is.na(ia)] <- 0

YES <- ia %>% 
  filter(ASSISTANCE == "YES") %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

YES <- as.numeric(YES)

NO <- ia %>% 
  filter(ASSISTANCE == "NO") %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

NO <- as.numeric(NO)

ia_contingency <- rbind(YES, NO)
colnames(ia_contingency) <- c("K1040", "K1041", "K1042", "K1043", "K1044", "K1045", "K1046", "K1047")

#row-wise proportions
prop <- prop.table(ia_contingency, margin = 1)
prop

# column-wise proportions
prop2 <- prop.table(ia_contingency, margin = 2)
prop2

# chi square -- YES!--There is an association between source of information and assistance
ia_chi <- chisq.test(ia_contingency)
ia_chi

####################################################################################################################
# barriers + information
####################################################################################################################
bi <- merge(b,i,by=c("GEOID", "YEAR", "STPOID", "STATE_ALPHA", "STATE_ABBR", "STATEFP", "COUNTY_ALPHA", "COUNTYFP", "statewij"))

# only interested in understanding 2018
bi <- bi %>% filter(YEAR == 2018)

bi[is.na(bi)] <- 0

# investing in improvements is not a priority at this time
K1070 <- bi %>% 
  filter(K1070 == 1) %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

K1070 <- as.numeric(K1070)

# Risk of reduced yield or poorer quality crop from not meeting water needs
K1071 <- bi %>% 
  filter(K1071 ==1) %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

K1071 <- as.numeric(K1071)

# Physical field or crop conditions limit system improvements
K1072 <- bi %>% 
  filter(K1072 ==1) %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

K1072 <- as.numeric(K1072)

# Improvements will reduce costs but not enough to recover implementation costs
K1073 <- bi %>% 
  filter(K1073 ==1) %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

K1073 <- as.numeric(K1073)

# Cannot finance improvements
K1074 <- bi %>% 
  filter(K1074 ==1) %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

K1074 <- as.numeric(K1074)

# Landlord(s) will not share cost of improvements
K1075 <- bi %>% 
  filter(K1075 ==1) %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

K1075 <- as.numeric(K1075)

# Uncertainty about future availability of water
K1076 <- bi %>% 
  filter(K1076 ==1) %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

K1076 <- as.numeric(K1076)

#Will not be farming this operation long enough to justify new improvements
K1077 <- bi %>% 
  filter(K1077 ==1) %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

K1077 <- as.numeric(K1077)

# Improvements will increase management time or costs
K1078 <- bi %>% 
  filter(K1078 ==1) %>% 
  summarise(K1040 = sum(K1040), 
            K1041 = sum(K1041), 
            K1042 = sum(K1042), 
            K1043 = sum(K1043), 
            K1044 = sum(K1044), 
            K1045 = sum(K1045), 
            K1046 = sum(K1046), 
            K1047 = sum(K1047))

K1078 <- as.numeric(K1078)

bi_contingency <- rbind(K1070, K1071, K1072, K1073, K1074, K1075, K1076, K1077, K1078)
colnames(bi_contingency) <- c("K1040", "K1041", "K1042", "K1043", "K1044", "K1045", "K1046", "K1047")

# row-wise proportions
prop <- prop.table(bi_contingency, margin = 1)
prop

# column-wise proportions
prop2 <- prop.table(bi_contingency, margin = 2)
prop2

# chi square -- YES!--There is an association between barriers and information
bi_chi <- chisq.test(bi_contingency)
bi_chi
```

# ANOVAs

**DATA IN**: full-panel-FRIS.RDS & all_transformed.RDS

Aim 7, Q. Is irrigation productivity in 2018 significantly different by *X*?


*Key variables Code Chunk N)*: K51, K52, (CORN GRAIN); K61, K62 (CORN SILAGE); K81, K82 (WHEAT); K141, K142 (ALFALFA); K151, K152, (OTHER HAY). I then found irrigation productivity measured as unit yield/unit water applied. These variables represent:

K51: yield (I converted all yield to pounds)
K52: water use (I took all water use in acre feet/acre).

Code chunk N):

1) Inspect histograms at the individual-level and transform irrigation productivity accordingly for ANOVA analyses (must meet assumption of normality!). These analyses will be at the individual-level, so we will not be multiplying by 'statewij'. We are simply testing differences in associations across individual irrigators.

```{r N) LOG TRANSFORM, warning = FALSE, message = FALSE, eval = FALSE}

fris <- readRDS("./data/full-panel-FRIS.RDS")

# Okay, so this data is pretty darn skewed -- for any analyses (ANOVA, etc.) that assume normality, we will have to use a log() transformed version; full disclosure, I tried using square root and cube root and they do not normalize the data like we want. log() works beautifully, see below.
df_all <- fris %>% mutate(log_IP_AF = log(IP_AF))
hist(df_all$IP_AF)
hist(df_all$log_IP_AF)

saveRDS(df_all, "./out/all_transformed.RDS")

# now let's plot the four major crops to see what their distributions look like: alfalfa, corn silage/grain, wheat, other hay

distribution <- function(df, crop) {
  
  y <- df %>% 
    filter(CROP == crop)
  
  viz <- hist(y$IP_AF)
  
  viz2 <- hist(y$log_IP_AF)
    
  viz3 <- hist(y$YIELD_LBS)
  
  viz4 <- hist(y$WATER_AF)
  
  return(viz)
  
}

# everything looks good and normal-ish!
alfalfa <- distribution(df_all, "ALFALFA") # yield/water use skewed tail right
hay <- distribution(df_all, "OTHER_HAY") # yield/water use skewed tail right
wheat <- distribution(df_all, "WHEAT") # yield somewhat normal, water use skewed tail right
corngrain <- distribution(df_all, "CORN_GRAIN") # yield/water use somewhat normal
cornsilage <- distribution(df_all, "CORN_SILAGE") # yield/water use somewhat normal
```

*Key variables Code Chunk O)*: K51, K52, (CORN GRAIN); K61, K62 (CORN SILAGE); K81, K82 (WHEAT); K141, K142 (ALFALFA); K151, K152, (OTHER HAY). I then found irrigation productivity measured as unit yield/unit water applied. 

Code chunk O):

**DATA IN**: all_transformed.RDS

1) Now that we have transformed the data, it looks relatively normal. The log() transform of IP satisfies ANOVAs 'normality' assumption. Now, we have to test for equal variances, using the Levene's Test, where:

Null: The variances are equal (no significant difference in variances).
Alternative: Variances are NOT equal.

```{r O) LEVENEs TEST, warning = FALSE, message = FALSE}
library(car)
fris_transformed <- readRDS("./out/all_transformed.RDS")
f18_transformed <- fris_transformed %>% filter(YEAR==2018)

anova_data <- f18_transformed %>% filter(CROP == "ALFALFA" | CROP == "CORN_GRAIN" | CROP == "CORN_SILAGE" | CROP == "OTHER_HAY" | CROP == "WHEAT")

var_test <- leveneTest(log_IP_AF ~ CROP, data = anova_data)
var_test

# p < 0.05, reject the null hypothesis and conclude that variances are NOT equal

# Guideline: If the largest sample variance is less than 4x greater than the smallest sample variance, these are considered "close enough" to assume equal variances for a one-way ANOVA

crop_vars <- aggregate(log_IP_AF ~ CROP, data = anova_data, FUN = var)
crop_vars # largest variance is ~3.5x the smallest variance; shows that even though a formal hypothesis test reveals significant differences in varainces, the variances are not so different that we can't use an ANOVA
```

(1) Normal distribution? YES
(2) Independent? YES
(3) Continuous data? YES
(4) Equal variances? YES

### Crops

*Key variables*: K51, K52, (CORN GRAIN); K61, K62 (CORN SILAGE); K81, K82 (WHEAT); K141, K142 (ALFALFA); K151, K152, (OTHER HAY). I then found irrigation productivity measured as unit yield/unit water applied. 

Code chunk P):

1) Is IP different between crops? If yes, then test differences!

```{r P) IP-2018, warning = FALSE, message = FALSE}
crop_aov <- aov(log_IP_AF ~ CROP, data = anova_data)
summary(crop_aov) # mean ip is NOT the same across crops.

post_hoc <- TukeyHSD(crop_aov)
post_hoc
```

**Conclusion**: Irrigation productivity differed significantly by crop (one-way ANOVA, F(4) = 117.8, *p* < 0.001.). Post-hoc testing by Tukey's HSD revealed significant differences between irrigation productivity across all crops.

### Primary Irrigation Source

Because we do not have a balanced factorial design, we will use Type III ANOVAs in the car() package to asses two-way ANOVA.

*Key variables Code Chunk Q)*: K254, K255, K256 (CORN GRAIN); K264, K265, K266 (CORN SILAGE); K284, K285, K286 (WHEAT); K344, K345, K346 (ALFALFA); K354, K355, K356 (OTHER HAY). These variables represent:

K254: acres irrigated using on-farm surface water
K255: acres irrigated using well water
K256: acres irrigated using water from off-farm suppliers.

I summarised data for each irrigator-crop-year and created a categorical variable for an irrigators primary irrigation source (the water source on any operation that irrigated the greatest number of acres) in irrigation-productivity.Rmd.

Code chunks Q-W):

1) Build Type III ANOVA using the car() package to run a two-way ANOVA on log() trasnformed irrgiation productivity by crop and variable;
2) Run individual one-way ANOVAs on log_IP_AF for each crop-variable combination--when significant, run a post-hoc Tukey's HSD test to figure out where the differences actually are.;
3) Run summary statistics on all significantly different crops to show magnitude of differences.

```{r Q) PRIM_SOURCE, warning = FALSE, message = FALSE}
# select only where prim_source is a single category
source <- anova_data %>% filter(PRIM_SOURCE == "OFF_FARM" | PRIM_SOURCE == "ONFARM_SURFACE" | PRIM_SOURCE == "WELL_WATER")

crop_lm <- lm(log_IP_AF ~ CROP + PRIM_SOURCE + CROP*PRIM_SOURCE, data = source)
plot(crop_lm)
summary(crop_lm)

cropsource_aov <- Anova(crop_lm, type = "III")
cropsource_aov

# Let's do multiple one-way ANOVAs; ie. an ANOVA for each crop-source combination
# Alfalfa
alf_data <- source %>% filter(CROP == "ALFALFA")
alf_source_aov <- aov(log_IP_AF ~ PRIM_SOURCE, data = alf_data)
summary(alf_source_aov) # mean ip is the same across souces in Alfalfa

# Corn grain
cg_data <- source %>% filter(CROP == "CORN_GRAIN")
cg_source_aov <- aov(log_IP_AF ~ PRIM_SOURCE, data = cg_data)
summary(cg_source_aov) # mean ip is NOT the same across souces in Corn Grain

cg_source_post_hoc <- TukeyHSD(cg_source_aov)
cg_source_post_hoc

cg_source_ss <- cg_data %>% 
  group_by(PRIM_SOURCE) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

cg_source_ss

# Corn silage
cs_data <- source %>% filter(CROP == "CORN_SILAGE")
cs_source_aov <- aov(log_IP_AF ~ PRIM_SOURCE, data = cs_data)
summary(cs_source_aov) # mean ip is the same across souces in Corn Silage


# Hay
hay_data <- source %>% filter(CROP == "OTHER_HAY")
hay_source_aov <- aov(log_IP_AF ~ PRIM_SOURCE, data = hay_data)
summary(hay_source_aov) # mean ip is NOT the same across souces in HAy

hay_source_post_hoc <- TukeyHSD(hay_source_aov)
hay_source_post_hoc

hay_source_ss <- hay_data %>% 
  group_by(PRIM_SOURCE) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

hay_source_ss

# Wheat
wh_data <- source %>% filter(CROP == "WHEAT")
wh_source_aov <- aov(log_IP_AF ~ PRIM_SOURCE, data = wh_data)
summary(wh_source_aov) # mean ip is the same across souces in Wheat

```

**Conclusion**: Analysis by two-way ANOVA revealed that irrigation productivity differs significantly by crop (F(4) = 70.0, *p* < 0.001), but not primary source (F(2) = 0.38, *p* = 0.68), and there is a significant interaction between crop and primary source (F(8) = 5.6, *p* < 0.001).

**Conclusion**: Irrigation productivity differed significantly by water source in corn grain (one-way ANOVA, F(2) = 11.6, *p* < 0.001) and hay (one-way ANOVA, F(2) = 20.3, *p* < 0.001) in 2018. Post-hoc testing by Tukey's HSD revealed significant differences in mean irrigation productivity between well water (mean ip = 10522 lbs/af) and off farm water (mean ip = 6466 lbs/af) (*p* < 0.001) and well water and on farm surface water (mean ip = 6495.9 lbs/af) (*p* < 0.001) in hay, and between well water (mean ip = 7749 lbs/af) and off farm water (6195 lbs/af) in corn grain.

### Water Distribution Type

*Key Variables Code Chunk R)*: K250, K260, K280, K340, K350. These variables represent:

K250: primary method of water distribution (by crop).

I took the 24 potential codes reported by irrigators and created a binary response: Gravity (15-24) or Pressure (1-14) in irrigation-productivity.Rmd.

```{r R) H2O_DIST_BIN, warning = FALSE, message = FALSE}
# H2O_DIST_BIN
distbin_lm <- lm(log_IP_AF ~ CROP + H2O_DIST_BIN + CROP*H2O_DIST_BIN, data = anova_data)
plot(distbin_lm)
summary(distbin_lm)

cropdistbin_aov <- Anova(distbin_lm, type = "III")
cropdistbin_aov

distbin_ss <- anova_data %>% 
  group_by(H2O_DIST_BIN) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

distbin_ss

# Let's do multiple one-way ANOVAs; ie. an ANOVA for each crop-distribution combination
# Alfalfa
alf_data <- anova_data %>% filter(CROP == "ALFALFA")
alf_dist_aov <- aov(log_IP_AF ~ H2O_DIST_BIN, data = alf_data)
summary(alf_dist_aov) # mean ip is NOT the same across distribution types in Alfalfa

alf_dist_ss <- alf_data %>% 
  group_by(H2O_DIST_BIN) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

alf_dist_ss

# Corn grain
cg_data <- anova_data %>% filter(CROP == "CORN_GRAIN")
cg_dist_aov <- aov(log_IP_AF ~ H2O_DIST_BIN, data = cg_data)
summary(cg_dist_aov) # mean ip is NOT the same across distribution types in Corn Grain

cg_dist_ss <- cg_data %>% 
  group_by(H2O_DIST_BIN) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

cg_dist_ss

# Corn silage
cs_data <- anova_data %>% filter(CROP == "CORN_SILAGE")
cs_dist_aov <- aov(log_IP_AF ~ H2O_DIST_BIN, data = cs_data)
summary(cs_dist_aov) # mean ip is NOT the same across distribution types in Corn Silage

cs_dist_ss <- cs_data %>% 
  group_by(H2O_DIST_BIN) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

cs_dist_ss

# Hay
hay_data <- anova_data %>% filter(CROP == "OTHER_HAY")
hay_dist_aov <- aov(log_IP_AF ~ H2O_DIST_BIN, data = hay_data)
summary(hay_dist_aov) # mean ip is NOT the same across distribution types in HAy

hay_dist_ss <- hay_data %>% 
  group_by(H2O_DIST_BIN) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

hay_dist_ss

# Wheat
wh_data <- anova_data %>% filter(CROP == "WHEAT")
wh_dist_aov <- aov(log_IP_AF ~ H2O_DIST_BIN, data = wh_data)
summary(wh_dist_aov) # mean ip is NOT the same across souces in Wheat

wh_dist_ss <- wh_data %>% 
  group_by(H2O_DIST_BIN) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

wh_dist_ss
```

**Conclusion**: Analysis by two-way ANOVA revealed that irrigation productivity differs significantly by crop (F(4) = 50.8, *p* < 0.001), and water distribution type (F(1) = 25.5, *p* < 0.001), and there is a significant interaction between crop and water distribution type (F(4) = 4.4, *p* = 0.001).

**Conclusion**: Irrigation productivity differed significantly by water distribution type across all crops (*p* < 0.001) in 2018. For all crops, irrigation productivity in pressure systems (mean ip ~ 7864 lbs/af) is higher than in gravity systems (mean ip ~ 6947 lbs/af).

### Water Use Scheduling Method

*Key Variables Code Chunk S)*: K1020-K1029 (2013 FRIS & 2018 IWMS) and K527 (2003 & 2008 FRIS); (responses from K527 were separated into K1020-K1029 to merge the four survey years in a different .Rmd).

```{r S) WU_SCHED, warning = FALSE, message = FALSE}
# select only where WU_SCHED is not NA
sched <- anova_data %>% filter(!is.na(WU_SCHED))

# WU_SCHED
sched_lm <- lm(log_IP_AF ~ CROP + WU_SCHED, data = sched) # no interaction term because of NAs with 2 interactions
plot(sched_lm)
summary(sched_lm)

cropsched_aov <- Anova(sched_lm, type = "III")
cropsched_aov

# Let's do multiple one-way ANOVAs; ie. an ANOVA for each crop-scheduling combination
# Alfalfa
alf_data <- sched %>% filter(CROP == "ALFALFA")
alf_sched_aov <- aov(log_IP_AF ~ WU_SCHED, data = alf_data)
summary(alf_sched_aov) # mean ip is NOT the same across scheduling method in Alfalfa

alf_sched_post_hoc <- TukeyHSD(alf_sched_aov)
alf_sched_post_hoc

alf_sched_ss <- alf_data %>% 
  group_by(WU_SCHED) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

alf_sched_ss

# Corn grain
cg_data <- sched %>% filter(CROP == "CORN_GRAIN")
cg_sched_aov <- aov(log_IP_AF ~ WU_SCHED, data = cg_data)
summary(cg_sched_aov) # mean ip is the same across scheduling methods in Corn Grain


# Corn silage
cs_data <- sched %>% filter(CROP == "CORN_SILAGE")
cs_sched_aov <- aov(log_IP_AF ~ WU_SCHED, data = cs_data)
summary(cs_sched_aov) # mean ip is the same across scheduling methods in Corn Silage

# Hay
hay_data <- sched %>% filter(CROP == "OTHER_HAY")
hay_sched_aov <- aov(log_IP_AF ~ WU_SCHED, data = hay_data)
summary(hay_sched_aov) # mean ip is the same across scheduling methods in HAy

# Wheat
wh_data <- sched %>% filter(CROP == "WHEAT")
wh_sched_aov <- aov(log_IP_AF ~ WU_SCHED, data = wh_data)
summary(wh_sched_aov) # mean ip is the same across scheduling methods in Wheat

```
**Conclusion**: Analysis by two-way ANOVA revealed that irrigation productivity differs significantly by crop (F(4) = 41.7, *p* < 0.001), and primary water scheduling practice (F(7) = 2.5, *p* = 0.015).

**Conclusion**: Irrigation productivity differed significantly by scheduling method in alfalfa (one-way ANOVA, F(2) = 2.0, *p* = 0.05) but in no other crops in 2018. Post-hoc testing by Tukey's HSD revealed significant differences in mean irrigation productivity between when water was delivered (mean ip = 7332 lbs/af) and the feel of the soil (mean ip = 8939 lbs/af) (*p* = 0.04) in alfalfa.

### Barriers to Implementing Irrigation Improvements

*Key Variables Code Chunk T)*: K1070-K1078 (2013 FRIS and 2018 IWMS) and K695 (2003 & 2008 FRIS); (responses from K695 were separated into K1070-K1078 to merge the four survey years in a different .Rmd).

```{r T) BARRIER, warning = FALSE, message = FALSE}
# select only where BARRIER is not NA
barrier <- anova_data %>% filter(!is.na(BARRIER))

# BARRIER
barr_lm <- lm(log_IP_AF ~ CROP + BARRIER + CROP*BARRIER, data = barrier) 
plot(barr_lm)
summary(barr_lm)

cropbarr_aov <- Anova(barr_lm, type = "III")
cropbarr_aov

# Let's do multiple one-way ANOVAs; ie. an ANOVA for each crop-scheduling combination
# Alfalfa
alf_data <- barrier %>% filter(CROP == "ALFALFA")
alf_liv_aov <- aov(log_IP_AF ~ BARRIER, data = alf_data)
summary(alf_liv_aov) # mean ip is the same across barriers in Alfalfa

# Corn grain
cg_data <- barrier %>% filter(CROP == "CORN_GRAIN")
cg_barr_aov <- aov(log_IP_AF ~ BARRIER, data = cg_data)
summary(cg_barr_aov) # mean ip is the same across barriers in Corn Grain

cg_barr_post_hoc <- TukeyHSD(cg_barr_aov)
cg_barr_post_hoc

cg_barr_ss <- cg_data %>% 
  group_by(BARRIER) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

cg_barr_ss

# Corn silage
cs_data <- barrier %>% filter(CROP == "CORN_SILAGE")
cs_barr_aov <- aov(log_IP_AF ~ BARRIER, data = cs_data)
summary(cs_barr_aov) # mean ip is the same across barriers in Corn Silage

# Hay
hay_data <- barrier %>% filter(CROP == "OTHER_HAY")
hay_barr_aov <- aov(log_IP_AF ~ BARRIER, data = hay_data)
summary(hay_barr_aov) # mean ip is the same across barriers in HAy

# Wheat
wh_data <- barrier %>% filter(CROP == "WHEAT")
wh_barr_aov <- aov(log_IP_AF ~ BARRIER, data = wh_data)
summary(wh_barr_aov) # mean ip is the same across scheduling methods in Wheat
```

**Conclusion**: Analysis by two-way ANOVA revealed that irrigation productivity differs significantly by crop (F(4) = 6.1, *p* < 0.001), but not by primary barrier to irrigation improvements (F(8) = 1.3, *p* = 0.25). Additionally, there is no significant interaction between crop and primary barrier (F(32) = 0.76, *p* = 0.84). 

**Conclusion**: Irrigation productivity differed significantly by barrier in corn grain (one-way ANOVA, F(8) = 2.3, *p* = 0.03) but in no other crops in 2018. Post-hoc testing by Tukey's HSD revealed significant differences in mean irrigation productivity between etc.etc.etc.

### Sources of Information

*Key Variables Code Chunk U)*: K1040-K1047 (2013 FRIS & 2018 IWMS) and K696 (2003 & 2008 FRIS); (responses from K695 were separated into K1040-K1047 to merge the four survey years in a different .Rmd).

```{r U) INFO, warning = FALSE, message = FALSE}
# select only where INFO is not NA
info <- anova_data %>% filter(!is.na(INFO))

# INFORMATION
info_lm <- lm(log_IP_AF ~ CROP + INFO, data = info) # no interaction term because of NAs with 1 interactions
plot(info_lm)
summary(info_lm)

cropinfo_aov <- Anova(info_lm, type = "III")
cropinfo_aov

# Let's do multiple one-way ANOVAs; ie. an ANOVA for each crop-information combination
# Alfalfa
alf_data <- info %>% filter(CROP == "ALFALFA")
alf_info_aov <- aov(log_IP_AF ~ INFO, data = alf_data)
summary(alf_info_aov) # mean ip is the same across information type in Alfalfa

# Corn grain
cg_data <- info %>% filter(CROP == "CORN_GRAIN")
cg_info_aov <- aov(log_IP_AF ~ INFO, data = cg_data)
summary(cg_info_aov) # mean ip is NOT the same across information sources in Corn Grain

cg_info_post_hoc <- TukeyHSD(cg_info_aov)
cg_info_post_hoc

cg_info_ss <- cg_data %>% 
  group_by(INFO) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

cg_info_ss

# Corn silage
cs_data <- info %>% filter(CROP == "CORN_SILAGE")
cs_info_aov <- aov(log_IP_AF ~ INFO, data = cs_data)
summary(cs_info_aov) # mean ip is the same across information sources in Corn Silage

# Hay
hay_data <- info %>% filter(CROP == "OTHER_HAY")
hay_info_aov <- aov(log_IP_AF ~ INFO, data = hay_data)
summary(hay_info_aov) # mean ip is NOT the same across information sources in HAy

hay_info_post_hoc <- TukeyHSD(hay_info_aov)
hay_info_post_hoc

hay_info_ss <- hay_data %>% 
  group_by(INFO) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

hay_info_ss

# Wheat
wh_data <- info %>% filter(CROP == "WHEAT")
wh_info_aov <- aov(log_IP_AF ~ INFO, data = wh_data)
summary(wh_info_aov) # mean ip is the same across information sources in HAy
```

**Conclusion**: Analysis by two-way ANOVA revealed that irrigation productivity differs significantly by crop (F(4) = 19.5, *p* < 0.001), and by sources of information (F(7) = 2.8, *p* = 0.007). 

**Conclusion**: Irrigation productivity differed significantly by information source in corn grain (one-way ANOVA, F(2) = 2.8, *p* = 0.02) and hay (one-way ANOVA, F(7) = 2.02, *p* = 0.05) in 2018. Post-hoc testing by Tukey's HSD revealed significant differences in mean irrigation productivity between information provided by private irrigation specialists or crop consultants (mean ip = 8272 lbs/af) and irrigation equipment dealers (mean ip = 5567 lbs/af) (*p* < 0.05) in corn grain, and private irrigation specialists or crop consultants (mean ip = 10234 lbs/af) and local irrigation employees (mean ip = 4847 lbs/af) in hay.

### Livestock

*Key Variables Code Chunk V)*: K982 changed from "What percent of total sales were from livestock sales?" in 2003 to "What percent of total sales were from non-irigated crop sales and/or livestock sales?" in 2008/2013/2018. Coded as a binary indicator in value-crop-sales.Rmd assuming that IF irrigators answered any percentage in K982, some part of their operation had livestock.

```{r V) LIVESTOCK, warning = FALSE, message = FALSE}
# select only where LIVESTOCK is not NA
livestock <- anova_data %>% filter(!is.na(LIVESTOCK))
livestock$LIVESTOCK <- ifelse(livestock$LIVESTOCK == 1, "Yes", "No")

# LIVESTOCK
livestock_lm <- lm(log_IP_AF ~ CROP + LIVESTOCK + CROP*LIVESTOCK, data = livestock) 
plot(livestock_lm)
summary(livestock_lm)

croplive_aov <- Anova(livestock_lm, type = "III")
croplive_aov

# Let's do multiple one-way ANOVAs; ie. an ANOVA for each crop-scheduling combination
# Alfalfa
alf_data <- livestock %>% filter(CROP == "ALFALFA")
alf_liv_aov <- aov(log_IP_AF ~ LIVESTOCK, data = alf_data)
summary(alf_liv_aov) # mean ip is the different across livestock holding in Alfalfa

alf_liv_post_hoc <- TukeyHSD(alf_liv_aov)
alf_liv_post_hoc

alf_liv_ss <- alf_data %>% 
  dplyr::group_by(LIVESTOCK) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

alf_liv_ss

# Corn grain
cg_data <- livestock %>% filter(CROP == "CORN_GRAIN")
cg_liv_aov <- aov(log_IP_AF ~ LIVESTOCK, data = cg_data)
summary(cg_liv_aov) # mean ip is the different across livestock holding in Corn Grain

cg_liv_post_hoc <- TukeyHSD(cg_liv_aov)
cg_liv_post_hoc

cg_liv_ss <- cg_data %>% 
  dplyr::group_by(LIVESTOCK) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

cg_liv_ss

# Corn silage
cs_data <- barrier %>% filter(CROP == "CORN_SILAGE")
cs_liv_aov <- aov(log_IP_AF ~ LIVESTOCK, data = cs_data)
summary(cs_liv_aov) # mean ip is the same across livestock holding in Corn Silage

# Hay
hay_data <- livestock %>% filter(CROP == "OTHER_HAY")
hay_liv_aov <- aov(log_IP_AF ~ LIVESTOCK, data = hay_data)
summary(hay_liv_aov) # mean ip is the different across livestock holding in HAy

hay_liv_post_hoc <- TukeyHSD(hay_liv_aov)
hay_liv_post_hoc

hay_liv_ss <- hay_data %>% 
  dplyr::group_by(LIVESTOCK) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

hay_liv_ss

# Wheat
wh_data <- livestock %>% filter(CROP == "WHEAT")
wh_liv_aov <- aov(log_IP_AF ~ LIVESTOCK, data = wh_data)
summary(wh_liv_aov) # mean ip is the same across livestock holding in Wheat

wh_liv_ss <- wh_data %>% 
  dplyr::group_by(LIVESTOCK) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

wh_liv_ss
```

**Conclusion**: Analysis by two-way ANOVA revealed that irrigation productivity differs significantly by crop (F(4) = 33.0, *p* < 0.001), and by livestock holding (F(1) = 5.4, *p* = 0.02). Additionally, there is a significant interaction between crop and livestock holding (F(4) = 9.9, *p* < 0.001). 

**Conclusion**: Irrigation productivity differed significantly by livestock ownership in alfalfa (one-way ANOVA, *p* = 0.03), corn grain (one-way ANOVA, *p* < 0.001), wheat (one-way ANOVA, *p* = 0.02) and hay (one-way ANOVA, *p* = 0.006) in 2018. Livestock owners had a mean irrigation productivity (mean ip = 7558 lbs/af) over 700 lbs/acre higher than irrigators with no livestock  (mean ip = 6818 lbs/af) in alfalfa, and over 1350 lbs/acre higher (mean ip = 7798 v. 6438) in corn grain. In hay, irrigators with livestock had an average irrigation productivity over 600 lbs/acre less (mean ip = 7080 v. 7725) than irrigators with no livestock. In wheat, irrigators with livestock had an average irrigation productivity over 600 lbs/acre more (mean ip = 4961 v. 4329) than irrigators with no livestock.


### Range of Gross Sales

*Key Variables Code Chunk W)*: K980. 

In the 2018 IWMS:
         K980-1 = $0-9,999;
         K980-2 = $10k-24,999;
         K980-3 = $25k-49,999;
         K980-4 = $50k-99,999;
         K980-5 = $100k-249,999;
         K980-6 = $250k-499,999;
         K980-7 = $500k-999,999;
         K980-8 = > $1 million.

```{r W) GROSS SALES, warning = FALSE, message =  FALSE}
# select only where SALES is not NA
sales <- anova_data %>% filter(!is.na(SALES))

# SALES
sales_lm <- lm(log_IP_AF ~ CROP + SALES + CROP*SALES, data = sales) 
plot(sales_lm)
summary(sales_lm)

cropsales_aov <- Anova(sales_lm, type = "III")
cropsales_aov

# Let's do multiple one-way ANOVAs; ie. an ANOVA for each crop-scheduling combination
# Alfalfa
alf_data <- sales %>% filter(CROP == "ALFALFA")
alf_sales_aov <- aov(log_IP_AF ~ SALES, data = alf_data)
summary(alf_sales_aov) # mean ip is the same across sales holding in Alfalfa

# Corn grain
cg_data <- sales %>% filter(CROP == "CORN_GRAIN")
cg_sales_aov <- aov(log_IP_AF ~ SALES, data = cg_data)
summary(cg_sales_aov) # mean ip is the different across sales holding in Corn Grain

# Corn silage
cs_data <- barrier %>% filter(CROP == "CORN_SILAGE")
cs_sales_aov <- aov(log_IP_AF ~ SALES, data = cs_data)
summary(cs_sales_aov) # mean ip is the same across sales holding in Corn Silage

# Hay
hay_data <- sales %>% filter(CROP == "OTHER_HAY")
hay_sales_aov <- aov(log_IP_AF ~ SALES, data = hay_data)
summary(hay_sales_aov) # mean ip is the different across sales holding in HAy

hay_sales_post_hoc <- TukeyHSD(hay_sales_aov)
hay_sales_post_hoc

hay_sales_ss <- hay_data %>% 
  dplyr::group_by(SALES) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

hay_sales_ss

# Wheat
wh_data <- sales %>% filter(CROP == "WHEAT")
wh_sales_aov <- aov(log_IP_AF ~ SALES, data = wh_data)
summary(wh_sales_aov) # mean ip is the same across sales holding in Wheat

wheat_sales_post_hoc <- TukeyHSD(wh_sales_aov)
wheat_sales_post_hoc

wheat_sales_ss <- wh_data %>% 
  dplyr::group_by(SALES) %>% 
  summarise(mean = mean(IP_AF),
            sd = sd(IP_AF),
            med = median(IP_AF))

wheat_sales_ss
```

**Conclusion**: Analysis by two-way ANOVA revealed that irrigation productivity differs significantly by crop (F(4) = 9.2, *p* < 0.001), but not by sales (F(7) = 1.3, *p* = 0.23). Additionally, there is not a significant interaction between crop and sales (F(28) = 1.4, *p* = 0.06). 

**Conclusion**: Irrigation productivity differed significantly by gross sales bin in hay (one-way ANOVA, f(7) = 2.3, *p* = 0.02) and wheat (one-way ANOVA, F(7) = 3.5, *p* < 0.001) in 2018. Post-hoc testing by Tukey's HSD revealed significant differences in mean irrigation productivity between farmers bringing in 500K-99,999 in gross sales (mean ip = 8409 lbs/af) and $100k-249,999 in gross sales (mean ip = 5324 lbs/af) (*p* < 0.05) in hay, 50k-99,999 and 10k-24,999 (mean ip = 10410 lbs/af), greater than 1m and 10k-24,999, etc. etc. etc. in hay.

# Correlation plots

**DATA IN**: all_transformed.RDS

Aim 6, Q: Is irrigation productivity across major crops related to land ownership, farm size, or irrigated proportion in 2018?

*Key Variables Code Chunk X)*: K51, K52, (CORN GRAIN); K61, K62 (CORN SILAGE); K81, K82 (WHEAT); K141, K142 (ALFALFA); K151, K152, (OTHER HAY). I then found irrigation productivity measured as unit yield/unit water applied. K25 = land owned, K26 = land rented or leased FROM others, K27 = land rented or leased TO others, K28 = total acres in this operation, K41 = total acres on "this operation", K42 = total irrigated acres

Code chunk X):

1) Bring in the transformed data (this data is log(RAW) for irrigation productivity and raw for all other variables). Build correlations and visualize correlations (removing outliers >95th percentile).

^ very little correlation between any of the variables we chose to look at and irrigation productivity.
```{r X) CORRELATIONS, warning = FALSE, message = FALSE}
library(plyr)

corr_data <- readRDS("./out/all_transformed.RDS")

corr_data <- corr_data %>% 
    filter(YEAR == 2018) %>% 
    select(CROP, IP_AF, IP_IN, log_IP_AF, own, rent_to, rent_from, total_acres, irrig_acres, prop_own, prop_ir) 

corr_data$total_acres <- as.numeric(corr_data$total_acres)
corr_data$irrig_acres <- as.numeric(corr_data$irrig_acres)

#####################################################################################################################
### ownership
#####################################################################################################################
require(plyr)
own_data <- corr_data %>% filter(!is.na(own))
corr_own <- function(own_data) {
  return(data.frame(CORR = cor(own_data$log_IP_AF, own_data$own)))
}

corr_own_table <- ddply(own_data, .(CROP), corr_own)
corr_own_table

#####################################################################################################################
### rent from
#####################################################################################################################
rent_data <- corr_data %>% filter(!is.na(rent_from))
corr_rent <- function(rent_data) {
  return(data.frame(CORR = cor(rent_data$log_IP_AF, rent_data$rent_from)))
}

corr_rent_table <- ddply(rent_data, .(CROP), corr_rent)
corr_rent_table

#####################################################################################################################
### rent to
#####################################################################################################################
rent2_data <- corr_data %>% filter(!is.na(rent_to))
corr_rent2 <- function(rent2_data) {
  return(data.frame(CORR = cor(rent2_data$log_IP_AF, rent2_data$rent_to)))
}

corr_rent2_table <- ddply(rent2_data, .(CROP), corr_rent2)
corr_rent2_table

#####################################################################################################################
### proportion owned
#####################################################################################################################
propown_data <- corr_data %>% filter(!is.na(prop_own))
corr_propown <- function(propown_data) {
  return(data.frame(CORR = cor(propown_data$log_IP_AF, propown_data$prop_own)))
}

corr_propown_table <- ddply(propown_data, .(CROP), corr_propown)
corr_propown_table

#####################################################################################################################
### irrigated acres
#####################################################################################################################
irrigacres_data <- corr_data %>% filter(!is.na(irrig_acres))
corr_irrigacres <- function(irrigacres_data) {
  return(data.frame(CORR = cor(irrigacres_data$log_IP_AF, irrigacres_data$irrig_acres)))
}

corr_irrigacres_table <- ddply(irrigacres_data, .(CROP), corr_irrigacres)
corr_irrigacres_table

#####################################################################################################################
### proportion irrigated
#####################################################################################################################
propirr_data <- corr_data %>% filter(!is.na(prop_ir))
corr_propirr <- function(propirr_data) {
  return(data.frame(CORR = cor(propirr_data$log_IP_AF, propirr_data$prop_ir)))
}

corr_propirr_table <- ddply(propirr_data, .(CROP), corr_propirr)
corr_propirr_table

```

# Box and Whisker Plots 

**DATA IN**: all_transformed.RDS

### Box and Whisker Plots of Irrigation Productivity by crop through time

*Key Variables Code Chunk Y)*: K51, K52, (CORN GRAIN); K61, K62 (CORN SILAGE); K81, K82 (WHEAT); K141, K142 (ALFALFA); K151, K152, (OTHER HAY). I then found irrigation productivity measured as unit yield/unit water applied.

Code chunk Y):

1) Build box and whisker plots of irrigation productivity by crop through time using the raw data (not transformed because not summarised to a county or state level). Ensure that all state-crop-year combinations have *n* >= 30 observations across which the box and whisker plots are being made. I have also suppressed all outliers and data < 0.95 percentile (removes maximum state-crop-year combination in all cases).

```{r Y) BOX AND WHISKER IP, warning = FALSE, message = FALSE, eval = FALSE}
detach("package:spdplyr")
detach("package:dplyr")
library(dplyr)
# bring in data
bw_data <- readRDS("./out/all_transformed.RDS")
write.csv(bw_data, "./out/DEER/transformed_rawdata_NOTFOREXPORT.csv")
bw_data <- bw_data %>% filter(CROP == "ALFALFA" | CROP == "WHEAT" | CROP == "OTHER_HAY" | CROP == "CORN_SILAGE" | CROP == "CORN_GRAIN")

bw_data$YEAR <- as.factor(bw_data$YEAR)

# rename elements for legend
bw_data$CROP[bw_data$CROP == "ALFALFA"] <- "Alfalfa"
bw_data$CROP[bw_data$CROP == "CORN_GRAIN"] <- "Corn Grain"
bw_data$CROP[bw_data$CROP == "CORN_SILAGE"] <- "Corn Silage"
bw_data$CROP[bw_data$CROP == "OTHER_HAY"] <- "Other Hay"
bw_data$CROP[bw_data$CROP == "WHEAT"] <- "Wheat"


# ensure at n >= 30 unique POIDs per state-crop-year
bw_data <- bw_data %>% 
    group_by(STATE_ALPHA, CROP, YEAR) %>% 
    filter(IP_AF < quantile(IP_AF, 0.95))

bw_data_fix <- bw_data %>% 
    group_by(STATE_ALPHA, CROP, YEAR) %>% 
    filter(IP_AF < quantile(IP_AF, 0.95))%>% 
    dplyr::summarise(unique = length(unique((STPOID)))) %>% 
    filter(unique >= 30)

bw_data_fix <- merge(bw_data, bw_data_fix, by = c("STATE_ALPHA", "CROP", "YEAR"))

write.csv(bw_data_fix, "./out/DEER/boxwhisker_rawdata_NOTFOREXPORT.csv")

# build factor levels
bw_data_fix$CROP <- factor(bw_data_fix$CROP, levels = c("Alfalfa", "Corn Grain", "Corn Silage", "Other Hay", "Wheat"))
colors=c("#74C476", "#6BAED6", "#9E9AC8", "#FB6A4A", "#FD8D3C")
names(colors) <- levels(bw_data_fix$CROP)

#############################################################################################################################
# facetted by state! Irrigation productivity (lbs/AF)
#############################################################################################################################

bw <- ggplot(bw_data_fix, aes(x = YEAR, y = IP_AF, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Irrigation Productivity (pounds/acre foot)") +
  scale_fill_manual(name ="CROP",  values=colors) +
  coord_cartesian(ylim = c(0,30000)) +
  facet_wrap(~ STATE_ALPHA)

bw
ggsave(filename = "./viz/box_whisker/ip_bw_lbsAF.png", plot = bw, dpi = 300, width = 7, height = 5)

#############################################################################################################################
# Individual states! Irrigation productivity (lbs/AF)
#############################################################################################################################

bw_func <- function(state,top) {
  
        df <- bw_data_fix %>% 
          filter(STATE_ALPHA == state)
        
bw <- ggplot(df, aes(x = YEAR, y = IP_AF, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Irrigation Productivity (pounds/acre foot)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,top)) 

bw
ggsave(paste0(filename = './viz/box_whisker/ip_bw_',state,'_lbsAF.png'), plot = bw, dpi = 300, width = 7, height = 5)

}

bw_func("ARIZONA",8000)
bw_func("CALIFORNIA",16000)
bw_func("COLORADO",25000)
bw_func("IDAHO",22000)
bw_func("NEVADA",8000)
bw_func("NEW MEXICO",40000)
bw_func("MONTANA",30000)
bw_func("OREGON",28000)
bw_func("UTAH",23000)
bw_func("WASHINGTON",18000)
bw_func("WYOMING",20000)

#############################################################################################################################
# facetted by state! Irrigation productivity (kg/m)
#############################################################################################################################

bw_M <- ggplot(bw_data_fix, aes(x = YEAR, y = IP_KGM, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Irrigation Productivity (kilograms/cubic meter)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,10)) +
  facet_wrap(~ STATE_ALPHA)

bw_M
ggsave(filename = "./viz/box_whisker/ip_bw_kgM.png", plot = bw_M, dpi = 300, width = 7, height = 5)

#############################################################################################################################
# facetted by state! Irrigation productivity (kg/m)
#############################################################################################################################

bw_func <- function(state,top) {
  
        df <- bw_data_fix %>% 
          filter(STATE_ALPHA == state)
        
bw <- ggplot(df, aes(x = YEAR, y = IP_KGM, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Irrigation Productivity (kilograms/cubic meter)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,top)) 

bw
ggsave(paste0(filename = './viz/box_whisker/ip_bw_',state,'_kgM.png'), plot = bw, dpi = 300, width = 7, height = 5)

}

bw_func("ARIZONA",3)
bw_func("CALIFORNIA",5.5)
bw_func("COLORADO",9)
bw_func("IDAHO",8)
bw_func("NEVADA",3)
bw_func("NEW MEXICO",17)
bw_func("MONTANA",11)
bw_func("OREGON",10)
bw_func("UTAH",9)
bw_func("WASHINGTON",7)
bw_func("WYOMING",7)

#############################################################################################################################
# facetted by crop! Irrigation productivity (kg/m)
#############################################################################################################################

bw_func <- function(crop,top) {
  
        df <- bw_data_fix %>% 
          filter(CROP == crop)
        
bw <- ggplot(df, aes(x = YEAR, y = IP_KGM, fill = crop)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Irrigation Productivity (kilograms/cubic meter)") +
  scale_fill_manual(name = "CROP", values=colors) +
  facet_wrap(~ STATE_ALPHA) +
  coord_cartesian(ylim = c(0,top)) 

bw
ggsave(paste0(filename = './viz/box_whisker/ip_bw_',crop,'_kgM.png'), plot = bw, dpi = 300, width = 7, height = 5)

}

bw_func("Alfalfa",8)
bw_func("Corn Grain",5)
bw_func("Corn Silage",11)
bw_func("Other Hay",17)
bw_func("Wheat",4.5)

#############################################################################################################################
# let's pull some statistics out of these boxplots!
#############################################################################################################################
sos <- as.data.frame(ggplot_build(bw_M)$data[[1]])
sos <- sos[c(1,3:5,12:13)]

sos$CROP <- ifelse(sos$fill == "#74C476", "Alfalfa",
                   ifelse(sos$fill == "#6BAED6", "Corn Grain",
                        ifelse(sos$fill == "#9E9AC8", "Corn Silage",
                            ifelse(sos$fill == "#FB6A4A", "Other Hay",
                                ifelse(sos$fill == "#FD8D3C", "Wheat", 
                   "NA")))))
sos$YEAR <- ifelse(sos$group %in% c(1,2,3,4,5), 2003,
                   ifelse(sos$group %in% c(6,7,8,9,10), 2008,
                          ifelse(sos$group %in% c(11,12,13,14,15), 2013,
                                 ifelse(sos$group %in% c(16,17,18,19,20), 2018,
                                        "NA"))))

#by crop-state
sd_sos <- sos %>% group_by(CROP, PANEL) %>% summarise(median = median(middle), sd = sd(middle))

sos_merge <- merge(sos,sd_sos, by = c("CROP", "PANEL"), all = T)

sos_merge <- sos_merge %>% mutate(variable = "IP_M") %>% select(-fill)

# by crop across panel
median_sos <- sos %>% group_by(CROP) %>% summarise(median = median(middle), sd = sd(middle)) %>% mutate(variable = "IP_M")

bw_all <- ggplot(bw_data_fix, aes(x = YEAR, y = IP_KGM, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Irrigation Productivity (kilograms/cubic meter)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,8.5)) 

bw_all
ggsave(filename = "./viz//box_Whisker/ip_bw_all_kgM.png", plot = bw_all, dpi = 300, width = 7, height = 5)

# by crop-year
cy <- as.data.frame(ggplot_build(bw_all)$data[[1]])
cy <- cy[c(1,3:5,12:13)]

cy$CROP <- ifelse(cy$fill == "#74C476", "Alfalfa",
                   ifelse(cy$fill == "#6BAED6", "Corn Grain",
                        ifelse(cy$fill == "#9E9AC8", "Corn Silage",
                            ifelse(cy$fill == "#FB6A4A", "Other Hay",
                                ifelse(cy$fill == "#FD8D3C", "Wheat", 
                   "NA")))))
cy$YEAR <- ifelse(cy$group %in% c(1,2,3,4,5), 2003,
                   ifelse(cy$group %in% c(6,7,8,9,10), 2008,
                          ifelse(cy$group %in% c(11,12,13,14,15), 2013,
                                 ifelse(cy$group %in% c(16,17,18,19,20), 2018,
                                        "NA"))))

x_cy <- cy %>% group_by(CROP) %>% summarise(median = median(middle), sd = sd(middle)) %>% mutate(variable = "IP_M")

cy_merge <- merge(cy,x_cy, by = c("CROP"), all = T)

```

### Box and Whisker plot of water applied by crop through time

*Key Variables Code Chunk Z)*: K52 (CORN GRAIN); K62 (CORN SILAGE); K82 (WHEAT); K142 (ALFALFA); K152, (OTHER HAY). 

Code chunk Z):

1) Build box and whisker plots of water use by crop and through time using the raw data (not transformed because not summarised to a county or state level). I have also suppressed all outliers and data < 0.95 percentile (removes maximum state-crop-year combination in all cases).

```{r Z) BOX AND WHISKER WU, warning = FALSE, message = FALSE, eval = FALSE}
# use data from code chunk above
#############################################################################################################################
# facetted by state! Water Use (AF)
#############################################################################################################################

bw_wa <- ggplot(bw_data_fix, aes(x = YEAR, y = WATER_AF, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) +
  theme(legend.position = "bottom") +
  ylab("Water Applied (acre feet/acre)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,11.5)) + 
  facet_wrap(~ STATE_ALPHA)

bw_wa
ggsave(filename = "./viz/box_whisker/wa_bw_AF.png", plot = bw_wa, dpi = 300, width = 7, height = 5)

#############################################################################################################################
# Individual states! Water use (AF)
#############################################################################################################################

bw_func <- function(state,top) {
  
        df <- bw_data_fix %>% 
          filter(STATE_ALPHA == state)
        
bw <- ggplot(df, aes(x = YEAR, y = WATER_AF, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Water Applied (acre feet/acre)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,top)) 

bw
ggsave(paste0(filename = './viz/box_whisker/wa_bw_',state,'_AF.png'), plot = bw, dpi = 300, width = 7, height = 5)

}

bw_func("ARIZONA",12)
bw_func("CALIFORNIA",8)
bw_func("COLORADO",5)
bw_func("IDAHO",6)
bw_func("NEVADA",7.5)
bw_func("NEW MEXICO",6)
bw_func("MONTANA",4.5)
bw_func("OREGON",5)
bw_func("UTAH",6)
bw_func("WASHINGTON",5.5)
bw_func("WYOMING",5.5)

#############################################################################################################################
# facetted by state! Water use (cubic meters)
#############################################################################################################################

bw_wa_m <- ggplot(bw_data_fix, aes(x = YEAR, y = WATER_M, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) +
  theme(legend.position = "bottom") +
  ylab("Water Applied (cubic meters/hectare)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,35000)) + 
  facet_wrap(~ STATE_ALPHA)

bw_wa_m
ggsave(filename = "./viz/box_whisker/wa_bw_M.png", plot = bw_wa_m, dpi = 300, width = 7, height = 5)


#############################################################################################################################
# Individual states! Water use (cubic meters)
#############################################################################################################################

bw_func <- function(state,top) {
  
        df <- bw_data_fix %>% 
          filter(STATE_ALPHA == state)
        
bw <- ggplot(df, aes(x = YEAR, y = WATER_M, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Water Applied (cubic meters/hectare)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,top)) 

bw
ggsave(paste0(filename = './viz/box_whisker/wa_bw_',state,'_M.png'), plot = bw, dpi = 300, width = 7, height = 5)

}

bw_func("ARIZONA",35000)
bw_func("CALIFORNIA",25000)
bw_func("COLORADO",15000)
bw_func("IDAHO",16000)
bw_func("NEVADA",20000)
bw_func("NEW MEXICO",16000)
bw_func("MONTANA",13000)
bw_func("OREGON",16000)
bw_func("UTAH",16000)
bw_func("WASHINGTON",16000)
bw_func("WYOMING",16000)

#############################################################################################################################
# facetted by crop! Water use (cubic meters)
#############################################################################################################################

bw_func <- function(crop,top) {
  
        df <- bw_data_fix %>% 
          filter(CROP == crop)
        
bw <- ggplot(df, aes(x = YEAR, y = WATER_M, fill = crop)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Water Applied (cubic meters/hectare)") +
  scale_fill_manual(name = "CROP", values=colors) +
  facet_wrap(~ STATE_ALPHA) +
  coord_cartesian(ylim = c(0,top)) 

bw
ggsave(paste0(filename = './viz/box_whisker/wa_bw_',crop,'_M.png'), plot = bw, dpi = 300, width = 7, height = 5)

}

bw_func("Alfalfa",35000)
bw_func("Corn Grain",17000)
bw_func("Corn Silage",25000)
bw_func("Other Hay",28000)
bw_func("Wheat",20000)

#############################################################################################################################
# let's pull some statistics out of these boxplots!
#############################################################################################################################

sos2 <- as.data.frame(ggplot_build(bw_wa_m)$data[[1]])
sos2 <- sos2[c(1,3:5,12:13)]

sos2$CROP <- ifelse(sos2$fill == "#74C476", "Alfalfa",
                   ifelse(sos2$fill == "#6BAED6", "Corn Grain",
                        ifelse(sos2$fill == "#9E9AC8", "Corn Silage",
                            ifelse(sos2$fill == "#FB6A4A", "Other Hay",
                                ifelse(sos2$fill == "#FD8D3C", "Wheat", 
                   "NA")))))

sos2$YEAR <- ifelse(sos2$group %in% c(1,2,3,4,5), 2003,
                   ifelse(sos2$group %in% c(6,7,8,9,10), 2008,
                          ifelse(sos2$group %in% c(11,12,13,14,15), 2013,
                                 ifelse(sos2$group %in% c(16,17,18,19,20), 2018,
                                        "NA"))))

#by crop-state
sd_sos2 <- sos2 %>% group_by(CROP, PANEL) %>% summarise(median = median(middle), sd = sd(middle))

sos2_merge <- merge(sos2,sd_sos2, by = c("CROP", "PANEL"), all = T)

sos2_merge <- sos2_merge %>% mutate(variable = "WA_M") %>% select(-fill)

# by crop
median_sos2 <- sos2 %>% group_by(CROP) %>% summarise(median = median(middle),
                                                     sd = sd(middle)) %>% mutate(variable = "WA_M")

bw_wa_all <- ggplot(bw_data_fix, aes(x = YEAR, y = WATER_M, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Water Applied (cubic meters/hectare)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,20000)) 

bw_wa_all
ggsave(filename = "./viz//box_Whisker/wa_bw_all_M.png", plot = bw_wa_all, dpi = 300, width = 7, height = 5)

# by crop-year
cy2 <- as.data.frame(ggplot_build(bw_wa_all)$data[[1]])
cy2 <- cy2[c(1,3:5,12:13)]

cy2$CROP <- ifelse(cy2$fill == "#74C476", "Alfalfa",
                   ifelse(cy2$fill == "#6BAED6", "Corn Grain",
                        ifelse(cy2$fill == "#9E9AC8", "Corn Silage",
                            ifelse(cy2$fill == "#FB6A4A", "Other Hay",
                                ifelse(cy2$fill == "#FD8D3C", "Wheat", 
                   "NA")))))
cy2$YEAR <- ifelse(cy2$group %in% c(1,2,3,4,5), 2003,
                   ifelse(cy2$group %in% c(6,7,8,9,10), 2008,
                          ifelse(cy2$group %in% c(11,12,13,14,15), 2013,
                                 ifelse(cy2$group %in% c(16,17,18,19,20), 2018,
                                        "NA"))))

x_cy2 <- cy2 %>% group_by(CROP) %>% summarise(median = median(middle), sd = sd(middle)) %>% mutate(variable = "WA_M")

cy2_merge <- merge(cy2,x_cy2, by = c("CROP"), all = T)

```


### Box and Whisker plot of yield by crop through time

*Key Variables Code Chunk AA)*: K51 (CORN GRAIN); K61 (CORN SILAGE); K81 (WHEAT); K141 (ALFALFA); K151, (OTHER HAY). 

Code chunk AA):

1) Build box and whisker plots of yield by crop and through time using the raw data (not transformed because not summarised to a county or state level). I have also suppressed all outliers and data < 0.95 percentile (removes maximum state-crop-year combination in all cases).

```{r AA) BOX AND WHISKER YIELD, warning = FALSE, message = FALSE, eval = FALSE}
# use data from code chunk above ensuring n >= 10
#############################################################################################################################
# facetted by state! Yield (pounds)
#############################################################################################################################

bw_y <- ggplot(bw_data_fix, aes(x = YEAR, y = YIELD_LBS, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) +
  theme(legend.position = "bottom") +
  ylab("Yield (pounds/acre)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,35000)) + 
  facet_wrap(~ STATE_ALPHA)

bw_y
ggsave(filename = "./viz/box_whisker/y_bw_LBS.png", plot = bw_y, dpi = 300, width = 7, height = 5)

#############################################################################################################################
# Individual states! Yield (pounds)
#############################################################################################################################

bw_func <- function(state,top) {
  
        df <- bw_data_fix %>% 
          filter(STATE_ALPHA == state)
        
bw <- ggplot(df, aes(x = YEAR, y = YIELD_LBS, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Yield (pounds/acre)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,top)) 

bw
ggsave(paste0(filename = './viz/box_whisker/y_bw_',state,'_LBS.png'), plot = bw, dpi = 300, width = 7, height = 5)

}

bw_func("ARIZONA",30000)
bw_func("CALIFORNIA",35000)
bw_func("COLORADO",30000)
bw_func("IDAHO",30000)
bw_func("NEVADA",20000)
bw_func("NEW MEXICO",30000)
bw_func("MONTANA",25000)
bw_func("OREGON",30000)
bw_func("UTAH",27500)
bw_func("WASHINGTON",30000)
bw_func("WYOMING",22000)

#############################################################################################################################
# Facetted by states! Yield (kilograms)
#############################################################################################################################

bw_y_kg <- ggplot(bw_data_fix, aes(x = YEAR, y = YIELD_KG, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) +
  theme(legend.position = "bottom") +
  ylab("Yield (kilograms/hectare)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,35000)) + 
  facet_wrap(~ STATE_ALPHA)

bw_y_kg
# save plot
ggsave(filename = "./viz/box_whisker/y_bw_KG.png", plot = bw_y_kg, dpi = 300, width = 7, height = 5)

#############################################################################################################################
# Individual states! Yield (kilograms)
#############################################################################################################################

bw_func <- function(state,top) {
  
        df <- bw_data_fix %>% 
          filter(STATE_ALPHA == state)
        
bw <- ggplot(df, aes(x = YEAR, y = YIELD_KG, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Yield (kilograms/hectare)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,top)) 

bw
ggsave(paste0(filename = './viz/box_whisker/y_bw_',state,'_KG.png'), plot = bw, dpi = 300, width = 7, height = 5)

}

bw_func("ARIZONA",32000)
bw_func("CALIFORNIA",35000)
bw_func("COLORADO",30000)
bw_func("IDAHO",30000)
bw_func("NEVADA",20000)
bw_func("NEW MEXICO",30000)
bw_func("MONTANA",27000)
bw_func("OREGON",32000)
bw_func("UTAH",30000)
bw_func("WASHINGTON",32000)
bw_func("WYOMING",24000)

#############################################################################################################################
# facetted by crop! Yield (kilograms)
#############################################################################################################################

bw_func <- function(crop,top) {
  
        df <- bw_data_fix %>% 
          filter(CROP == crop)
        
bw <- ggplot(df, aes(x = YEAR, y = YIELD_KG, fill = crop)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Yield (kilograms/hectare)") +
  scale_fill_manual(name = "CROP", values=colors) +
  facet_wrap(~ STATE_ALPHA) +
  coord_cartesian(ylim = c(0,top)) 

bw
ggsave(paste0(filename = './viz/box_whisker/y_bw_',crop,'_KG.png'), plot = bw, dpi = 300, width = 7, height = 5)

}

bw_func("Alfalfa",30000)
bw_func("Corn Grain",20000)
bw_func("Corn Silage",34000)
bw_func("Other Hay",37000)
bw_func("Wheat",12000)

#############################################################################################################################
# let's pull some statistics out of these boxplots!
#############################################################################################################################

sos3 <- as.data.frame(ggplot_build(bw_y_kg)$data[[1]])
sos3 <- sos3[c(1,3:5,12:13)]

sos3$CROP <- ifelse(sos3$fill == "#74C476", "Alfalfa",
                   ifelse(sos3$fill == "#6BAED6", "Corn Grain",
                        ifelse(sos3$fill == "#9E9AC8", "Corn Silage",
                            ifelse(sos3$fill == "#FB6A4A", "Other Hay",
                                ifelse(sos3$fill == "#FD8D3C", "Wheat", 
                   "NA")))))
sos3$YEAR <- ifelse(sos3$group %in% c(1,2,3,4,5), 2003,
                   ifelse(sos3$group %in% c(6,7,8,9,10), 2008,
                          ifelse(sos3$group %in% c(11,12,13,14,15), 2013,
                                 ifelse(sos3$group %in% c(16,17,18,19,20), 2018,
                                        "NA"))))

# by crop=state
sd_sos3 <- sos3 %>% group_by(CROP, PANEL) %>% summarise(median = median(middle), sd = sd(middle))

sos3_merge <- merge(sos3,sd_sos3, by = c("CROP", "PANEL"), all = T)

sos3_merge <- sos3_merge %>% mutate(variable = "Y_KG") %>% select(-fill)

# by crop
median_sos3 <- sos3 %>% group_by(CROP) %>% summarise(median = median(middle),
                                            sd = sd(middle)) %>% mutate(variable = "Y_KG") 


bw_y_all <- ggplot(bw_data_fix, aes(x = YEAR, y = YIELD_KG, fill = CROP)) +
  geom_boxplot(outlier.shape = NA) +
  labs(fill = "Crop") +
  theme_classic() +
  theme(axis.title=element_text(size=14),
        axis.text=element_text(size=11)) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_text(size=9, angle=45, hjust = 1)) +
  theme(strip.background = element_blank()) + 
  theme(legend.position = "bottom") +
  ylab("Yield (kilograms/hectare)") +
  scale_fill_manual(name = "CROP", values=colors) +
  coord_cartesian(ylim = c(0,35000)) 

bw_y_all
ggsave(filename = "./viz//box_Whisker/y_bw_all_KG.png", plot = bw_y_all, dpi = 300, width = 7, height = 5)
# by crop-year
cy3 <- as.data.frame(ggplot_build(bw_wa_all)$data[[1]])
cy3 <- cy3[c(1,3:5,12:13)]

cy3$CROP <- ifelse(cy3$fill == "#74C476", "Alfalfa",
                   ifelse(cy3$fill == "#6BAED6", "Corn Grain",
                        ifelse(cy3$fill == "#9E9AC8", "Corn Silage",
                            ifelse(cy3$fill == "#FB6A4A", "Other Hay",
                                ifelse(cy3$fill == "#FD8D3C", "Wheat", 
                   "NA")))))
cy3$YEAR <- ifelse(cy3$group %in% c(1,2,3,4,5), 2003,
                   ifelse(cy3$group %in% c(6,7,8,9,10), 2008,
                          ifelse(cy3$group %in% c(11,12,13,14,15), 2013,
                                 ifelse(cy3$group %in% c(16,17,18,19,20), 2018,
                                        "NA"))))

x_cy3 <- cy3 %>% group_by(CROP) %>% summarise(median = median(middle), sd = sd(middle)) %>% mutate(variable = "Y_KG")

cy3_merge <- merge(cy3,x_cy3, by = c("CROP"), all = T)

# cbind all datasets out to .RDSs
bw_st_stats <- rbind(sos_merge, sos2_merge, sos3_merge) # by CROP and STATE

bw_panel_stats <- rbind(median_sos, median_sos2, median_sos3) # by CROP

bw_yr_stats <- rbind(cy_merge, cy2_merge, cy3_merge) # by CROP and YEAR
#####################################################################################################################
### build maximum and unique STPOID into each dataframe above
#####################################################################################################################
#### st
st <- bw_data_fix %>% 
  group_by(CROP, YEAR, STATE_ALPHA) %>% 
  summarise(unique = length(unique((STPOID))),
              max_wu = max(WATER_M),
              max_yield = max(YIELD_KG))

# merge
STATE_ALPHA <- c("ARIZONA", "CALIFORNIA", "COLORADO", "IDAHO", "MONTANA", "NEVADA", "NEW MEXICO", "OREGON", "UTAH", "WASHINGTON", "WYOMING")
PANEL <- c(1:11)

foo <- as.data.frame(cbind(STATE_ALPHA, PANEL))

foo <- merge(bw_st_stats, foo, by = "PANEL")

bw_st_stats <- merge(st, foo, by = c("CROP", "YEAR", "STATE_ALPHA"))
            
#### panel
panel <- bw_data_fix %>% 
  group_by(CROP) %>% 
  summarise(unique = length(unique((STPOID))),
              max_wu = max(WATER_M),
              max_yield = max(YIELD_KG))
# merge
bw_panel_stats <- merge(bw_panel_stats, panel, by = "CROP")

#### year
yr <- bw_data_fix %>% 
  group_by(CROP, YEAR) %>% 
  summarise(unique = length(unique((STPOID))),
              max_wu = max(WATER_M),
              max_yield = max(YIELD_KG))
# merge
bw_yr_stats <- merge(bw_yr_stats, panel, by = "CROP")

# save
write.csv(bw_st_stats, "./out/DEER/bw_st_stats.csv")
write.csv(bw_panel_stats, "./out/DEER/bw_panel_stats.csv")
write.csv(bw_yr_stats, "./out/DEER/bw_yr_stats.csv")

```