NASIS_data_to_EDIT_multiple.Rmd

---
output: 
  html_document:
    code_folding: hide
params:
  ecosites_of_interest: 'placeholder'
  mapunit_boundaries: 'placeholder'
  ecoregion_boundaries: 'placeholder'
  mlra_boundary: 'placeholder'
---
Author: Nathan Roe

Contact: nateroe@gmail.com

Date report was produced: `r format(Sys.time(), '%d %B, %Y')`

Data generated using NASIS Reports > MLRA02_Davis > EXPORT - Ecological site concept data by MUKEY list v4

Visit the [GitHub repository](https://github.com/natearoe/EDIT_ecosite_data) for version control information and instructions for producing these reports yourself. 


<!-- This is the final step in the 'EDIT ecosite data' workflow. It will produce reports for all the ecosites that you provide -->
<!-- to the corresponding for_loop_file.R.  -->

<!-- Note: You will need to change the file path on line 28. The easiest way to do that is to find the file that you created in -->
<!-- Step 2, your NASIS ecosite report. It should be in the folder with all your other 'EDIT ecosite data' files. Shift + right click > copy as path. Then paste the path in, replacing the existing one. That path will have backslashes \ these all need to be changed to forwardslashes /.  -->

<!-- If you want this file to produce maps of ecosites, you will need to remove the "#" before each line in the last two code blocks and follow the   -->

<!-- #### Packages -->
```{r class.source = "fold-show", message=FALSE, warning=FALSE, echo=FALSE }

```

```{r include = FALSE, echo=FALSE}
knitr::opts_chunk$set(fig.height = 7, fig.width = 9)
```

<!-- **Read in NASIS report** -->
```{r class.source = "fold-show", echo=FALSE}
EDIT_updates <- read.csv("./ecosite_report.csv")
```

<!-- **Remove miscellaneous areas** -->
```{r  class.source = "fold-show", echo=FALSE}
EDIT_updates <- EDIT_updates |>  dplyr::filter(compkind != "miscellaneous area")
```

<!-- **Define ecosite of interest** -->
```{r  class.source = "fold-show", echo=FALSE}
EDIT_updates_filtered <- EDIT_updates |>  dplyr::filter(ecosite_id == ecosites_of_interest)
```

<!-- **Determine component acreage unit** -->
```{r  class.source = "fold-show", echo=FALSE}

comp_acreage_unit <- ifelse(min(EDIT_updates_filtered$comp_acres[EDIT_updates_filtered$comp_acres != 0], na.rm = TRUE) < 100, 1, 
               ifelse(min(EDIT_updates_filtered$comp_acres[EDIT_updates_filtered$comp_acres != 0], na.rm = TRUE) < 1000, 10, 
                      ifelse(min(EDIT_updates_filtered$comp_acres[EDIT_updates_filtered$comp_acres != 0], na.rm = TRUE) < 10000, 100, 
                             ifelse(min(EDIT_updates_filtered$comp_acres[EDIT_updates_filtered$comp_acres != 0], na.rm = TRUE) < 
                                      100000, 1000, 10000))))

```

---
title: <center> `r EDIT_updates_filtered$ecosite_id[1]` ecosite report </center>
---


#### **Important Discussion:**

**Purpose of this report:**

This report summarizes data from components correlated to an ecosite. It could be used to populate EDIT, to QC component-ecosite correlations, or as documentation of ecosite properties. Of particular importance to QC, the user can hover their cursor over figures

this document allows you to hover your cursor over figures to see what component ids (coiids) are being represented. This allows you to hover over values that look like outliers and identify which coiid(s) might be improperly correlated.

**# of components vs. component acreage:**

The following analysis provides two ways of looking at data - (1) by component and (2) by component acreage. Looking at components (1) means that we treat all components equally, regardless of whether one component has a total acreage of 5 acres and another component has a total acreage of 50,000 acres. Looking at data by component acreage (2) acknowledges that components may differ in their importance based on their acreage. Therefore, you will see figures below that summarize data both by component (1) and by component acreage (2). 

I would suggest that looking at components (1) is the best way to QC your data. The correlation of components to ecosites should be assessed based on the *similarity of properties* between the ecosite and component. A bad component-ecosite correlation can just as easily happen for a component with low acreage as it can for a component with high acreage. For the QC process, you should think of component acreage as a metric of how important it is to have a good component-ecosite correlation. A bad component-ecosite correlation for a high acreage component will be more disruptive to the ecosite concept than a bad component-ecosite correlation for a low acreage component. Why? Continue reading...

**Standardized calculations:**

This report provides a standardized way of populating the range of properties for an ecosite. This is accomplished by weighting components by their acreage, hence why good correlations for high acreage components are important. Component properties are replicated by their associated acres. After replication, the following calculations are used:

Low (actual) - 5th percentile of the component low values
High (actual) - 95th percentile of the component high values
Low (representative) - 20th percentile of the combined component high and component low values
High (representative) - 80th percentile of the combined component high and component low values

As an example of the replication process, let's say we have two components - Component A and Component B. Component A has an acreage of 1. Component B has an acreage of 5. The elevation of Component A is 100 ft. The elevation of Component B is 500 ft. The resulting data set would be: 100, 500, 500, 500, 500, 500. Those values would be used to calculate percentiles. 

NA values mean that there are missing values. Figures do not show NA values but summary statistics do, so be sure to observe the summary statistics below figures to understand how many missing values there are. For water table and restriction depth, NAs can be interpreted as an ordinal value that, for most purposes, are greater than quantitative values. NAs for ponding depth can be interpretted as an ordinal value that, for most purposes, are less than quantitative values. Therefore, NAs are included in statistical summaries for these variables, including calculations of percentiles. 

Values should be rounded appropriate to the data before being entered into EDIT for all numeric fields except pH. 

**EDIT does not accept non-numerics:** 

EDIT does not accept non-numeric values being entered in numeric fields. If you enter non-numeric values, they will be changed to zeros after you save and re-open the site. If a field is not-applicable you should uncheck the box rather than entering NAs. For example, if there is no observed restriction, you should uncheck the depth to restriction box. If you enter NAs those values will be converted to zeros, suggesting that there is a restriction at the soil surface. 

**Components with no acreage value populated or zero acreage**

This report relies on component acreage for calculating standardized values and visualizing the distribution of data. If possible, correct acreages should be entered in NASIS from components with a missing component acreage or zero component acreage. Components with missing acreage or zero acreage will only be represented in the figures showing the distribution of data by the number of components. Components with missing acreage or zero acreage will not be represented in the acreage weighted figures, summary statistics, or standardized values. 

```{r}
NA_acreage <- paste(EDIT_updates_filtered |> dplyr::filter(is.na(comp_acres)) |> dplyr::pull(coiid))
NA_acreage <- ifelse(length(NA_acreage) == 0, "No components missing acreage.")

acreage_zero <- paste(EDIT_updates_filtered |> dplyr::filter(comp_acres == 0) |> dplyr::pull(coiid))
acreage_zero <- ifelse(length(acreage_zero) == 0, "No components with zero acreage.")

data.frame(`Comp missing acreage` = NA_acreage,
           `Comp acreage zero` = acreage_zero, check.names = FALSE) |> knitr::kable()

```

```{r, }
EDIT_updates_filtered$comp_acres <- ifelse(is.na(EDIT_updates_filtered$comp_acres), 0, EDIT_updates_filtered$comp_acres)
```

<center><h3> **Physiographic features tab:** </h3></center>

#### **Landform** 

Several properties in EDIT have the option to include multiple values, including landform. Other properties that allow for multiple values to be populated include parent material and surface texture. If there are multiple distinctly different properties that are associated with considerable acreage, list them. For most situations, one to two values is ideal. More than three values is most likely excessive. 

**Full landform string table by acreage and number of components**
```{r}
if (all(is.na(EDIT_updates_filtered$landform)) | 
    all(EDIT_updates_filtered$landform == "")) {print("No data populated in NASIS for this field.")
  
} else{
  
    landform.df <- EDIT_updates_filtered |> dplyr::select(landform, comp_acres, coiid) |> dplyr::group_by(landform) |> 
    dplyr::summarise(acres = sum(comp_acres), `# components` = dplyr::n(), coiids = toString(coiid)) |> 
    dplyr::arrange(dplyr::desc(acres)) 
  
    landform.df |> knitr::kable()

}
```

Wordcloud of acreage weighted landform descriptions. The wordcloud should give an overview of the words used to describe the landform, please view the full landform strings above for explicit details. 
```{r, message=FALSE, warning=FALSE}
if (all(is.na(EDIT_updates_filtered$landform))) {print("No data populated in NASIS for this field.")
  
} else{
  
  acreage_unit <- ifelse(min(landform.df$acres[landform.df$acres != 0], na.rm = TRUE) < 100, 1, 
               ifelse(min(landform.df$acres[landform.df$acres != 0], na.rm = TRUE) < 1000, 10, 
                      ifelse(min(landform.df$acres[landform.df$acres != 0], na.rm = TRUE) < 10000, 100, 
                             ifelse(min(landform.df$acres[landform.df$acres != 0], na.rm = TRUE) < 
                                      100000, 1000, 10000))))
  
  
 landform_words <-  rep(landform.df$landform, plyr::round_any(landform.df$acres/acreage_unit, accuracy = 1, f = round)) |> strsplit(" ") |> unlist() |>   
              stringr::str_replace(",", "")
 
 wordcloud::wordcloud(landform_words)

}
```

**Landform wordcloud table**
```{r, message=FALSE, warning=FALSE}
if (all(is.na(EDIT_updates_filtered$landform))) {print("No data populated in NASIS for this field.")
  
} else{
  
  table(landform_words) |>  as.data.frame() |>  dplyr::arrange(desc(Freq)) |>  dplyr::rename(`Landform words` =
                                                                    landform_words) |> dplyr::top_n(15)
  
}

```

#### **Runoff class**

This field is obsolete and will be removed from EDIT. Enter medium for all values.


#### **Flooding freq.**

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

EDIT_updates_filtered$flood_freq <- ifelse(EDIT_updates_filtered$flood_freq == "", "none", EDIT_updates_filtered$flood_freq)

if (all(is.na(EDIT_updates_filtered$flood_freq))) {print("No data populated in NASIS for this field.")
  
} else{
  
    my_levels <- c(
    "none",
    "very rare",
    "rare",
    "occasional",
    "frequent",
    "very frequent"
  )
  
  flood.freq.df <- EDIT_updates_filtered |>  dplyr::select(coiid, flood_freq, comp_acres)  
  
  flood.freq.df$flood_freq <- factor(flood.freq.df$flood_freq, levels = my_levels)
    
# Component histogram 
flood.freq1 <- ggplot2::ggplot(flood.freq.df, ggplot2::aes(flood_freq)) + 
  ggplot2::geom_bar(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Flood frequency") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Flood frequency")

# Convert to plotly
flood.freq1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage histogram
flood.freq2 <- ggplot2::ggplot(flood.freq.df, ggplot2::aes(flood_freq, comp_acres)) +
  ggplot2::geom_col(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Flood frequency") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Flood frequency")

#Convert to plotly
flood.freq2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures
manipulateWidget::combineWidgets(flood.freq1,
                                 flood.freq2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))
  
}

```

**Flooding frequency acreage weighted summary statistics**
```{r}
if (all(is.na(EDIT_updates_filtered$flood_freq))) {print("No data populated in NASIS for this field.")
  
} else{
  
  EDIT_updates_filtered |>  dplyr::select(coiid, flood_freq, comp_acres) |>  dplyr::group_by(flood_freq) |>  
    dplyr::summarise(acres = sum(comp_acres),  '# components' = dplyr::n()) |>  dplyr::arrange(desc(acres)) |>  as.data.frame()

}
```

**Flooding frequency acreage weighted percentiles.**

Actual low (5%) and actual high (95%). Representative low (20%) and representative high (80%).
```{r}
if (all(is.na(EDIT_updates_filtered$flood_freq))) {print("No data populated in NASIS for this field.")
  
} else{
  
      expanded_flood_freq <- rep(EDIT_updates_filtered$flood_freq, 
                               plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit,
                                               accuracy = 1, f = round))
    
    expanded_flood_freq <- expanded_flood_freq[order(match(expanded_flood_freq, my_levels))]
  
    data.frame("5%" = expanded_flood_freq[ceiling(0.05 * length(expanded_flood_freq))],
               "20%" = expanded_flood_freq[ceiling(0.20 * length(expanded_flood_freq))],
               "80%" = expanded_flood_freq[ceiling(0.80 * length(expanded_flood_freq))],
             "95%" = expanded_flood_freq[ceiling(0.95 * length(expanded_flood_freq))],
             check.names = FALSE)

}

```

#### **Flooding duration**

Put "--" for duration if frequency is "none"

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

EDIT_updates_filtered$flood_dur <- ifelse(EDIT_updates_filtered$flood_dur == "", "none", EDIT_updates_filtered$flood_dur)


if (all(is.na(EDIT_updates_filtered$flood_dur))) {print("No data populated in NASIS for this field.")
  
} else{
  
    my_levels <- c(
    "none",
    "extremely brief",
    "very brief",
    "brief",
    "long",
    "very long"
  )
  
  flood.dur.df <- EDIT_updates_filtered |>  dplyr::select(coiid, flood_dur, comp_acres)  
  
  flood.dur.df$flood_dur <- factor(flood.dur.df$flood_dur, levels = my_levels)
    
# Component histogram 
flood.dur1 <- ggplot2::ggplot(flood.dur.df, ggplot2::aes(flood_dur)) + 
  ggplot2::geom_bar(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Flood duration") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Flood duration")

# Convert to plotly
flood.dur1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage histogram
flood.dur2 <- ggplot2::ggplot(flood.dur.df, ggplot2::aes(flood_dur, comp_acres)) +
  ggplot2::geom_col(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Flood duration") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Flood duration")

#Convert to plotly
flood.dur2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures
manipulateWidget::combineWidgets(flood.dur1,
                                 flood.dur2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))
  
}

```

**Flooding duration acreage weighted summary statistics**
```{r}
if (all(is.na(EDIT_updates_filtered$flood_dur))) {print("No data populated in NASIS for this field.")
  
} else{
  
  EDIT_updates_filtered |>  dplyr::select(coiid, flood_dur, comp_acres) |>  dplyr::group_by(flood_dur) |>  
    dplyr::summarise(acres = sum(comp_acres),  '# components' = dplyr::n()) |>  dplyr::arrange(desc(acres)) |>  as.data.frame()

}
```

**Flooding duration acreage weighted statistics.**

Actual low (5%) and actual high (95%). Representative low (20%) and representative high (80%).
```{r}
if (all(is.na(EDIT_updates_filtered$flood_dur))) {print("No data populated in NASIS for this field.")
  
} else{
  
    expanded_flood_dur <- rep(EDIT_updates_filtered$flood_dur, 
                               plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit,
                                               accuracy = 1, f = round))
    
    expanded_flood_dur <- expanded_flood_dur[order(match(expanded_flood_dur, my_levels))]
  
    data.frame("5%" = expanded_flood_dur[ceiling(0.05 * length(expanded_flood_dur))],
               "20%" = expanded_flood_dur[ceiling(0.20 * length(expanded_flood_dur))],
               "80%" = expanded_flood_dur[ceiling(0.80 * length(expanded_flood_dur))],
             "95%" = expanded_flood_dur[ceiling(0.95 * length(expanded_flood_dur))],
             check.names = FALSE)

}

```

#### **Ponding freq.**

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

EDIT_updates_filtered$pond_freq <- ifelse(EDIT_updates_filtered$pond_freq == "", "none", EDIT_updates_filtered$pond_freq)


if (all(is.na(EDIT_updates_filtered$pond_freq))) {print("No data populated in NASIS for this field.")
  
} else{
  
  my_levels <- c(
    "none",
    "rare",
    "occasional",
    "frequent",
    "common"
  )
  
  pond.freq.df <- EDIT_updates_filtered |>  dplyr::select(coiid, pond_freq, comp_acres)  
  
  pond.freq.df$pond_freq <- factor(pond.freq.df$pond_freq, levels = my_levels)
    
# Component histogram 
pond.freq1 <- ggplot2::ggplot(pond.freq.df, ggplot2::aes(pond_freq)) + 
  ggplot2::geom_bar(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Ponding freq.") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Ponding freq.")

# Convert to plotly
pond.freq1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage histogram
pond.freq2 <- ggplot2::ggplot(pond.freq.df, ggplot2::aes(pond_freq, comp_acres)) +
  ggplot2::geom_col(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Ponding freq.") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Ponding freq.")

#Convert to plotly
pond.freq2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures
manipulateWidget::combineWidgets(pond.freq1,
                                 pond.freq2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))
  
}

```

**Ponding frequency acreage weighted summary statistics**
```{r}
if (all(is.na(EDIT_updates_filtered$pond_freq))) {print("No data populated in NASIS for this field.")
  
} else{
  
  EDIT_updates_filtered |>  dplyr::select(coiid, pond_freq, comp_acres) |>  dplyr::group_by(pond_freq) |>  
    dplyr::summarise(acres = sum(comp_acres),  '# components' = dplyr::n()) |>  dplyr::arrange(desc(acres)) |>  as.data.frame()

}
```

**Ponding frequency acreage weighted statistics.**

Actual low (5%) and actual high (95%). Representative low (20%) and representative high (80%).
```{r}
if (all(is.na(EDIT_updates_filtered$pond_freq))) {print("No data populated in NASIS for this field.")
  
} else{
  
    expanded_pond_freq <- rep(EDIT_updates_filtered$pond_freq, 
                               plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit,
                                               accuracy = 1, f = round))
    
    expanded_pond_freq <- expanded_pond_freq[order(match(expanded_pond_freq, my_levels))]
  
    data.frame("5%" = expanded_pond_freq[ceiling(0.05 * length(expanded_pond_freq))],
               "20%" = expanded_pond_freq[ceiling(0.20 * length(expanded_pond_freq))],
               "80%" = expanded_pond_freq[ceiling(0.80 * length(expanded_pond_freq))],
             "95%" = expanded_pond_freq[ceiling(0.95 * length(expanded_pond_freq))],
             check.names = FALSE)

}

```

#### **Ponding duration**

Put "--" for duration if frequency is "none"

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

EDIT_updates_filtered$pond_dur <- ifelse(EDIT_updates_filtered$pond_dur == "", "none", EDIT_updates_filtered$pond_dur)


if (all(is.na(EDIT_updates_filtered$pond_dur))) {print("No data populated in NASIS for this field.")
  
} else{
  
  my_levels <- c(
    "none",
    "very brief",
    "brief",
    "long",
    "very long"
  )
  
  pond.dur.df <- EDIT_updates_filtered |>  dplyr::select(coiid, pond_dur, comp_acres)  
  
  pond.dur.df$pond_dur <- factor(pond.dur.df$pond_dur, levels = my_levels)
    
# Component histogram 
pond.dur1 <- ggplot2::ggplot(pond.dur.df, ggplot2::aes(pond_dur)) + 
  ggplot2::geom_bar(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Ponding duration") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Ponding duration")

# Convert to plotly
pond.dur1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage histogram
pond.dur2 <- ggplot2::ggplot(pond.dur.df, ggplot2::aes(pond_dur, comp_acres)) +
  ggplot2::geom_col(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Ponding duration") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Ponding duration")

#Convert to plotly
pond.dur2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures
manipulateWidget::combineWidgets(pond.dur1,
                                 pond.dur2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))
  
}

```

**Ponding duration acreage weighted summary statistics**
```{r}
if (all(is.na(EDIT_updates_filtered$pond_dur))) {print("No data populated in NASIS for this field.")
  
} else{
  
  EDIT_updates_filtered |>  dplyr::select(coiid, pond_dur, comp_acres) |>  dplyr::group_by(pond_dur) |>  
    dplyr::summarise(acres = sum(comp_acres),  '# components' = dplyr::n()) |>  dplyr::arrange(desc(acres)) |>  as.data.frame()

}
```

**Ponding duration acreage weighted statistics.**

5th and 95th percentiles should be used for 'low (actual)' and 'high (actual)'. The 20th and 80th percentiles should be used for the 'low (representative)' and 'high (representative)'. 
```{r}
if (all(is.na(EDIT_updates_filtered$pond_dur))) {print("No data populated in NASIS for this field.")
  
} else{
  
    expanded_pond_dur <- rep(EDIT_updates_filtered$pond_dur, 
                               plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit,
                                               accuracy = 1, f = round))
    
    expanded_pond_dur <- expanded_pond_dur[order(match(expanded_pond_dur, my_levels))]
  
    data.frame("5%" = expanded_pond_dur[ceiling(0.05 * length(expanded_pond_dur))],
               "20%" = expanded_pond_dur[ceiling(0.20 * length(expanded_pond_dur))],
               "80%" = expanded_pond_dur[ceiling(0.80 * length(expanded_pond_dur))],
             "95%" = expanded_pond_dur[ceiling(0.95 * length(expanded_pond_dur))],
             check.names = FALSE)

}

```

#### **Elevation**

Elevation prep code
```{r, message=FALSE}
if (all(is.na(EDIT_updates_filtered$elev_h_ft)) &
    all(is.na(EDIT_updates_filtered$elev_l_ft))) {print("No data populated in NASIS for this field.")
  
} else{

# Define bin width and round by measures
  
my_vector <- c(50, 100, 250, 500, 1000, 2000, 5000)

my_factor <- diff(range(c(EDIT_updates_filtered$elev_h_ft, 
             EDIT_updates_filtered$elev_l_ft), na.rm = TRUE))/my_vector



reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe
  
elev.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, elev_h_ft, elev_l_ft, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "elev",
    values_to = "value",
    cols = c("elev_h_ft", "elev_l_ft")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$elev_l_ft, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$elev_h_ft, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$elev_l_ft, na.rm = TRUE) &
             min(EDIT_updates_filtered$elev_l_ft, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$elev_l_ft, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$elev_h_ft, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq  
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

}

```


```{r, warning=FALSE, message=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$elev_h_ft))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

el_h1 <- ggplot2::ggplot(elev.df |>  dplyr::filter(elev == "elev_h_ft"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2, color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Elevation high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Elevation (ft.)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

el_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage histogram

el_h2 <- ggplot2::ggplot(elev.df |>  dplyr::filter(elev == "elev_h_ft"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Elevation high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Elevation (ft.)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)]
  )

# Convert acreage histogram to plotly

el_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(el_h1,
                                 el_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))


}

```

**Elevation high acreage weighted 95th percentile**  -  use as 'high (actual)'

Round elevation in a manner consistent with your data. For many situations, it will likely be to the hundreds place. 
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$elev_h_ft))) {print("No data populated in NASIS for this field.")
  
} else{
  
  elev_h_stats <- rep(EDIT_updates_filtered$elev_h_ft, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                accuracy = 1, f = round)) |> sort(na.last = NA)
  
  elev_h_stats |> quantile(probs = c(0.95), na.rm = TRUE, type = 1)

}

```

**Elevation high acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$elev_h_ft))) {print("No data populated in NASIS for this field.")
  
} else{

  summary(elev_h_stats)

}
```

```{r, warning=FALSE, message=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$elev_l_ft))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

el_l1 <- ggplot2::ggplot(elev.df |>  dplyr::filter(elev == "elev_l_ft"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Elevation low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Elevation (ft.)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

el_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

el_l2 <- ggplot2::ggplot(elev.df |>  dplyr::filter(elev == "elev_l_ft"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Elevation low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Elevation (ft.)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

el_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(el_l1,
                                 el_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```
 
**Elevation low acreage weighted 5th percentile** - use as 'low (actual)'
```{r}

if (all(is.na(EDIT_updates_filtered$elev_l_ft))) {print("No data populated in NASIS for this field.")
  
} else{
  
    elev_l_stats <- rep(EDIT_updates_filtered$elev_l_ft, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                          accuracy = 1, f = round)) |> sort(na.last = NA)
    
    elev_l_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```

**Elevation low acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$elev_l_ft))) {print("No data populated in NASIS for this field.")
  
} else{

  summary(elev_l_stats)

}

```


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

if (all(is.na(EDIT_updates_filtered$elev_l_ft & EDIT_updates_filtered$elev_h_ft))) {print("No data populated in NASIS for this field.")
  
} else{
  
  
el_r1 <- ggplot2::ggplot(elev.df |>  dplyr::filter(elev %in% c("elev_h_ft", "elev_l_ft")) |> 
  dplyr::mutate(highlow = ifelse(elev == "elev_h_ft", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")),
                         ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Elevation representative") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Elevation (ft.)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert histogram to plotly

el_r1 <- plotly::ggplotly(tooltip = c("text", "value"))


# Acreage weighted histogram

el_r2 <- ggplot2::ggplot(elev.df |>  dplyr::filter(elev %in% c("elev_h_ft", "elev_l_ft")) |> 
  dplyr::mutate(highlow = ifelse(elev == "elev_h_ft", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")),
                         ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Elevation representative") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Elevation (ft.)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

el_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(el_r1,
                                 el_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))


}

```

**Elevation representative - combined component high/low 20th and 80th percentile** - use as 'low (representative)' & 'high (representative)' 
```{r}
if (all(is.na(EDIT_updates_filtered$elev_l_ft & EDIT_updates_filtered$elev_h_ft))) {print("No data populated in NASIS for this field.")
  
} else{
  
  elev_r_stats <- rep(elev.df$value, plyr::round_any(elev.df$acres/comp_acreage_unit, 
                            accuracy = 1, f = round)) |> sort(na.last = NA)
  
  elev_r_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE, type = 1)

}

```

**Elevation representative acreage weighted summary statistics**
```{r}
if (all(is.na(EDIT_updates_filtered$elev_l_ft & EDIT_updates_filtered$elev_h_ft))) {print("No data populated in NASIS for this field.")
  
} else{

  summary(elev_r_stats)

}

```

#### **Slope**

Slope prep code
```{r, message=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$slope_h)) & 
    all(is.na(EDIT_updates_filtered$slope_l))) {print("No data populated in NASIS for this field.")
  
} else{

# Define bin width and round by measures
  
my_vector <- c(1, 2, 5, 10, 15, 20)

my_factor <- diff(range(c(EDIT_updates_filtered$slope_h, 
             EDIT_updates_filtered$slope_l), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

slope.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, slope_h, slope_l, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "slope",
    values_to = "value",
    cols = c("slope_h", "slope_l")) |> 
  dplyr::rename(acres = comp_acres) 

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$slope_l, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$slope_h, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$slope_l, na.rm = TRUE) &
             min(EDIT_updates_filtered$slope_l, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$slope_l, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$slope_h, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq   
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$slope_h))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

slope_h1 <- ggplot2::ggplot(slope.df |>  dplyr::filter(slope == "slope_h"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Slope high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Slope (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

slope_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

slope_h2 <- ggplot2::ggplot(slope.df |>  dplyr::filter(slope == "slope_h"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Slope high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Slope (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

slope_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(slope_h1,
                                 slope_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Slope high acreage weighted 95th percentile** -  use as 'high (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$slope_h))) {print("No data populated in NASIS for this field.")
  
} else{

  slope_h_sum_stats <- rep(EDIT_updates_filtered$slope_h, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                               accuracy = 1, f = round)) |> sort(na.last = NA)
  
  slope_h_sum_stats |>  quantile(probs = c(0.95), na.rm = TRUE, type = 1)

}

```

**Slope high acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$slope_h))) {print("No data populated in NASIS for this field.")
  
} else{
  
summary(slope_h_sum_stats)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$slope_l))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

slope_l1 <- ggplot2::ggplot(slope.df |>  dplyr::filter(slope == "slope_l"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Slope low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Slope (ft.)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

slope_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

slope_l2 <- ggplot2::ggplot(slope.df |>  dplyr::filter(slope == "slope_l"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Slope low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Slope (ft.)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert to plotly

slope_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(slope_l1,
                                 slope_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Slope low acreage weighted 5th percentile** -  use as 'low (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$slope_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
  slope_l_sum_stats <- rep(EDIT_updates_filtered$slope_l, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                            accuracy = 1, f = round)) |> sort(na.last = NA)
  
  
  slope_l_sum_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```

**Slope low acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$slope_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
  summary(slope_l_sum_stats)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$slope_h)) & all(is.na(EDIT_updates_filtered$slope_l))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram
  
slope_r1 <- ggplot2::ggplot(slope.df |>  dplyr::filter(slope %in% c("slope_h", "slope_l")) |> 
  dplyr::mutate(highlow = ifelse(slope == "slope_h", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Slope representative") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Slope (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

slope_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

slope_r2 <- ggplot2::ggplot(slope.df |>  dplyr::filter(slope %in% c("slope_h", "slope_l")) |> 
  dplyr::mutate(highlow = ifelse(slope == "slope_h", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Slope representative") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Slope (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

slope_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(slope_r1,
                                 slope_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Slope representative - combined component high/low 20th and 80th percentile** - use as 'low (representative)' & 'high (representative)' 
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$slope_h)) & all(is.na(EDIT_updates_filtered$slope_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
  slope_r_sum_stats <- rep(slope.df$value, plyr::round_any(slope.df$acres/comp_acreage_unit, 
                                accuracy = 1, f = round)) |> sort(na.last = NA)
  
  slope_r_sum_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE)

}

```

**Slope representative acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$slope_h)) & all(is.na(EDIT_updates_filtered$slope_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
  summary(slope_r_sum_stats)

}
```
#### **Slope shape**

Acreage weighted slope shape. If there are multiple comma separated slope shapes listed for a component, half of the component acreage is attributed to each slope shape. To see the raw tabular data, see the next section. 
```{r, message=FALSE, warning=FALSE}
if (all(is.na(EDIT_updates_filtered$slope_shp_across_down))) {print("No data populated in NASIS for this field.")
  
} else{

### all possible slope shape combinations
possible_shapes <- expand.grid(c("convex", "linear", "concave"), c("convex", "linear", "concave")) |> 
  dplyr::rename(across = Var1, down = Var2)

possible_shapes$acres <- NA
possible_shapes$combined <- paste0(possible_shapes$across, "/", possible_shapes$down)

### now let's look at the slope shapes that exist in our data
slope_shape <- EDIT_updates_filtered |> dplyr::select(comp_acres, slope_shp_across_down)

numb_term <-  slope_shape$slope_shp_across_down |> stringr::str_count(",") + 1

slope_shape$comp_acres2 <- slope_shape$comp_acres/numb_term

slope_shape_new <- slope_shape |> dplyr::mutate(slope_shape_new = strsplit(slope_shape$slope_shp_across_down, ",")) |> tidyr::unnest(slope_shape_new)

slope_shape_new$slope_shape_new <- gsub(" ", "", slope_shape_new$slope_shape_new)

# now combine them and summarize

slope_shape_combined <-  slope_shape_new |> 
  dplyr::group_by(slope_shp_across_down) |> dplyr::summarise(sum = sum(comp_acres2))

# merge all possible with actual slope shapes
slope_shape_combined_all <- dplyr::full_join(possible_shapes, slope_shape_combined |> dplyr::rename(combined = slope_shp_across_down)) |> 
  dplyr::select(-acres) |> dplyr::rename(acres = sum)
slope_shape_combined_all$acres[is.na(slope_shape_combined_all$acres)] <- 0

# plot
ggplot2::ggplot(slope_shape_combined_all |> dplyr::select(-combined), ggplot2::aes(x = across, y = down)) +
  ggplot2::geom_raster(ggplot2::aes(fill = acres)) + 
  ggplot2::scale_fill_gradient(low = "grey90", high = "red") +
  ggplot2::labs(x = "across", y = "down", title = "Slope shape")

}
```

Raw slope shape with coiids
```{r, message=FALSE, warning=FALSE}
if (all(is.na(EDIT_updates_filtered$slope_shp_across_down))) {print("No data populated in NASIS for this field.")
  
} else{
    
    EDIT_updates_filtered |> dplyr::select(slope_shp_across_down, comp_acres, coiid) |> dplyr::group_by(slope_shp_across_down) |> 
    dplyr::summarise(acres = sum(comp_acres), `# components` = dplyr::n(), coiids = toString(coiid)) |> 
    dplyr::arrange(dplyr::desc(acres)) |> knitr::kable()
  

}
```

#### **Hillslope profile**

Acreage weighted hillslope profile. If there are multiple comma separated hillslope profiles listed for a component, half of the component acreage is attributed to each hillslope profile. To see the raw tabular data, see the next section. 
```{r, warning=FALSE, message=FALSE}
if (all(is.na(EDIT_updates_filtered$hillslope_profile))) {print("No data populated in NASIS for this field.")
  
} else{

hillslope <- EDIT_updates_filtered |> dplyr::select(hillslope_profile, coiid, comp_acres)

# some components have multiple slope shapes described. for those we will 
# divide the acres by the number of terms used to describe it 

numb_term <-  hillslope$hillslope_profile |> stringr::str_count(",") + 1

hillslope$comp_acres2 <- hillslope$comp_acres/numb_term

hillslope_new <- hillslope |> dplyr::mutate(hillslope_new = strsplit(hillslope$hillslope_profile, ",")) |> tidyr::unnest(hillslope_new)

hillslope_new$hillslope_new <- gsub(" ", "", hillslope_new$hillslope_new)

# now combine them and summarize
# hillslope_combined <- hillslope_new |> 
#   dplyr::group_by(hillslope_new) |> dplyr::summarise(acres = sum(comp_acres2), coiids = toString(coiid))

hillslope_combined <- hillslope_new |> 
  dplyr::group_by(hillslope_new) |> dplyr::summarise(acres = sum(comp_acres2), coiids = toString(coiid))


# plotting

my_levels <- c(
    "summit",
    "shoulder",
    "backslope",
    "footslope",
    "toeslope"
  )
  
  
hillslope_new$hillslope_new <- factor(hillslope_new$hillslope_new, levels = my_levels)
    
# Component histogram 
hillslope1 <- ggplot2::ggplot(hillslope_new, ggplot2::aes(hillslope_new)) + 
  ggplot2::geom_bar(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Hillslope profile") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Hillslope profile")

# Convert to plotly
hillslope1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage histogram
hillslope2 <- ggplot2::ggplot(hillslope_new, ggplot2::aes(hillslope_new, comp_acres2)) +
  ggplot2::geom_col(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Hillslope profile") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Hillslope profile")

#Convert to plotly
hillslope2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures
manipulateWidget::combineWidgets(hillslope1,
                                 hillslope2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}
  
```

Raw hillslope profile with coiids
```{r, warning=FALSE, message=FALSE}
if (all(is.na(EDIT_updates_filtered$hillslope_profile))) {print("No data populated in NASIS for this field.")
  
} else{
    
    EDIT_updates_filtered |> dplyr::select(hillslope_profile, comp_acres, coiid) |> dplyr::group_by(hillslope_profile) |> 
    dplyr::summarise(acres = sum(comp_acres), `# components` = dplyr::n(), coiids = toString(coiid)) |> 
    dplyr::arrange(dplyr::desc(acres)) |> knitr::kable()
  

}
```

#### **Ponding depth**

One thing to clarify about ponding depth - this is a measure of how deep the ponding is above the soil surface. If ecosite does not have a ponding depth, uncheck the box. 

Ponding depth prep code
```{r, message=FALSE}
if (all(is.na(EDIT_updates_filtered$pond_depth_h))  & all(is.na(EDIT_updates_filtered$pond_depth_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
# Define bin width and round by measures

my_vector <- c(1, 2, 5, 10, 15, 20)

my_factor <- diff(range(c(EDIT_updates_filtered$pond_depth_h, 
             EDIT_updates_filtered$pond_depth_l), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

pond_depth.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, pond_depth_l, pond_depth_h, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "pond_depth",
    values_to = "value",
    cols = c("pond_depth_l", "pond_depth_h")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$pond_depth_l, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$pond_depth_h, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$pond_depth_l, na.rm = TRUE) &
             min(EDIT_updates_filtered$pond_depth_l, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$pond_depth_l, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$pond_depth_h, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq   
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$pond_depth_h))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

pond_depth_h1 <- ggplot2::ggplot(pond_depth.df |>  dplyr::filter(pond_depth == "pond_depth_h"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Ponding depth high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Ponding depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  n.breaks = length(my_seq))

# Convert histogram to plotly

pond_depth_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

pond_depth_h2 <- ggplot2::ggplot(pond_depth.df |>  dplyr::filter(pond_depth == "pond_depth_h"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Ponding depth high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Ponding depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  n.breaks = length(my_seq))

# Convert to plotly

pond_depth_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(pond_depth_h1,
                                 pond_depth_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

NAs in ponding depth likely mean that no ponding was observed. NA's are not represented in the figures above. Please observe the table below to determine if NAs exist in your dataset.
```{r}
if (all(is.na(EDIT_updates_filtered$pond_depth_h))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# components NA' = sum(is.na(EDIT_updates_filtered$pond_depth_h)),
             '% components NA' = 100*sum(is.na(EDIT_updates_filtered$pond_depth_h))/length(EDIT_updates_filtered$pond_depth_h),
             check.names = FALSE)

  
}
```

**Ponding depth high acreage weighted 95th percentile** -  use as 'high (actual)'. This value may be NA.  

The percentile calculation for ponding depth is performed differently because NAs in the ponding depth likely mean no ponding was observed. For that reason, we will consider NAs to be an ordinal value that is lower than all other quantitative values. 
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$pond_depth_h))) {print("No data populated in NASIS for this field.")
  
} else{
  
    pond_depth_na_h_rem_expanded <- c(rep(EDIT_updates_filtered$pond_depth_h, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, accuracy = 1, f = round))) |> sort(na.last = FALSE)

   data.frame("95%" = pond_depth_na_h_rem_expanded[ceiling(0.95 * length(pond_depth_na_h_rem_expanded))], 
             check.names = FALSE)
  
}

```

**Ponding depth high acreage weighted summary statistics** - Summary statistics are calculated using numeric values (NA values omitted). The combined acreage of components with an NA value is listed. 
```{r}

if (all(is.na(EDIT_updates_filtered$pond_depth_h))) {print("No data populated in NASIS for this field. For ponding, this likely means no ponding was observed.")
  
} else{
  
  summary(pond_depth_na_h_rem_expanded, na.rm = TRUE)

}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$pond_depth_l))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

pond_depth_l1 <- ggplot2::ggplot(pond_depth.df |>  dplyr::filter(pond_depth == "pond_depth_l"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Ponding depth low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Ponding depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  n.breaks = length(my_seq))

# Convert histogram to plotly

pond_depth_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

pond_depth_l2 <- ggplot2::ggplot(pond_depth.df |>  dplyr::filter(pond_depth == "pond_depth_l"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Ponding depth low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Ponding depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  n.breaks = length(my_seq))

# Convert to plotly

pond_depth_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(pond_depth_l1,
                                 pond_depth_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

NAs in ponding depth likely mean that no ponding was observed. NA's are not represented in the figures above. Please observe the table below to determine if NAs exist in your dataset.
```{r}
if (all(is.na(EDIT_updates_filtered$pond_depth_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# components NA' = sum(is.na(EDIT_updates_filtered$pond_depth_l)),
             '% components NA' = 100*sum(is.na(EDIT_updates_filtered$pond_depth_l))/length(EDIT_updates_filtered$pond_depth_l),
             check.names = FALSE)
}

```

**Ponding depth low acreage weighted 5th percentile** -  use as 'low (actual)'. This value may be NA.  

The percentile calculation for ponding depth is performed differently because NAs in ponding depth likely mean no ponding was observed. For that reason, we will consider NAs to be an ordinal value that is less than all other quantitative values.
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$pond_depth_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
    pond_depth_na_l_rem_expanded <- c(rep(EDIT_updates_filtered$pond_depth_l, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, accuracy = 1, f = round))) |> sort(na.last = FALSE)

   data.frame("5%" = pond_depth_na_l_rem_expanded[ceiling(0.05 * length(pond_depth_na_l_rem_expanded))], 
             check.names = FALSE)
  
}

```

**Ponding depth low acreage weighted summary statistics** - Summary statistics are calculated using numeric values (NA values omitted). The combined acreage of components with an NA value is listed.
```{r}

if (all(is.na(EDIT_updates_filtered$pond_depth_l))) {print("No data populated in NASIS for this field. For ponding, this likely means no ponding was observed.")
  
} else{
  
summary(pond_depth_na_l_rem_expanded)

}
```



```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$pond_depth_l & EDIT_updates_filtered$pond_depth_h))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

pond_depth_r1 <- ggplot2::ggplot(pond_depth.df |>  dplyr::filter(pond_depth %in% c("pond_depth_h", "pond_depth_l")) |> 
  dplyr::mutate(highlow = ifelse(pond_depth == "pond_depth_h", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Ponding depth representative") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Ponding depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  n.breaks = length(my_seq))

# Convert histogram to plotly

pond_depth_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

pond_depth_r2 <- ggplot2::ggplot(pond_depth.df |>  dplyr::filter(pond_depth %in% c("pond_depth_h", "pond_depth_l")) |> 
  dplyr::mutate(highlow = ifelse(pond_depth == "pond_depth_h", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Ponding depth representative") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Ponding depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  n.breaks = length(my_seq))

# Convert to plotly

pond_depth_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(pond_depth_r1,
                                 pond_depth_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

NAs in ponding depth likely mean that no ponding was observed. NA's are not represented in the figures above. Please observe the table below to determine if NAs exist in your dataset.
```{r}
if (all(is.na(EDIT_updates_filtered$pond_depth_l & EDIT_updates_filtered$pond_depth_h))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# values NA' = sum(is.na(EDIT_updates_filtered$pond_depth_l)) + sum(is.na(EDIT_updates_filtered$pond_depth_h)),
             '% values NA' = 100*(sum(is.na(EDIT_updates_filtered$pond_depth_h)) + sum(is.na(EDIT_updates_filtered$pond_depth_l)))/
               (length(EDIT_updates_filtered$pond_depth_h) * 2), check.names = FALSE)
  
}

```

**Ponding depth representative - combined component high/low 20th and 80th percentile** - use as 'low (representative)' & 'high (representative)'. This value may be NA.  

The percentile calculation for ponding depth is performed differently because NAs in ponding depth likely mean no ponding was observed. For that reason, we will consider NAs to be an ordinal value that is less than all other quantitative values.
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$pond_depth_l & EDIT_updates_filtered$pond_depth_h))) {print("No data populated in NASIS for this field.")
  
} else{
  
  pond_depth_r_sum_stats <- rep(pond_depth.df$value, plyr::round_any(pond_depth.df$acres/comp_acreage_unit, accuracy = 1, f = round)) |> sort(na.last = FALSE)

  data.frame("20%" = pond_depth_r_sum_stats[ceiling(0.2 * length(pond_depth_r_sum_stats))],
             "80%" = pond_depth_r_sum_stats[ceiling(0.8 * length(pond_depth_r_sum_stats))],
             check.names = FALSE)
  
}

```

**Ponding depth representative acreage weighted summary statistics** - Summary statistics are calculated using numeric values (NA values omitted). The combined acreage of components with an NA value is listed.
```{r}

if (all(is.na(EDIT_updates_filtered$pond_depth_l & EDIT_updates_filtered$pond_depth_h))) {print("No data populated in NASIS for this field. For ponding, this likely means no ponding was observed.")
  
} else{
  
summary(pond_depth_r_sum_stats)

}
```

#### **Water table depth**

Water table depth prep code
```{r, message=FALSE}
# Is variable NA?



if (all(is.na(EDIT_updates_filtered$water_table_h_in)) & all(is.na(EDIT_updates_filtered$water_table_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
# Define bin width and round by measures

my_vector <- c(1, 2, 5, 10, 15, 20)

my_factor <- diff(range(c(EDIT_updates_filtered$water_table_h_in, 
             EDIT_updates_filtered$water_table_l_in), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

water_table.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, water_table_l_in, water_table_h_in, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "water_table",
    values_to = "value",
    cols = c("water_table_l_in", "water_table_h_in")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$water_table_l_in, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$water_table_h_in, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$water_table_l_in, na.rm = TRUE) &
             min(EDIT_updates_filtered$water_table_l_in, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$water_table_l_in, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$water_table_h_in, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq   
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$water_table_h_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

water_table_h1 <- ggplot2::ggplot(water_table.df |>  dplyr::filter(water_table == "water_table_h_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Water table high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Water table (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  n.breaks = length(my_seq))
  
# Convert histogram to plotly

water_table_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

water_table_h2 <- ggplot2::ggplot(water_table.df |>  dplyr::filter(water_table == "water_table_h_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Water table high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Water table (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  n.breaks = length(my_seq))

# Convert to plotly

water_table_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(water_table_h1,
                                 water_table_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

NAs in water table likely mean that no water table was observed. NA's are not represented in the figures above. Please observe the table below to understand the relative acreage classified as NA. 
```{r}
if (all(is.na(EDIT_updates_filtered$water_table_h_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# components NA' = sum(is.na(EDIT_updates_filtered$water_table_h_in)),
             '% components NA' = 100*sum(is.na(EDIT_updates_filtered$water_table_h_in))/length(EDIT_updates_filtered$water_table_h_in),
             check.names = FALSE) 
  
}
```

**Water table high acreage weighted 95th percentile** -  use as 'high (actual)'. This value may be NA.  

The percentile calculation for water table is performed differently because NAs in the water table likely mean no water table was observed. For that reason, we will consider NAs to be an ordinal value that is deeper than all other quantitative values. 
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$water_table_h_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
    water_table_na_h_rem_expanded <- c(rep(EDIT_updates_filtered$water_table_h_in, EDIT_updates_filtered$comp_acres/comp_acreage_unit)) |> sort(na.last = TRUE)

   data.frame("95%" = water_table_na_h_rem_expanded[ceiling(0.95 * length(water_table_na_h_rem_expanded))], 
             check.names = FALSE)

}

```

**Water table high acreage weighted summary statistics** - Summary statistics are calculated using numeric values (NA values omitted). The combined acreage of components with an NA value is listed.
```{r}

if (all(is.na(EDIT_updates_filtered$water_table_h_in))) {print("No data populated in NASIS for this field. For water table, this likely means no water table was observed.")
  
} else{
  
 summary(water_table_na_h_rem_expanded, na.rm = TRUE)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$water_table_l_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

water_table_l1 <- ggplot2::ggplot(water_table.df |>  dplyr::filter(water_table == "water_table_l_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Water table low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Water table (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  n.breaks = length(my_seq))

# Convert histogram to plotly

water_table_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

water_table_l2 <- ggplot2::ggplot(water_table.df |>  dplyr::filter(water_table == "water_table_l_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Water table low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Water table (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  n.breaks = length(my_seq))

# Convert to plotly

water_table_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(water_table_l1,
                                 water_table_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

NA's in water table likely mean that no water table was observed. NA's are not represented in the figures above. Please observe the table below to understand the relative acreage classified as NA. 
```{r}
if (all(is.na(EDIT_updates_filtered$water_table_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# components NA' = sum(is.na(EDIT_updates_filtered$water_table_l_in)),
             '% components NA' = 100*sum(is.na(EDIT_updates_filtered$water_table_l_in))/length(EDIT_updates_filtered$water_table_l_in),
             check.names = FALSE)

}

```

**Water table low acreage weighted 5th percentile** -  use as 'low (actual)'. This value may be NA.  

The percentile calculation for water table is performed differently because NAs in the water table likely mean no water table was observed. For that reason, we will consider NAs to be an ordinal value that is deeper than all other quantitative values.
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$water_table_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
    water_table_na_l_rem_expanded <- rep(EDIT_updates_filtered$water_table_l_in, EDIT_updates_filtered$comp_acres/comp_acreage_unit) |> sort(na.last = TRUE)

    data.frame("5%" = water_table_na_l_rem_expanded[ceiling(0.05 * length(water_table_na_l_rem_expanded))], 
             check.names = FALSE)
  
}

```

**Water table low acreage weighted summary statistics** - Summary statistics are calculated using numeric values (NA values omitted). The combined acreage of components with an NA value is listed.
```{r}

if (all(is.na(EDIT_updates_filtered$water_table_l_in))) {print("No data populated in NASIS for this field. For water table, this likely means no water table was observed.")
  
} else{
  
summary(water_table_na_l_rem_expanded)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$water_table_l_in)) & all(is.na(EDIT_updates_filtered$water_table_h_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

  water_table_r1 <-
    ggplot2::ggplot(water_table.df |>  dplyr::filter(water_table %in% c("water_table_h_in", "water_table_l_in")) |> 
  dplyr::mutate(highlow = ifelse(water_table == "water_table_h_in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(
                      ggplot2::aes(text = paste("coiid(s):", coiid2)),
                      binwidth = bin_width,
                      center = min(my_seq) + bin_width / 2,
                      color = "black",
                      fill = "grey", closed = "left"
                    ) +
    ggplot2::ggtitle("Water table representative") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
    ggplot2::ylab("# of components") + ggplot2::xlab("Water table (inches)") +
    ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                                n.breaks = length(my_seq))

# Convert histogram to plotly

water_table_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

  water_table_r2 <-
    ggplot2::ggplot(water_table.df |>  dplyr::filter(water_table %in% c("water_table_h_in", "water_table_l_in")) |> 
  dplyr::mutate(highlow = ifelse(water_table == "water_table_h_in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(
                      ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres),
                      binwidth = bin_width,
                      center = min(my_seq) + bin_width / 2,
                      color = "black",
                      fill = "grey", closed = "left"
                    ) +
    ggplot2::ggtitle("Water table representative") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
    ggplot2::ylab("# of acres") + ggplot2::xlab("Water table (inches)") +
    ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                                n.breaks = length(my_seq))

# Convert to plotly

water_table_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(water_table_r1,
                                 water_table_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

NA's in water table likely mean that no water table was observed. NA's are not represented in the figures above. Please observe the table below to understand the relative acreage classified as NA. 
```{r}
if (all(is.na(EDIT_updates_filtered$water_table_l_in)) & all(is.na(EDIT_updates_filtered$water_table_h_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# values NA' = sum(is.na(EDIT_updates_filtered$water_table_l_in)) + sum(is.na(EDIT_updates_filtered$water_table_h_in)),
             '% values NA' = 100*(sum(is.na(EDIT_updates_filtered$water_table_h_in)) +
                                        sum(is.na(EDIT_updates_filtered$water_table_l_in)))/
               (length(EDIT_updates_filtered$water_table_l_in) * 2), check.names = FALSE)

}

```

**Water table representative - combined component high/low 20th and 80th percentile** - use as 'low (representative)' & 'high (representative)'. This value may be NA. 

The percentile calculation for water table is performed differently because NAs in the water table likely mean no water table was observed. For that reason, we will consider NAs to be an ordinal value that is deeper than all other quantitative values.
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$water_table_l_in)) & all(is.na(EDIT_updates_filtered$water_table_h_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  water_table_r_sum_stats <- rep(water_table.df$value, plyr::round_any(water_table.df$acres/comp_acreage_unit, accuracy = 1, f = round)) |> sort(na.last = TRUE)

  data.frame("20%" = water_table_r_sum_stats[ceiling(0.2 * length(water_table_r_sum_stats))],
             "80%" = water_table_r_sum_stats[ceiling(0.8 * length(water_table_r_sum_stats))],
             check.names = FALSE)
  
}

```

**Water table representative acreage weighted summary statistics** - Summary statistics are calculated using numeric values (NA values omitted). The combined acreage of components with an NA value is listed.
```{r}

if (all(is.na(EDIT_updates_filtered$water_table_l_in))) {print("No data populated in NASIS for this field. For water table, this likely means no water table was observed.")
  
} else{
  
summary(water_table_r_sum_stats)

}
```

#### **Aspect**

If site is not strongly associated with a particular aspect, select "not applicable"

I would suggest that we change aspect to be northness [cos(aspect)] and eastness [sin(aspect)]. These are more ecologically meaningful and would allow the statistics to be summarized. 

I’m not aware of specific guidance on populating aspect data in EDIT. I personally would choose the aspects that you feel are representative. That would mean trimming values that don’t seem like a core part of the distribution/concept. I am fairly liberal in saying that aspect is ‘not applicable’ if I don’t think I have the data to make a meaningful call. If there are only 2-3 components correlated to the ecosite, I would consider it pretty challenging to make a strong call unless they are tightly clustered. 
```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$aspect_rep))) {print("No data populated in NASIS for this field.")
  
} else{

# Define bin width and round by measures

bin_width <- 30
round_by <- 30

# Default dataframe

aspect.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, aspect_rep, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "aspect",
    values_to = "value",
    cols = c("aspect_rep")) |> 
  dplyr::rename(acres = comp_acres)

# Create component histogram

aspect_r1 <- ggplot2::ggplot(aspect.df |>  dplyr::filter(aspect == "aspect_rep"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Aspect") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Aspect (degrees)") +
  ggplot2::scale_x_continuous(limits = c(0, 360),
  breaks = seq(0,360, 30))

# Convert histogram to plotly

aspect_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$aspect_rep, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$aspect_rep, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$aspect_rep, na.rm = TRUE) &
             min(EDIT_updates_filtered$aspect_rep, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$aspect_rep, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$aspect_rep, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq  
  
}

# Acreage weighted histogram

aspect_r2 <- ggplot2::ggplot(aspect.df |>  dplyr::filter(aspect == "aspect_rep"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Aspect") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Aspect (degrees)") +
  ggplot2::scale_x_continuous(limits = c(0, 360),
  breaks = seq(0,360, 30))

# Convert to plotly

aspect_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(aspect_r1,
                                 aspect_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```


<center><h3>**Climate features tab:** </h3></center>

Populate climate variables in EDIT by selecting appropriate climate stations. If appropriate climate stations are not available, Mean Annual Precip. and Frost Free Days can be accessed from data mapunits.  

#### **Mean annual precipitation**

MAP prep code
```{r, message=FALSE}
if (all(is.na(EDIT_updates_filtered$map_h_in)) & all(is.na(EDIT_updates_filtered$map_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
# Define bin width and round by measures

my_vector <- c(1, 2, 5, 10, 15, 20)

my_factor <- diff(range(c(EDIT_updates_filtered$map_h_in, 
             EDIT_updates_filtered$map_l_in), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

map.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, map_h_in, map_l_in, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "map",
    values_to = "value",
    cols = c("map_h_in", "map_l_in")) |> 
  dplyr::rename(acres = comp_acres)
  
# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$map_l_in, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$map_h_in, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$map_l_in, na.rm = TRUE) &
             min(EDIT_updates_filtered$map_l_in, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$map_l_in, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$map_h_in, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq  
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

}
```


```{r, warning=FALSE, message=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$map_h_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

map_h1 <- ggplot2::ggplot(map.df |>  dplyr::filter(map == "map_h_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Mean Annual Precip. high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Mean Annual Precip (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

map_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

map_h2 <- ggplot2::ggplot(map.df |>  dplyr::filter(map == "map_h_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Mean Annual Precip. high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Mean Annual Precip (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

map_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(map_h1,
                                 map_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**MAP high acreage weighted 95th percentile**  -  use as 'high (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$map_h_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  map_h_stats <- rep(EDIT_updates_filtered$map_h_in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                            accuracy = 1, f = round)) |> sort(na.last = NA)
  
  map_h_stats |> quantile(probs = c(0.95), na.rm = TRUE, type = 1)

}

```

**MAP high acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$map_h_in))) {print("No data populated in NASIS for this field.")
  
} else{

summary(map_h_stats)

}
```

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

if (all(is.na(EDIT_updates_filtered$map_l_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram
  
map_l1 <- ggplot2::ggplot(map.df |>  dplyr::filter(map == "map_l_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Mean Annual Precip. low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Mean Annual Precip. (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

map_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

map_l2 <- ggplot2::ggplot(map.df |>  dplyr::filter(map == "map_l_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Mean Annual Precip. low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Mean Annual Precip. (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

map_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(map_l1,
                                 map_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**MAP low acreage weighted 5th percentile** - use as 'low (actual)'
```{r}

if (all(is.na(EDIT_updates_filtered$map_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  

    map_l_stats <- rep(EDIT_updates_filtered$map_l_in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                            accuracy = 1, f = round)) |> sort(na.last = NA)
    
    map_l_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```

**MAP low acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$map_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  

summary(map_l_stats)

}

```


```{r, warning=FALSE, message=FALSE}
if (all(is.na(EDIT_updates_filtered$map_h_in)) & 
    all(is.na(EDIT_updates_filtered$map_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
# Create component histogram

map_r1 <- ggplot2::ggplot(map.df |>  dplyr::filter(map %in% c("map_h_in", "map_l_in")) |> 
  dplyr::mutate(highlow = ifelse(map == "map_h_in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Mean Annual Precip. representative") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Mean Annual Precip. (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

map_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

map_r2 <- ggplot2::ggplot(map.df |>  dplyr::filter(map %in% c("map_h_in", "map_l_in")) |> 
  dplyr::mutate(highlow = ifelse(map == "map_h_in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Mean Annual Precip. representative") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Mean Annual Precip. (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

map_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(map_r1,
                                 map_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**MAP representative - combined component high/low acreage weighted 20th and 80th percentile** - use as 'low (representative)' & 'high (representative)' 
```{r}
if (all(is.na(EDIT_updates_filtered$map_h_in)) &
    all(is.na(EDIT_updates_filtered$map_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  map_r_stats <- rep(map.df$value, plyr::round_any(map.df$acres/comp_acreage_unit, 
                      accuracy = 1, f = round)) |> sort(na.last = NA)
  
  map_r_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE, type = 1)

}

```

**MAP representative acreage weighted summary statistics**
```{r}
if (all(is.na(EDIT_updates_filtered$map_h_in)) &
    all(is.na(EDIT_updates_filtered$map_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
summary(map_r_stats)

}

```

#### **Frost Free Days**

Frost Free Days prep code
```{r, message=FALSE}

if (all(is.na(EDIT_updates_filtered$ffd_h)) & all(is.na(EDIT_updates_filtered$ffd_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
# Define bin width and round by measures

my_vector <- c(5, 10, 15, 25, 30, 50, 75)

my_factor <- diff(range(c(EDIT_updates_filtered$ffd_h, 
             EDIT_updates_filtered$ffd_l), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

ffd.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, ffd_h, ffd_l, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "ffd",
    values_to = "value",
    cols = c("ffd_h", "ffd_l")) |> 
  dplyr::rename(acres = comp_acres)
  

# Sequence of axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$ffd_l, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$ffd_h, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$ffd_l, na.rm = TRUE) &
             min(EDIT_updates_filtered$ffd_l, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$ffd_l, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$ffd_h, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq     
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

}
```


```{r, warning=FALSE, message=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$ffd_h))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

ffd_h1 <- ggplot2::ggplot(ffd.df |>  dplyr::filter(ffd == "ffd_h"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Frost Free Days high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 45)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Frost Free Days") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

ffd_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

ffd_h2 <- ggplot2::ggplot(ffd.df |>  dplyr::filter(ffd == "ffd_h"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Frost Free Days high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 45)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Frost Free Days") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

ffd_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(ffd_h1,
                                 ffd_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Frost Free Days high acreage weighted 95th percentile**  -  use as 'high (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$ffd_h))) {print("No data populated in NASIS for this field.")
  
} else{
  
  ffd_h_stats <- rep(EDIT_updates_filtered$ffd_h, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                          accuracy = 1, f = round)) |> sort(na.last = NA)
  
  ffd_h_stats |> quantile(probs = c(0.95), na.rm = TRUE, type = 1)

}

```

**Frost Free Days high acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$ffd_h))) {print("No data populated in NASIS for this field.")
  
} else{

summary(ffd_h_stats)
  
}
```

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

if (all(is.na(EDIT_updates_filtered$ffd_l))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

ffd_l1 <- ggplot2::ggplot(ffd.df |>  dplyr::filter(ffd == "ffd_l"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Frost Free Days low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 45)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Frost Free Days") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

ffd_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

ffd_l2 <- ggplot2::ggplot(ffd.df |>  dplyr::filter(ffd == "ffd_l"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Frost Free Days low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 45)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Frost Free Days") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

ffd_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(ffd_l1,
                                 ffd_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Frost Free Days low acreage weighted 5th percentile** - use as 'low (actual)'
```{r}

if (all(is.na(EDIT_updates_filtered$ffd_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
    ffd_l_stats <- rep(EDIT_updates_filtered$ffd_l, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                    accuracy = 1, f = round)) |> sort(na.last = NA)
    
    ffd_l_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```

**Frost Free Days low acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$ffd_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
summary(ffd_l_stats)

}

```

```{r, warning=FALSE, message=FALSE}
if (all(is.na(EDIT_updates_filtered$ffd_h)) &
    all(is.na(EDIT_updates_filtered$ffd_l))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

ffd_r1 <- ggplot2::ggplot(ffd.df |>  dplyr::filter(ffd %in% c("ffd_h", "ffd_l")) |> 
  dplyr::mutate(highlow = ifelse(ffd == "ffd_h", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Frost Free Days representitive") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 45)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Frost Free Days") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

ffd_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

ffd_r2 <- ggplot2::ggplot(ffd.df |>  dplyr::filter(ffd %in% c("ffd_h", "ffd_l")) |> 
  dplyr::mutate(highlow = ifelse(ffd == "ffd_h", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Frost Free Days representitive") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 45)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Frost Free Days") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

ffd_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(ffd_r1,
                                 ffd_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Frost Free Days representative - combined component high/low acreage weighted 20th and 80th percentile** - use as 'low (representative)' & 'high (representative)' 
```{r}
if (all(is.na(EDIT_updates_filtered$ffd_h)) &
    all(is.na(EDIT_updates_filtered$ffd_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
  ffd_r_stats <- rep(ffd.df$value, plyr::round_any(ffd.df$acres/comp_acreage_unit, 
                              accuracy = 1, f = round)) |> sort(na.last = NA)
  
  ffd_r_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE, type = 1)

}

```

**Frost Free Days representative acreage weighted summary statistics**
```{r}
if (all(is.na(EDIT_updates_filtered$ffd_h)) &
    all(is.na(EDIT_updates_filtered$ffd_l))) {print("No data populated in NASIS for this field.")
  
} else{
  
  summary(ffd_r_stats)
  
}

```

<center><h3>**Water features tab:** </h3></center>

Canned statement for both 'Water features narrative' and 'Wetland description':
Due to the topographic position, this site does not have water features or wetlands. 

<center><h3>**Soil features tab:** </h3></center>

#### **Parent material**

**Parent material description frequency**
```{r}
if (all(is.na(EDIT_updates_filtered$parent_material))) {print("No data populated in NASIS for this field.")
  
} else{

    EDIT_updates_filtered |> dplyr::select(parent_material, comp_acres, coiid) |> dplyr::group_by(parent_material) |> 
    dplyr::summarise(acres = sum(comp_acres), `# components` = dplyr::n(), coiids = toString(coiid)) |> 
    dplyr::arrange(dplyr::desc(acres)) |> dplyr::rename(`parent material` = parent_material)|> knitr::kable()
 
}
```

Wordcloud of acreage weighted parent material descriptions. The wordcloud should give an overview of the words used to describe the parent material, please view the full parent material strings above for explicit details. 
```{r, warning=FALSE}

if (all(is.na(EDIT_updates_filtered$parent_material))) {print("No data populated in NASIS for this field.")
  
} else{
    
  parent_material_df <- EDIT_updates_filtered |>  dplyr::select(coiid, parent_material, comp_acres) |>  dplyr::group_by(parent_material) |>  
    dplyr::summarise(acres = sum(comp_acres),  '# components' = dplyr::n()) |>  dplyr::arrange(desc(acres)) |>  as.data.frame()
  
  parent_material_df
  
  acreage_unit <- ifelse(min(parent_material_df$acres[parent_material_df$acres != 0], na.rm = TRUE) < 100, 1, 
               ifelse(min(parent_material_df$acres[parent_material_df$acres != 0], na.rm = TRUE) < 1000, 10, 
                      ifelse(min(parent_material_df$acres[parent_material_df$acres != 0], na.rm = TRUE) < 10000, 100, 
                             ifelse(min(parent_material_df$acres[parent_material_df$acres != 0], na.rm = TRUE) < 
                                      100000, 1000, 10000))))
  
  
 parent_material_words <-  rep(parent_material_df$parent_material, plyr::round_any(parent_material_df$acres/acreage_unit, accuracy = 1, f = round)) |> strsplit(" ") |> unlist() |>   
              stringr::str_replace(",", "")
 
 wordcloud::wordcloud(parent_material_words)

}

```

**Parent material wordcloud table**
```{r, message=FALSE, warning=FALSE}
if (all(is.na(EDIT_updates_filtered$parent_material))) {print("No data populated in NASIS for this field.")
  
} else{
  
  table(parent_material_words) |>  as.data.frame() |>  dplyr::arrange(desc(Freq)) |>  dplyr::rename(`Parent material words` =
                                                                    parent_material_words) |> dplyr::top_n(15)
  
}

```

#### **Surface texture**

```{r, warning=FALSE, message=FALSE}
if (all(is.na(EDIT_updates_filtered$surf_texture))) {print("No data populated in NASIS for this field.")
  
} else{
  
  if (length(unique(EDIT_updates_filtered$surf_texture <= 15))) {
    
    comp.df <- EDIT_updates_filtered |>  dplyr::select(coiid, surf_texture, comp_acres)  
    
# Component histogram 
comp1 <- ggplot2::ggplot(comp.df, ggplot2::aes(surf_texture)) + 
  ggplot2::geom_bar(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface texture") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 45)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Surface texture")

# Convert to plotly
comp1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage histogram
comp2 <- ggplot2::ggplot(comp.df, ggplot2::aes(surf_texture, comp_acres)) +
  ggplot2::geom_col(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface texture") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 45)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Surface texture")

#Convert to plotly
comp2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures
manipulateWidget::combineWidgets(comp1,
                                 comp2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))


  } else {
    
    EDIT_updates_filtered |> dplyr::select(surf_texture, comp_acres, coiid) |> dplyr::group_by(surf_texture) |> 
    dplyr::summarise(acres = sum(comp_acres), `# components` = dplyr::n(), coiids = toString(coiid)) |> 
    dplyr::arrange(dplyr::desc(acres)) |> knitr::kable()
    
  }

}
```


**Full texture string by frequency**
```{r}
if (all(is.na(EDIT_updates_filtered$surf_texture))) {print("No data populated in NASIS for this field.")
  
} else{
  
  EDIT_updates_filtered |>  dplyr::select(coiid, surf_texture, comp_acres) |>  dplyr::group_by(surf_texture) |>  
    dplyr::summarise(acres = sum(comp_acres),  '# components' = dplyr::n()) |>  dplyr::arrange(desc(acres)) |>  as.data.frame()

}
```


#### **Drainage class**

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

EDIT_updates_filtered$drainage_class <- ifelse(EDIT_updates_filtered$drainage_class == "", "none", EDIT_updates_filtered$drainage_class)


if (all(is.na(EDIT_updates_filtered$drainage_class))) {print("No data populated in NASIS for this field.")
  
} else{
  
my_levels <-  c(
    "subaqueous",
    "very poorly",
    "poorly",
    "somewhat poorly",
    "moderately well",
    "well",
    "somewhat excessively",
    "excessively"
  )
  
  drainage.class.df <- EDIT_updates_filtered |>  dplyr::select(coiid, drainage_class, comp_acres)  
  
  drainage.class.df$drainage_class <- factor(drainage.class.df$drainage_class, levels = my_levels)
    
# Component histogram 
drainage.class1 <- ggplot2::ggplot(drainage.class.df, ggplot2::aes(drainage_class)) + 
  ggplot2::geom_bar(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Drainage class") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Drainage class")

# Convert to plotly
drainage.class1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage histogram
drainage.class2 <- ggplot2::ggplot(drainage.class.df, ggplot2::aes(drainage_class, comp_acres)) +
  ggplot2::geom_col(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Drainage class") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Drainage class")

#Convert to plotly
drainage.class2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures
manipulateWidget::combineWidgets(drainage.class1,
                                 drainage.class2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))
  
}

```

**Drainage class table**

```{r}
if (all(is.na(EDIT_updates_filtered$drainage_class))) {print("No data populated in NASIS for this field.")
  
} else{
  
  EDIT_updates_filtered |>  dplyr::select(coiid, drainage_class, comp_acres) |>  dplyr::group_by(drainage_class) |>  
    dplyr::summarise(acres = sum(comp_acres),  '# components' = dplyr::n()) |>  dplyr::arrange(desc(acres)) |>  as.data.frame()

}
```
**Drainage class acreage weighted percentiles**

Actual low (5%) and actual high (95%). Representative low (20%) and representative high (80%).
```{r}
if (all(is.na(EDIT_updates_filtered$drainage_class))) {print("No data populated in NASIS for this field.")
  
} else{
  
    expanded_drainage_class <- rep(EDIT_updates_filtered$drainage_class, 
                                   times = plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit,
                                                           accuracy = 1, f = round))
    
    expanded_drainage_class <- expanded_drainage_class[order(match(expanded_drainage_class, my_levels))]
    
    data.frame("5%" = expanded_drainage_class[ceiling(0.05 * length(expanded_drainage_class))],
               "20%" = expanded_drainage_class[ceiling(0.20 * length(expanded_drainage_class))],
               "80%" = expanded_drainage_class[ceiling(0.80 * length(expanded_drainage_class))],
             "95%" = expanded_drainage_class[ceiling(0.95 * length(expanded_drainage_class))],
             check.names = FALSE)
  
}
```

#### **Permeability class**

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

if (all(is.na(EDIT_updates_filtered$ksat_l_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  my_levels <-  c(
    "impermeable",
    "very slow",
    "slow",
    "moderately slow",
    "moderately rapid",
    "rapid",
    "very rapid"
  )
  
  EDIT_updates_filtered$ksat_l_0_40_in_mod <-
    plyr::mapvalues(
      EDIT_updates_filtered$ksat_l_0_40_in,
      from = c(
        "very low",
        "low",
        "moderately low",
        "moderately high",
        "high",
        "very high"
      ),
      to = c(
        "very slow",
        "slow",
        "moderately slow",
        "moderately rapid",
        "rapid",
        "very rapid"
      )
    )   
  
  perm.class.df <- EDIT_updates_filtered |>  dplyr::select(coiid, ksat_l_0_40_in, comp_acres)  
  
# Component histogram 
perm.class1 <- ggplot2::ggplot(perm.class.df, ggplot2::aes(ksat_l_0_40_in)) + 
  ggplot2::geom_bar(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Permeability class") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Permeability class")

# Convert to plotly
perm.class1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage histogram
perm.class2 <- ggplot2::ggplot(perm.class.df, ggplot2::aes(ksat_l_0_40_in, comp_acres)) +
  ggplot2::geom_col(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Permeability class") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Permeability class")

#Convert to plotly
perm.class2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures
manipulateWidget::combineWidgets(perm.class1,
                                 perm.class2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))
  
}

```

**Permeability table**

```{r}
if (all(is.na(EDIT_updates_filtered$ksat_l_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  EDIT_updates_filtered |>  dplyr::select(coiid, ksat_l_0_40_in_mod, comp_acres) |>  dplyr::group_by(ksat_l_0_40_in_mod) |>  
    dplyr::summarise(acres = sum(comp_acres),  '# components' = dplyr::n()) |>  dplyr::arrange(desc(acres)) |>  as.data.frame()

}
```

**Permeability acreage weighted percentiles**

Actual low (5%) and actual high (95%). Representative low (20%) and representative high (80%).
```{r}
if (all(is.na(EDIT_updates_filtered$ksat_l_0_40_in_mod))) {print("No data populated in NASIS for this field.")
  
} else{
  
  
    expanded_perm_class <- rep(EDIT_updates_filtered$ksat_l_0_40_in_mod, EDIT_updates_filtered$comp_acres)
    expanded_perm_class <- expanded_perm_class[order(match(expanded_perm_class, my_levels))]
    
    data.frame("5%" = expanded_perm_class[ceiling(0.05 * length(expanded_perm_class))],
      "20%" = expanded_perm_class[ceiling(0.20 * length(expanded_perm_class))],
             "80%" = expanded_perm_class[ceiling(0.80 * length(expanded_perm_class))],
      "95%" = expanded_perm_class[ceiling(0.95 * length(expanded_perm_class))],
             check.names = FALSE)
  
  
}
```

#### **Depth to restriction**

Depth to restriction prep code - this is the top depth of the uppermost restriction in the component restriction table. All restriction types qualify. 

Note: Consider whether restriction data is detailed enough to justify using values beyond depth classes. For many situations, it may be best to use depth class limits for restriction depth and soil depth. Depth classes are represented with red vertical lines in figures for reference. 
```{r, message=FALSE}

depth_classes <- c(0, 20, 40, 60, 80)

if (all(is.na(EDIT_updates_filtered$restrict_depth_h_in)) & 
    all(is.na(EDIT_updates_filtered$restrict_depth_l_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Define bin width and round by measures

my_vector <- c(5, 10, 15, 25, 30, 50, 75)

my_factor <- diff(range(c(EDIT_updates_filtered$restrict_depth_h_in, 
             EDIT_updates_filtered$restrict_depth_l_in), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

restriction.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, restrict_depth_h_in, restrict_depth_l_in, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "restriction",
    values_to = "value",
    cols = c("restrict_depth_h_in", "restrict_depth_l_in")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$restrict_depth_l_in, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$restrict_depth_h_in, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$restrict_depth_l_in, na.rm = TRUE) &
             min(EDIT_updates_filtered$restrict_depth_l_in, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$restrict_depth_l_in, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$restrict_depth_h_in, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq  
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

}

```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$restrict_depth_h_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

restriction_h1 <- ggplot2::ggplot(restriction.df |>  dplyr::filter(restriction == "restrict_depth_h_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey") + 
  ggplot2::ggtitle("Restriction depth high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Restriction depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")


# Convert histogram to plotly


restriction_h1 <- plotly::ggplotly() |> plotly::layout(title = list(text = paste0('Restriction depth high',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Acreage weighted histogram

restriction_h2 <- ggplot2::ggplot(restriction.df |>  dplyr::filter(restriction == "restrict_depth_h_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey") + 
  ggplot2::ggtitle("Restriction depth high") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Restriction depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")




# Convert to plotly

restriction_h2 <- plotly::ggplotly() |> plotly::layout(title = list(text = paste0('Restriction depth high',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Join both figures

manipulateWidget::combineWidgets(restriction_h1,
                                 restriction_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Restriction depth high table**

NAs in restriction depth likely mean that no restriction was observed. NA's are not represented in the figures above. Please observe the table below to understand the relative acreage classified as NA.
```{r}
if (all(is.na(EDIT_updates_filtered$restrict_depth_h_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# components NA' = sum(is.na(EDIT_updates_filtered$restrict_depth_h_in)),
             '% components NA' = 100*sum(is.na(EDIT_updates_filtered$restrict_depth_h_in))/
               length(EDIT_updates_filtered$restrict_depth_h_in),
             check.names = FALSE)
  
  }
```

**Restriction high acreage weighted 95th percentile** -  use as 'high (actual)'. This value may be NA.  

The percentile calculation for restriction is performed differently because NAs in the restriction likely mean no restriction was observed. For that reason, we will consider NAs to be an ordinal value that is deeper than all other quantitative values.
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$restrict_depth_h_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
    restriction_na_h_rem_expanded <- c(rep(EDIT_updates_filtered$restrict_depth_h_in,
                                             EDIT_updates_filtered$comp_acres/comp_acreage_unit)) |> sort(na.last = TRUE)

    data.frame("95%" = restriction_na_h_rem_expanded[ceiling(0.95 * length(restriction_na_h_rem_expanded))], 
             check.names = FALSE)
  
}

```

**Restriction high acreage weighted summary statistics** - Summary statistics are calculated using numeric values (NA values omitted). The combined acreage of components with an NA value is listed.
```{r}

if (all(is.na(EDIT_updates_filtered$restrict_depth_h_in))) {print("No data populated in NASIS for this field. For restriction, this likely means no restriction was observed.")
  
} else{

  summary(restriction_na_h_rem_expanded)
  
}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$restrict_depth_l_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

restriction_l1 <- ggplot2::ggplot(restriction.df |>  dplyr::filter(restriction == "restrict_depth_l_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey") + 
  ggplot2::ggtitle("Restriction depth low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Restriction depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")


# Convert histogram to plotly

restriction_l1 <- plotly::ggplotly() |>  plotly::layout(title = list(text = paste0('Restriction depth low',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Acreage weighted histogram

restriction_l2 <- ggplot2::ggplot(restriction.df |>  dplyr::filter(restriction == "restrict_depth_l_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey") + 
  ggplot2::ggtitle("Restriction depth low") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Restriction depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")

# Convert to plotly

restriction_l2 <- plotly::ggplotly() |>  plotly::layout(title = list(text = paste0('Restriction depth low',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Join both figures

manipulateWidget::combineWidgets(restriction_l1,
                                 restriction_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Restriction depth low table**

NAs in restriction likely mean that no restriction was observed. NA's are not represented in the figures above. Please observe the table below to understand the acreage and number of component associated with a restriction of NA. 
```{r}
if (all(is.na(EDIT_updates_filtered$restrict_depth_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# components NA' = sum(is.na(EDIT_updates_filtered$restrict_depth_l_in)),
             '% components NA' = 100*sum(is.na(EDIT_updates_filtered$restrict_depth_l_in))/
               length(EDIT_updates_filtered$restrict_depth_l_in),
             check.names = FALSE)
  
}
```

**Restriction low: acreage weighted 5th percentile** -  use as 'low (actual)'. This value may be NA. 

The percentile calculation for restriction is performed differently because NAs in the restriction likely mean no restriction was observed. For that reason, we will consider NAs to be an ordinal value that is deeper than all other quantitative values.
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$restrict_depth_l_in))) {print("No data populated in NASIS for this field.")
  
} else{

    restriction_na_l_rem_expanded <- rep(EDIT_updates_filtered$restrict_depth_l_in, EDIT_updates_filtered$comp_acres/comp_acreage_unit)           |> sort(na.last = TRUE)

    data.frame("5%" = restriction_na_l_rem_expanded[ceiling(0.05 * length(restriction_na_l_rem_expanded))], 
             check.names = FALSE)
  
}

```

**Restriction low acreage weighted summary statistics** - Summary statistics are calculated using numeric values (NA values omitted). The combined acreage of components with an NA value is listed.
```{r}

if (all(is.na(EDIT_updates_filtered$restrict_depth_l_in))) {print("No data populated in NASIS for this field. For restriction, this likely means no restriction was observed.")
  
} else{
  
summary(restriction_na_l_rem_expanded)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$restrict_depth_h_in)) &
    all(is.na(EDIT_updates_filtered$restrict_depth_l_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

restriction_r1 <- ggplot2::ggplot(restriction.df |>  dplyr::filter(restriction %in% c("restrict_depth_h_in", "restrict_depth_l_in")) |> 
  dplyr::mutate(highlow = ifelse(restriction == "restrict_depth_h_in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey") + 
  ggplot2::ggtitle("Restriction depth representative") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Restriction depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")

# Convert histogram to plotly

restriction_r1 <- plotly::ggplotly() |>  plotly::layout(title = list(text = paste0('Restriction depth representative',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Acreage weighted histogram

restriction_r2 <- ggplot2::ggplot(restriction.df |>  dplyr::filter(restriction %in% c("restrict_depth_h_in", "restrict_depth_l_in")) |> 
  dplyr::mutate(highlow = ifelse(restriction == "restrict_depth_h_in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey") + 
  ggplot2::ggtitle("Restriction depth representative") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Restriction depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
                              breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")


# Convert to plotly

restriction_r2 <- plotly::ggplotly() |>  plotly::layout(title = list(text = paste0('Restriction depth representative',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Join both figures

manipulateWidget::combineWidgets(restriction_r1,
                                 restriction_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Restriction depth representative table**

NAs in restriction likely mean that no restriction was observed. NA's are not represented in the figures above. Please observe the table below to understand the acreage and number of component associated with a restriction of NA. 
```{r}
if (all(is.na(EDIT_updates_filtered$restrict_depth_h_in)) &
    all(is.na(EDIT_updates_filtered$restrict_depth_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# values NA' = sum(is.na(EDIT_updates_filtered$restrict_depth_h_in)) +
               sum(is.na(EDIT_updates_filtered$restrict_depth_l_in)),
             '% values NA' = 100*(sum(is.na(EDIT_updates_filtered$restrict_depth_h_in)) +
                                        sum(is.na(EDIT_updates_filtered$restrict_depth_l_in)))/
               (length(EDIT_updates_filtered$restrict_depth_h_in) * 2), check.names = FALSE)
  
}
```


**Restriction representative - combined component high/low acreage weighted 20th and 80th percentile** -  use as 'low (representative)' and 'high (representative)'. These values may be NA.  

The percentile calculation for restriction is performed differently because NAs in the restriction likely mean no restriction was observed. For that reason, we will consider NAs to be an ordinal value that is deeper than all other quantitative values.
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$restrict_depth_h_in)) &
    all(is.na(EDIT_updates_filtered$restrict_depth_l_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  restriction_r_sum_stats <- rep(restriction.df$value, plyr::round_any(restriction.df$acres/comp_acreage_unit, accuracy = 1, f = round)) |> sort(na.last = TRUE)

  data.frame("20%" = restriction_r_sum_stats[ceiling(0.2 * length(restriction_r_sum_stats))],
             "80%" = restriction_r_sum_stats[ceiling(0.8 * length(restriction_r_sum_stats))],
             check.names = FALSE)
  
}

```

**Restriction representative acreage weighted summary statistics** - Summary statistics are calculated using numeric values (NA values omitted). The combined acreage of components with an NA value is listed.
```{r}

if (all(is.na(EDIT_updates_filtered$restrict_depth_h_in)) &
    all(is.na(EDIT_updates_filtered$restrict_depth_l_in))) {print("No data populated in NASIS for this field. For restriction, this likely means no restriction was observed.")
  
} else{
  
summary(restriction_r_sum_stats)

}
```

#### **Restriction kind**

EDIT does not ask for restriction kind, but it is useful to know for ecological reasons. Additionally, in some situations the depth to restriction is different than the soil depth - see discussion in the next section, "Soil Depth." 

```{r, warning=FALSE, message=FALSE}
if (all(is.na(EDIT_updates_filtered$restrict_kind))) {print("No data populated in NASIS for this field.")
  
} else{
  
  if(all(EDIT_updates_filtered$restrict_kind == "")) {print("Restriction type blank for all components. This likely means no restriction was observed.")}
  else{
  
  
    
  restrict.kind.df <- EDIT_updates_filtered |>  dplyr::select(coiid, restrict_kind, comp_acres)  
  
  EDIT_updates_filtered |> dplyr::select(restrict_kind, comp_acres, coiid) |> dplyr::group_by(restrict_kind) |> 
    dplyr::summarise(acres = sum(comp_acres), `# components` = dplyr::n(), coiids = toString(coiid)) |> 
    dplyr::arrange(dplyr::desc(acres)) |> knitr::kable()


  }
  
}
```

#### **Soil depth**

For most situations, soil depth is the same as restriction depth. In some situations, the depth to restriction is different than the soil depth. For example, permafrost, abrupt textural changes, and salts constitute restrictions beyond which soil would continue. In these cases, restriction depth and soil depth would differ. In contrast, soil would end at a bedrock restriction and therefore restriction depth and soil depth would be the same. 

The vast majority of NASIS component horizons have RV depth values but do not have high and low depth values. Given that this report relies on component high/low values (see Important Discussion section at the beginning of the document for details), soil depth cannot be summarized in the same way as other sections in this document. 

If you have reason to believe that soil depth is different than restriction depth, you will have to use your own judgement in deciding how to utilize the information summarized for this section. The component high/low section below will report high/low values, but the vast majority of values (if not all) will be missing. Representative values will not be reported as the combined component high/low as is the case in all other sections of this report. Instead, component RV values will be reported in the Soil depth RV section. No standardized values will be calculated for populating EDIT. 

Soil depth prep code
```{r, message=FALSE}
if(all(is.na(EDIT_updates_filtered$soil_depth_l_in)) &
   all(is.na(EDIT_updates_filtered$soil_depth_h_in)) &
   all(is.na(EDIT_updates_filtered$soil_depth_r_in))) {print("No data populated in NASIS for this field.")
  
} else {
  if(all(is.na(EDIT_updates_filtered$soil_depth_l_in)) &
   all(is.na(EDIT_updates_filtered$soil_depth_h_in))) {
    
    
# Define bin width and round by measures
  
my_vector <- c(1, 5, 10, 15, 20, 25, 40, 50)

my_factor <- diff(range(EDIT_updates_filtered$soil_depth_r_in, na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe
  
soil.depth.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, soil_depth_r_in, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "soil.depth",
    values_to = "value",
    cols = c("soil_depth_r_in")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$soil_depth_r_in, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$soil_depth_r_in, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$soil_depth_r_in, na.rm = TRUE) &
             min(EDIT_updates_filtered$soil_depth_r_in, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$soil_depth_r_in, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$soil_depth_r_in, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq  
  
}



my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

    } else {
      
    
    # Define bin width and round by measures
  
my_vector <- c(1, 5, 10, 15, 20, 25, 40, 50)

my_factor <- diff(range(c(EDIT_updates_filtered$soil_depth_h_in, 
             EDIT_updates_filtered$soil_depth_l_in, EDIT_updates_filtered$soil_depth_r_in), 
             na.rm = TRUE))/my_vector



reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe
  
soil.depth.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, soil_depth_h_in, soil_depth_l_in,  soil_depth_r_in,
                                          comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "soil.depth",
    values_to = "value",
    cols = c("soil_depth_h_in", "soil_depth_l_in", "soil_depth_r_in")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes

my_seq <- seq(plyr::round_any(min(c(EDIT_updates_filtered$soil_depth_l_in, 
                                    EDIT_updates_filtered$soil_depth_r_in), na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(c(EDIT_updates_filtered$soil_depth_h_in, EDIT_updates_filtered$soil_depth_r_in), 
                      na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(c(EDIT_updates_filtered$soil_depth_l_in, 
                                    EDIT_updates_filtered$soil_depth_r_in), na.rm = TRUE) &
             min(c(EDIT_updates_filtered$soil_depth_l_in, EDIT_updates_filtered$soil_depth_r_in), 
                 na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(c(EDIT_updates_filtered$soil_depth_l_in, EDIT_updates_filtered$soil_depth_r_in), 
                          na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(c(EDIT_updates_filtered$soil_depth_h_in, 
                                    EDIT_updates_filtered$soil_depth_r_in), na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq  
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}
    
  }
}

```

Soil depth high
```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$soil_depth_h_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

soil_depth1 <- ggplot2::ggplot(soil.depth.df |>  dplyr::filter(soil.depth ==  "soil_depth_h_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Soil depth high") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Soil depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")

# Convert histogram to plotly

soil_depth1 <- plotly::ggplotly() |>  plotly::layout(title = list(text = paste0('Soil depth high',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Acreage weighted histogram

soil_depth2 <- ggplot2::ggplot(soil.depth.df |>  dplyr::filter(soil.depth ==  "soil_depth_h_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Soil depth high") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Soil depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")

# Convert to plotly

soil_depth2 <- plotly::ggplotly() |>  plotly::layout(title = list(text = paste0('Soil depth high',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Join both figures

manipulateWidget::combineWidgets(soil_depth1,
                                 soil_depth2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Soil depth high acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$soil_depth_h_in))) {print("No data populated in NASIS for this field.")

} else{
  
  soil_depth_stats <- rep(EDIT_updates_filtered$soil_depth_h_in,
                             plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, accuracy = 1, 
                        f = round)) |> sort(na.last = NA)
    

  summary(soil_depth_stats)

}
```

Soil depth low
```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$soil_depth_l_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

soil_depth1 <- ggplot2::ggplot(soil.depth.df |>  dplyr::filter(soil.depth ==  "soil_depth_l_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Soil depth low") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Soil depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")

# Convert histogram to plotly

soil_depth1 <- plotly::ggplotly() |>  plotly::layout(title = list(text = paste0('Soil depth low',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Acreage weighted histogram

soil_depth2 <- ggplot2::ggplot(soil.depth.df |>  dplyr::filter(soil.depth ==  "soil_depth_l_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Soil depth low") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Soil depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")

# Convert to plotly

soil_depth2 <- plotly::ggplotly() |>  plotly::layout(title = list(text = paste0('Soil depth low',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Join both figures

manipulateWidget::combineWidgets(soil_depth1,
                                 soil_depth2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Soil depth low acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$soil_depth_l_in))) {print("No data populated in NASIS for this field.")

} else{
  
  soil_depth_stats <- rep(EDIT_updates_filtered$soil_depth_l_in,
                             plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, accuracy = 1, 
                        f = round)) |> sort(na.last = NA)
    

  summary(soil_depth_stats)

}
```

**Soil depth RV** - <span style = "color:red"> Make sure to read the Soil Depth section description above. Soil Depth calculations differ from other sections. </span>
```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$soil_depth_r_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

soil_depth1 <- ggplot2::ggplot(soil.depth.df |>  dplyr::filter(soil.depth ==  "soil_depth_r_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Soil depth RV") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Soil depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")

# Convert histogram to plotly

soil_depth1 <- plotly::ggplotly() |>  plotly::layout(title = list(text = paste0('Soil depth representative',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Acreage weighted histogram

soil_depth2 <- ggplot2::ggplot(soil.depth.df |>  dplyr::filter(soil.depth ==  "soil_depth_r_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Soil depth RV") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Soil depth (inches)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)]) + ggplot2::geom_vline(xintercept = depth_classes, colour="red")

# Convert to plotly

soil_depth2 <- plotly::ggplotly() |>  plotly::layout(title = list(text = paste0('Soil depth representative',
                                    '<br>',
                                    '<sup>',
                                    'Red line(s) represent depth class boundaries',
                                    '</sup>')))

# Join both figures

manipulateWidget::combineWidgets(soil_depth1,
                                 soil_depth2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Soil depth RV acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$soil_depth_r_in))) {print("No data populated in NASIS for this field.")

} else{
  
  soil_depth_stats <- rep(EDIT_updates_filtered$soil_depth_r_in,
                             plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, accuracy = 1, 
                        f = round)) |> sort(na.last = NA)
    

  summary(soil_depth_stats)

}
```

#### **Surface frag (< 3 inches)**

**Please consider the values provided for surface fragments to be estimates.** There is not an exact calculation for surface fragments. The primary challenge is that fragments of different sizes/shapes all have low and high values. If you add the high values together, you will likely end up with an unrealistically high value - possibly over 100%. A method more restrictive than simply adding is used but high values can still be unrealistically high. Please do your own research for populating this field. 

Surface frag (< 3 inches) prep code
```{r, message=FALSE}

# NAs in surface frag generally means zero surface frags

EDIT_updates_filtered$frags_hi_2mm_3in_on_surf <- ifelse(is.na(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf), 0, EDIT_updates_filtered$frags_hi_2mm_3in_on_surf)

EDIT_updates_filtered$frags_lo_2mm_3in_on_surf <- ifelse(is.na(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf), 0, EDIT_updates_filtered$frags_lo_2mm_3in_on_surf)

# Is variable NA?

if (all(is.na(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{

# Define bin width and round by measures

my_vector <- c(1, 5, 10, 20, 25)

my_factor <- diff(range(c(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf, 
             EDIT_updates_filtered$frags_lo_2mm_3in_on_surf), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

surf_frag.df <-
  EDIT_updates_filtered |>  
  dplyr::select(coiid, frags_hi_2mm_3in_on_surf, frags_lo_2mm_3in_on_surf, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "surf_frag",
    values_to = "value",
    cols = c("frags_hi_2mm_3in_on_surf", "frags_lo_2mm_3in_on_surf")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf, na.rm = TRUE) &
                   min(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf, na.rm = TRUE), round_by, f = floor) - bin_width,
  plyr::round_any(max(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
         my_seq            
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

}

```

Surface frag high
```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

surf_frag1 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag == "frags_hi_2mm_3in_on_surf"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag high (< 3 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

surf_frag1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

surf_frag2 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag == "frags_hi_2mm_3in_on_surf"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag high (< 3 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

surf_frag2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(surf_frag1,
                                 surf_frag2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Surface frag high (< 3 inches) acreage weighted 95th percentile** - use as 'high (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{
  
  surf_frag_stats <- rep(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf, 
                        plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                        accuracy = 1, f = round)) |> sort(na.last = NA)
  
  surf_frag_stats |> quantile(probs = c(0.95), na.rm = TRUE, type = 1)


}

```

**Surface frag high (< 3 inches) acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf))) {print("No data populated in NASIS for this field.")

} else{

  summary(surf_frag_stats)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

surf_frag1 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag == "frags_lo_2mm_3in_on_surf"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag low (< 3 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

surf_frag1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

surf_frag2 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag == "frags_lo_2mm_3in_on_surf"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag low (< 3 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

surf_frag2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(surf_frag1,
                                 surf_frag2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Surface frag low (< 3 inches) acreage weighted 5th percentile** - use as 'low (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{
  
    surf_frag_stats <- rep(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, accuracy = 1, 
                        f = round)) |> sort(na.last = NA)
    
    surf_frag_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```

**Surface frag low (< 3 inches) acreage weighted summary statistics**

```{r}

if (all(is.na(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf))) {print("No data populated in NASIS for this field.")

} else{

  summary(surf_frag_stats)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf)) &
    all(is.na(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

surf_frag1 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag %in% c("frags_hi_2mm_3in_on_surf", "frags_lo_2mm_3in_on_surf")) |> 
  dplyr::mutate(highlow = ifelse(surf_frag == "frags_hi_2mm_3in_on_surf", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")),
                         ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag rep. (< 3 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

surf_frag1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

surf_frag2 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag %in% c("frags_hi_2mm_3in_on_surf", "frags_lo_2mm_3in_on_surf")) |> 
  dplyr::mutate(highlow = ifelse(surf_frag == "frags_hi_2mm_3in_on_surf", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")),
                         ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag rep. (< 3 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

surf_frag2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(surf_frag1,
                                 surf_frag2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Surface frag representative (< 3 inches) acreage weighted 20th and 80th percentiles** - use as 'low (representative)' and 'high (representative)' 
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf)) &
    all(is.na(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{
  
  surf_frag_stats <- rep(surf_frag.df$value, plyr::round_any(surf_frag.df$acres/comp_acreage_unit, 
                            accuracy = 1, f = round)) |> sort(na.last = NA)
  
  surf_frag_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE, type = 1)

}

```

**Surface frag representative (< 3 inches) acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$frags_lo_2mm_3in_on_surf)) &
    all(is.na(EDIT_updates_filtered$frags_hi_2mm_3in_on_surf))) {print("No data populated in NASIS for this field.")

} else{

  summary(surf_frag_stats)

}
```

#### **Surface frag (3+ inches)** 

**Please consider the values provided for surface fragments to be estimates.** There is not an exact calculation for surface fragments. The primary challenge is that fragments of different sizes/shapes all have low and high values. If you add the high values together, you will likely end up with an unrealistically high value - possibly over 100%. A method more restrictive than simply adding is used but high values can still be unrealistically high. Please do your own research for populating this field. 

Surface frag (3+ inches) prep code
```{r, message=FALSE}
# NAs in surface frag generally mean zero surface frags

EDIT_updates_filtered$frags_hi_gt_3in_on_surf <- ifelse(is.na(EDIT_updates_filtered$frags_hi_gt_3in_on_surf), 0, EDIT_updates_filtered$frags_hi_gt_3in_on_surf)

EDIT_updates_filtered$frags_lo_gt_3in_on_surf <- ifelse(is.na(EDIT_updates_filtered$frags_lo_gt_3in_on_surf), 0, EDIT_updates_filtered$frags_lo_gt_3in_on_surf)

# Is variable NA?

if (all(is.na(EDIT_updates_filtered$frags_hi_gt_3in_on_surf)) & 
    all(is.na(EDIT_updates_filtered$frags_lo_gt_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{
  
# Define bin width and round by measures

my_vector <- c(1, 5, 10, 15, 25)

my_factor <- diff(range(c(EDIT_updates_filtered$frags_hi_gt_3in_on_surf, 
             EDIT_updates_filtered$frags_lo_gt_3in_on_surf), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

surf_frag.df <-
  EDIT_updates_filtered |>
  dplyr::select(coiid,
                frags_lo_gt_3in_on_surf,
                frags_hi_gt_3in_on_surf,
                comp_acres) |>
  tidyr::pivot_longer(
    names_to = "surf_frag",
    values_to = "value",
    cols = c("frags_lo_gt_3in_on_surf", "frags_hi_gt_3in_on_surf")
  ) |>
  dplyr::rename(acres = comp_acres)

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$frags_lo_gt_3in_on_surf, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$frags_hi_gt_3in_on_surf, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$frags_lo_gt_3in_on_surf, na.rm = TRUE) &
                   min(EDIT_updates_filtered$frags_lo_gt_3in_on_surf, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$frags_lo_gt_3in_on_surf, na.rm = TRUE), round_by, f = floor) - bin_width,
  plyr::round_any(max(EDIT_updates_filtered$frags_hi_gt_3in_on_surf, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
         my_seq            
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}
  
}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$frags_hi_gt_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

surf_frag1 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag == "frags_hi_gt_3in_on_surf"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag high 3 inches +") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

surf_frag1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

surf_frag2 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag == "frags_hi_gt_3in_on_surf"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag high 3 inches +") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

surf_frag2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(surf_frag1,
                                 surf_frag2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Surface frag high (3+ inches) acreage weighted 5th percentile** - use as 'high (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$frags_hi_gt_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{
  
  surf_frag_stats <- rep(EDIT_updates_filtered$frags_hi_gt_3in_on_surf, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                                    accuracy = 1, f = round)) |> sort(na.last = NA)
  
  surf_frag_stats |> quantile(probs = c(0.95), na.rm = TRUE, type = 1)

}

```

**Surface frag high (3+ inches) acreage weighted summary statistics**

```{r}

if (all(is.na(EDIT_updates_filtered$frags_hi_gt_3in_on_surf))) {print("No data populated in NASIS for this field.")

} else{

  summary(surf_frag_stats)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$frags_lo_gt_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

surf_frag1 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag == "frags_lo_gt_3in_on_surf"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag low (3 inches +)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

surf_frag1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

surf_frag2 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag == "frags_lo_gt_3in_on_surf"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag low (3 inches +)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

surf_frag2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(surf_frag1,
                                 surf_frag2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Surface frag low (3+ inches) acreage weighted 5th percentile** - use as 'low (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$frags_lo_gt_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{
  
    surf_frag_stats <- rep(EDIT_updates_filtered$frags_lo_gt_3in_on_surf, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                                    accuracy = 1, f = round)) |> sort(na.last = NA)
    
    surf_frag_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```

**Surface frag low (3+ inches) acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$frags_lo_gt_3in_on_surf))) {print("No data populated in NASIS for this field.")

} else{

  summary(surf_frag_stats)

}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$frags_hi_gt_3in_on_surf)) &
    all(is.na(EDIT_updates_filtered$frags_lo_gt_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

surf_frag1 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag %in% c("frags_hi_gt_3in_on_surf", "frags_lo_gt_3in_on_surf")) |> 
  dplyr::mutate(highlow = ifelse(surf_frag == "frags_hi_gt_3in_on_surf", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")),
                         ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag representative (3 inches +)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

surf_frag1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

surf_frag2 <- ggplot2::ggplot(surf_frag.df |>  dplyr::filter(surf_frag %in% c("frags_hi_gt_3in_on_surf", "frags_lo_gt_3in_on_surf")) |> 
  dplyr::mutate(highlow = ifelse(surf_frag == "frags_hi_gt_3in_on_surf", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")),
                         ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Surface frag representative (3 inches +)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Surface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

surf_frag2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(surf_frag1,
                                 surf_frag2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Surface frag representative (3+ inches) - combined component high/low acreage weighted 20th and 80th percentiles** - use as 'low (representative)' and 'high (representative)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$frags_hi_gt_3in_on_surf)) &
    all(is.na(EDIT_updates_filtered$frags_lo_gt_3in_on_surf))) {print("No data populated in NASIS for this field.")
  
} else{
  
  surf_frag_stats <- rep(surf_frag.df$value, plyr::round_any(elev.df$acres/comp_acreage_unit, accuracy = 1, f = round)) |> 
    sort(na.last = NA)
  
  surf_frag_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE, type = 1)

}

```

**Surface frag representative (3+ inches) acreage weighted summary statistics**

```{r}

if (all(is.na(EDIT_updates_filtered$frags_hi_gt_3in_on_surf)) &
    all(is.na(EDIT_updates_filtered$frags_lo_gt_3in_on_surf))) {print("No data populated in NASIS for this field.")

} else{

  summary(surf_frag_stats)

}
```

#### **Available water capacity (0 - 40 inches)**

AWC (0 - 40 inches) prep code
```{r, message=FALSE}
if (all(is.na(EDIT_updates_filtered$awc_h_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
# Define bin width and round by measures

my_vector <- c(0.25, 0.5, 1, 2, 5, 10)

my_factor <- diff(range(c(EDIT_updates_filtered$awc_h_0_40_in, 
             EDIT_updates_filtered$awc_l_0_40_in), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

awc.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, awc_l_0_40_in, awc_h_0_40_in, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "awc",
    values_to = "value",
    cols = c("awc_l_0_40_in", "awc_h_0_40_in")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$awc_l_0_40_in, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$awc_h_0_40_in, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$awc_l_0_40_in, na.rm = TRUE) &
             min(EDIT_updates_filtered$awc_l_0_40_in, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$awc_l_0_40_in, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$awc_h_0_40_in, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq  
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}
  
}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$awc_h_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

awc_h1 <- ggplot2::ggplot(awc.df |>  dplyr::filter(awc == "awc_h_0_40_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("AWC high (0 - 40 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("AWC") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert histogram to plotly

awc_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

awc_h2 <- ggplot2::ggplot(awc.df |>  dplyr::filter(awc == "awc_h_0_40_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weighted = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("AWC high (0 - 40 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("AWC") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert to plotly

awc_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(awc_h1,
                                 awc_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**AWC high acreage weighted 95th percentile** -  use as 'high (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$awc_h_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  awc_stats <- rep(EDIT_updates_filtered$awc_h_0_40_in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                accuracy = 1, f = round)) |> sort(na.last = NA)
  
  awc_stats |> quantile(probs = c(0.95), na.rm = TRUE, type = 1)

}

```

**AWC high acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$awc_h_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
summary(awc_stats)

}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$awc_l_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

awc_l1 <- ggplot2::ggplot(awc.df |>  dplyr::filter(awc == "awc_l_0_40_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("AWC low (0 - 40 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("AWC") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert histogram to plotly

awc_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

awc_l2 <- ggplot2::ggplot(awc.df |>  dplyr::filter(awc == "awc_l_0_40_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("AWC low (0 - 40 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("AWC") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert to plotly

awc_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(awc_l1,
                                 awc_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**AWC low acreage weighted 5th percentile** -  use as 'low (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$awc_l_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
    awc_stats <- rep(EDIT_updates_filtered$awc_l_0_40_in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                          accuracy = 1, f = round)) |> sort(na.last = NA)
    
    awc_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```

**AWC low acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$awc_l_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
summary(awc_stats)

}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$awc_h_0_40_in)) &
    all(is.na(EDIT_updates_filtered$awc_l_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

awc_r1 <- ggplot2::ggplot(awc.df |>  dplyr::filter(awc %in% c("awc_h_0_40_in", "awc_l_0_40_in")) |> 
  dplyr::mutate(highlow = ifelse(awc == "awc_h_0_40_in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("AWC representative (0 - 40 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("AWC") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert histogram to plotly

awc_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

awc_r2 <- ggplot2::ggplot(awc.df |>  dplyr::filter(awc %in% c("awc_h_0_40_in", "awc_l_0_40_in")) |> 
  dplyr::mutate(highlow = ifelse(awc == "awc_h_0_40_in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("AWC representative (0 - 40 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 35)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("AWC") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert to plotly

awc_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(awc_r1,
                                 awc_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**AWC representative - combined component high/low acreage weighted 20th and 80th percentile** -  use as 'high (representative)' and 'low (representative)'.
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$awc_h_0_40_in)) &
    all(is.na(EDIT_updates_filtered$awc_l_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  awc_stats <- rep(awc.df$value, plyr::round_any(awc.df$acres/comp_acreage_unit, 
                            accuracy = 1, f = round)) |> sort(na.last = NA)
  
  awc_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE, type = 1)

}

```

**AWC representative acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$awc_h_0_40_in)) &
    all(is.na(EDIT_updates_filtered$awc_l_0_40_in))) {print("No data populated in NASIS for this field.")
  
} else{
  
summary(awc_stats)

}
```

#### **pH (0-10 inches)**

pH (0-10 inches) prep code
```{r, message=FALSE}

if (all(is.na(EDIT_updates_filtered$pH_h_0_10in)) &
    all(is.na(EDIT_updates_filtered$pH_l_0_10in))) {print("No data populated in NASIS for this field.")
  
} else{

# Define bin width and round by measures

my_vector <- c(.5, 1, 2, 3)

my_factor <- diff(range(c(EDIT_updates_filtered$pH_h_0_10in, 
             EDIT_updates_filtered$pH_l_0_10in), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

ph.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, pH_h_0_10in, pH_l_0_10in, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "ph",
    values_to = "value",
    cols = c("pH_l_0_10in", "pH_h_0_10in")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$pH_l_0_10in, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$pH_h_0_10in, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$pH_l_0_10in, na.rm = TRUE) &
             min(EDIT_updates_filtered$pH_l_0_10in, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$pH_l_0_10in, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$pH_h_0_10in, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq 
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

}

```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$pH_h_0_10in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

ph_h1 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph == "pH_h_0_10in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("pH high (0 - 10 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert histogram to plotly

ph_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

ph_h2 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph == "pH_h_0_10in"), ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("pH high (0 - 10 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

ph_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(ph_h1,
                                 ph_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**pH high acreage weighted 95th percentile** -  use as 'high (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$pH_h_0_10in))) {print("No data populated in NASIS for this field.")
  
} else{
  

  ph_stats <- rep(EDIT_updates_filtered$pH_h_0_10in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                accuracy = 1, f = round)) |> sort(na.last = NA)
  
  ph_stats |> quantile(probs = c(0.95), na.rm = TRUE, type = 1)

}

```

**pH high acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$pH_h_0_10in))) {print("No data populated in NASIS for this field.")
  
} else{
  
summary(ph_stats)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$pH_l_0_10in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

ph_l1 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph == "pH_l_0_10in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("pH low (0 - 10 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert histogram to plotly

ph_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

ph_l2 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph == "pH_l_0_10in"), ggplot2::aes(x = value)) +
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("pH low (0 - 10 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert to plotly

ph_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(ph_l1,
                                 ph_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**pH low acreage weighted 5th percentile** -  use as 'low (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$pH_l_0_10in))) {print("No data populated in NASIS for this field.")
  
} else{
  
    ph_stats <- rep(EDIT_updates_filtered$pH_l_0_10in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                          accuracy = 1, f = round)) |> sort(na.last = NA)
    
    ph_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```

**pH low acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$pH_l_0_10in))) {print("No data populated in NASIS for this field.")
  
} else{
  
summary(ph_stats)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$pH_l_0_10in)) &
    all(is.na(EDIT_updates_filtered$pH_h_0_10in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

ph_r1 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph %in% c("pH_h_0_10in", "pH_l_0_10in")) |> 
  dplyr::mutate(highlow = ifelse(ph == "pH_h_0_10in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("pH representative (0 - 10 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

ph_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

ph_r2 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph %in% c("pH_h_0_10in", "pH_l_0_10in")) |> 
  dplyr::mutate(highlow = ifelse(ph == "pH_h_0_10in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) +
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("pH representative (0 - 10 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])


# Convert to plotly

ph_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(ph_r1,
                                 ph_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```


**pH representative - combined component high/low acreage weighted 20th and 80th percentile** - use as 'low (representative)' & 'high (representative)' 
```{r}
if (all(is.na(EDIT_updates_filtered$pH_l_0_10in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  ph_stats <- rep(ph.df$value, plyr::round_any(ph.df$acres/comp_acreage_unit, 
                            accuracy = 1, f = round)) |> sort(na.last = NA)
  
  ph_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE, type = 1)
  
}
```

**pH representative acreage weighted summary statistics**
```{r}
if (all(is.na(EDIT_updates_filtered$pH_l_0_10in))) {print("No data populated in NASIS for this field.")
  
} else{
  
summary(ph_stats)
  
}
```


#### **pH (10-40 inches)**

pH (10 - 40 inches) prep code
```{r, message=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$pH_h_10_40in)) &
    all(is.na(EDIT_updates_filtered$pH_l_10_40in))) {print("No data populated in NASIS for this field.")

} else{
  
# Define bin width and round by measures

my_vector <- c(.5, 1, 2, 3)

my_factor <- diff(range(c(EDIT_updates_filtered$pH_h_10_40in, 
             EDIT_updates_filtered$pH_l_10_40in), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

ph.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, pH_h_10_40in, pH_l_10_40in, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "ph",
    values_to = "value",
    cols = c("pH_l_10_40in", "pH_h_10_40in")) |> 
  dplyr::rename(acres = comp_acres)
  
# Sequence for axes

my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$pH_l_10_40in, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$pH_h_10_40in, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$pH_l_10_40in, na.rm = TRUE) &
             min(EDIT_updates_filtered$pH_l_10_40in, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$pH_l_10_40in, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$pH_h_10_40in, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq  
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}

}

```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$pH_h_10_40in))) {print("No data populated in NASIS for this field.")

} else{

# Create component histogram

ph_h1 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph == "pH_h_10_40in"), ggplot2::aes(x = value)) +
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres),
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") +
  ggplot2::ggtitle("pH high (10 - 40 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

ph_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

ph_h2 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph == "pH_h_10_40in"), ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres),
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") +
  #ggplot2::theme_minimal() +
  ggplot2::ggtitle("pH high (10 - 40 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

ph_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(ph_h1,
                                 ph_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

NAs in pH (10-40) likely mean that soils are less than 10 inches deep. NA's are not represented in the figures above. Please observe the table below to determine if NAs exist in your dataset.
```{r}
if (all(is.na(EDIT_updates_filtered$pH_h_10_40in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# components NA' = sum(is.na(EDIT_updates_filtered$pH_h_10_40in)),
             '% components NA' = 100*sum(is.na(EDIT_updates_filtered$pH_h_10_40in))/length(EDIT_updates_filtered$pH_h_10_40in),
             check.names = FALSE)

  
}
```

**pH high (10-40 inches) acreage weighted 95th percentile** -  use as 'high (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$pH_h_10_40in))) {print("No data populated in NASIS for this field.")

} else{

  ph_stats <- rep(EDIT_updates_filtered$pH_h_10_40in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                accuracy = 1, f = round)) |> sort(na.last = NA)
  
  ph_stats |> quantile(probs = c(0.95), na.rm = TRUE, type = 1)

}

```

**pH high acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$pH_h_10_40in))) {print("No data populated in NASIS for this field.")

} else{

summary(ph_stats)
  
}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$pH_l_10_40in))) {print("No data populated in NASIS for this field.")

} else{

# Create component histogram

ph_l1 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph == "pH_l_10_40in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)),
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") +
  #ggplot2::theme_minimal() +
  ggplot2::ggtitle("pH low (10 - 40 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

ph_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

ph_l2 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph == "pH_l_10_40in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)),
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") +
  #ggplot2::theme_minimal() +
  ggplot2::ggtitle("pH low (10 - 40 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

ph_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(ph_l1,
                                 ph_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

NAs in pH (10-40) likely mean that soils are less than 10 inches deep. NA's are not represented in the figures above. Please observe the table below to determine if NAs exist in your dataset.
```{r}
if (all(is.na(EDIT_updates_filtered$pH_l_10_40in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# components NA' = sum(is.na(EDIT_updates_filtered$pH_l_10_40in)),
             '% components NA' = 100*sum(is.na(EDIT_updates_filtered$pH_l_10_40in))/length(EDIT_updates_filtered$pH_l_10_40in),
             check.names = FALSE)
}

```

**pH low (10-40 inches) acreage weighted 5th percentile** -  use as 'low (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$pH_l_10_40in))) {print("No data populated in NASIS for this field.")

} else{

    ph_stats <- rep(EDIT_updates_filtered$pH_l_10_40in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                          accuracy = 1, f = round)) |> sort(na.last = NA)
    
    ph_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```

**pH low (10-40 inches) acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$pH_l_10_40in))) {print("No data populated in NASIS for this field.")

} else{

summary(ph_stats)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$pH_l_10_40in)) &
    all(is.na(EDIT_updates_filtered$pH_h_10_40in))) {print("No data populated in NASIS for this field.")

} else{

# Create component histogram

ph_r1 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph %in% c("pH_h_10_40in", "pH_l_10_40in")) |> 
  dplyr::mutate(highlow = ifelse(ph == "pH_h_10_40in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)),
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") +
  ggplot2::ggtitle("pH representative (10 - 40 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

ph_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

ph_r2 <- ggplot2::ggplot(ph.df |>  dplyr::filter(ph %in% c("pH_h_10_40in", "pH_l_10_40in")) |> 
  dplyr::mutate(highlow = ifelse(ph == "pH_h_10_40in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres),
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") +
  ggplot2::ggtitle("pH representative (10 - 40 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("pH") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

ph_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(ph_r1,
                                 ph_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

NAs in pH (10-40) likely mean that soils are less than 10 inches deep. NA's are not represented in the figures above. Please observe the table below to determine if NAs exist in your dataset.
```{r}
if (all(is.na(EDIT_updates_filtered$pH_l_10_40in)) &
    all(is.na(EDIT_updates_filtered$pH_h_10_40in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  data.frame('# values NA' = sum(is.na(EDIT_updates_filtered$pH_l_10_40in)) + sum(is.na(EDIT_updates_filtered$pH_h_10_40in)),
             '% values NA' = 100*(sum(is.na(EDIT_updates_filtered$pH_l_10_40in)) + sum(is.na(EDIT_updates_filtered$pH_h_10_40in)))/
               (length(EDIT_updates_filtered$pH_l_10_40in) * 2), check.names = FALSE)
  
}

```


**pH (10-40 inches) representative - combined component high/low 20th and 80th percentile** - use as 'low (representative)' & 'high (representative)' 
```{r}
if (all(is.na(EDIT_updates_filtered$pH_l_10_40in)) &
    all(is.na(EDIT_updates_filtered$pH_h_10_40in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  ph_stats <- rep(ph.df$value, plyr::round_any(ph.df$acres/comp_acreage_unit, 
                            accuracy = 1, f = round)) |> sort(na.last = NA)
  
  ph_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE, type = 1)

}

```

**pH (10-40 inches) representative acreage weighted summary statistics**
```{r}
if (all(is.na(EDIT_updates_filtered$pH_l_10_40in)) &
    all(is.na(EDIT_updates_filtered$pH_h_10_40in))) {print("No data populated in NASIS for this field.")
  
} else{
  
summary(ph_stats)

}

```


```{r, message=FALSE, echo=FALSE}
  
# #### **Subsurface fragments top depth (pseudo-organic matter depth)**
# 
# Subsurface fragments top depth is equivalent to organic matter except for Histosols. Organic matter depth is not a field that needs to be populated in EDIT. It is included here for your own information and to be used in determining the depth range that subsurface fragments are being reported for. Rock fragments are not calculated for
# 
# Subsurface fragments top depth


# if(all(is.na(EDIT_updates_filtered$o_horizon_depth_in))) {
#     
#     
# # Define bin width and round by measures
#   
# my_vector <- c(1, 5, 10, 15, 20, 25, 40, 50)
# 
# my_factor <- diff(range(EDIT_updates_filtered$o_horizon_depth_in, na.rm = TRUE))/my_vector
# 
# reducer <- my_factor < 25
# 
# bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]
# 
# round_by <- bin_width
# 
# # Default dataframe
#   
# o.horizon.df <-
#   EDIT_updates_filtered |>  dplyr::select(coiid, o_horizon_depth_in, comp_acres) |> 
#   tidyr::pivot_longer(
#     names_to = "o.horizon",
#     values_to = "value",
#     cols = c("o_horizon_depth_in")) |> 
#   dplyr::rename(acres = comp_acres)
# 
# # Sequence for axes
# 
# my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$o_horizon_depth_in, na.rm = TRUE), round_by, f = floor),
#   plyr::round_any(max(EDIT_updates_filtered$o_horizon_depth_in, na.rm = TRUE), round_by, f = ceiling), bin_width)
# 
# my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$o_horizon_depth_in, na.rm = TRUE) &
#              min(EDIT_updates_filtered$o_horizon_depth_in, na.rm = TRUE) - bin_width >= 0){
#   seq(plyr::round_any(min(EDIT_updates_filtered$o_horizon_depth_in, na.rm = TRUE), round_by, f = floor) - bin_width,
#       plyr::round_any(max(EDIT_updates_filtered$o_horizon_depth_in, na.rm = TRUE), round_by, f = ceiling), bin_width)
# } else{
#   
#   my_seq  
#   
# }
# 
# my_seq <- if(length(my_seq) == 1) {
#   my_seq <- c(my_seq[1], my_seq[1] + bin_width)
# } else {
#   my_seq
# }
#     } else {
#       
#     
#     # Define bin width and round by measures
#   
# my_vector <- c(1, 5, 10, 15, 20, 25, 40, 50)
# 
# my_factor <- diff(range(c(EDIT_updates_filtered$soil_depth_h_in, 
#              EDIT_updates_filtered$soil_depth_l_in, EDIT_updates_filtered$o_horizon_depth_in), 
#              na.rm = TRUE))/my_vector
# 
# 
# 
# reducer <- my_factor < 25
# 
# bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]
# 
# round_by <- bin_width
# 
# # Default dataframe
#   
# o.horizon.df <-
#   EDIT_updates_filtered |>  dplyr::select(coiid, soil_depth_h_in, soil_depth_l_in,  o_horizon_depth_in,
#                                           comp_acres) |> 
#   tidyr::pivot_longer(
#     names_to = "o.horizon",
#     values_to = "value",
#     cols = c("o_horizon_depth_in")) |> 
#   dplyr::rename(acres = comp_acres)
# 
# # Sequence for axes
# 
# my_seq <- seq(plyr::round_any(min(c(EDIT_updates_filtered$soil_depth_l_in, 
#                                     EDIT_updates_filtered$o_horizon_depth_in), na.rm = TRUE), round_by, f = floor),
#   plyr::round_any(max(c(EDIT_updates_filtered$soil_depth_h_in, EDIT_updates_filtered$o_horizon_depth_in), 
#                       na.rm = TRUE), round_by, f = ceiling), bin_width)
# 
# my_seq <- if(min(my_seq) == min(c(EDIT_updates_filtered$soil_depth_l_in, 
#                                     EDIT_updates_filtered$o_horizon_depth_in), na.rm = TRUE) &
#              min(c(EDIT_updates_filtered$soil_depth_l_in, EDIT_updates_filtered$o_horizon_depth_in), 
#                  na.rm = TRUE) - bin_width >= 0){
#   seq(plyr::round_any(min(c(EDIT_updates_filtered$soil_depth_l_in, EDIT_updates_filtered$o_horizon_depth_in), 
#                           na.rm = TRUE), round_by, f = floor) - bin_width,
#       plyr::round_any(max(c(EDIT_updates_filtered$soil_depth_h_in, 
#                                     EDIT_updates_filtered$o_horizon_depth_in), na.rm = TRUE), round_by, f = ceiling), bin_width)
# } else{
#   
#   my_seq  
#   
# }
# 
# my_seq <- if(length(my_seq) == 1) {
#   my_seq <- c(my_seq[1], my_seq[1] + bin_width)
# } else {
#   my_seq
# }
#     
#     
#   }

```

```{r, message=FALSE, warning=FALSE, echo=FALSE}
# Is variable NA?

# if (all(is.na(EDIT_updates_filtered$o_horizon_depth_in))) {print("No data populated in NASIS for this field.")
#   
# } else{

# Create component histogram

# o_horizon1 <- ggplot2::ggplot(o.horizon.df |>  dplyr::filter(o.horizon ==  "o_horizon_depth_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
#   binwidth = bin_width,
#   center = min(my_seq) + bin_width/2,
#   color = "black",
#   fill = "grey", closed = "left") + 
#   ggplot2::ggtitle("Organic horizon(s) depth RV") + 
#   ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
#   ggplot2::ylab("# of components") + ggplot2::xlab("Organic horizon(s) depth (inches)") +
#   ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
#   breaks = my_seq[seq(1, length(my_seq), by = 2)])
# 
# # Convert histogram to plotly
# 
# o_horizon1 <- plotly::ggplotly(tooltip = c("text", "value"))
# 
# # Acreage weighted histogram
# 
# o_horizon2 <- ggplot2::ggplot(o.horizon.df |>  dplyr::filter(o.horizon ==  "o_horizon_depth_in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
#   binwidth = bin_width,
#   center = min(my_seq) + bin_width/2,
#   color = "black",
#   fill = "grey", closed = "left") + 
#   ggplot2::ggtitle("Organic horizon(s) depth RV") + 
#   ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
#   ggplot2::ylab("# of acres") + ggplot2::xlab("Organic horizon(s) depth (inches)") +
#   ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
#   breaks = my_seq[seq(1, length(my_seq), by = 2)])
# 
# # Convert to plotly
# 
# o_horizon2 <- plotly::ggplotly(tooltip = c("text", "value"))
# 
# # Join both figures
# 
# manipulateWidget::combineWidgets(o_horizon1,
#                                  o_horizon2,
#                                  nrow = 1,
#                                  ncol = 2,
#                                  colsize = c(1, 1))

# }

```

**Organic horizon(s) depth RV acreage weighted summary statistics**
```{r, echo=FALSE}

if (all(is.na(EDIT_updates_filtered$o_horizon_depth_in))) {print("No data populated in NASIS for this field.")

} else{
  
  o_horizon_stats <- rep(EDIT_updates_filtered$o_horizon_depth_in,
                             plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, accuracy = 1, 
                        f = round)) |> sort(na.last = NA)
    

  summary(o_horizon_stats)

}
```

#### **Subsurface frag (< 3 inches)**

The calculation of subsurface fragments is more difficult than other properties in this report for several reasons. Firstly, adding component high and low values is not as simple we would like. If you add all the component high values for gravels, cobbles, stones, and boulders, you will very likely end up with an unrealistically high amount of fragments - potentially over 100%. Similarly, adding component low values will likely result in unrealistically low amounts of fragments. For this reason, there is not an exact calculation of subsurface rock fragments. **Please consider the values provided to be estimates.** 

The depth range for subsurface fragments needs to be user-defined. Subsurface frags are calculated from below the O horizon(s) to top of the first restriction, except in Histosols, where rock fragments in the O horizon(s) are included. Subsurface fragments are omitted from non-Histosol soils because rock fragments in those O horizons can be incompletely populated and/or inaccurate. 

Median RV values of O horizon(s) depth and depth to restriction are provided here to be used for the depth range of subsurface fragments. The soil orders are also provided so that you can determine whether O horizons from Histosols have been included. 

Note: Below are proposed values to. Consider whether your data has sufficient detail to use a depth range beyond the depth class. 

Use the following value as depth top - this value is the median of the acreage weighted component RV depth of O horizon(s). Why is the component RV used instead of the combined component high/low as is done in the majority of this script? Organic horizon data comes from the component horizon table, and high/low values from the component horizon table are missing for many surveys.
```{r}

if (all(is.na(EDIT_updates_filtered$o_horizon_depth_in))) {print("No organic horizon depth populated in NASIS.")
  
} else{
  
median(o_horizon_stats, na.rm = TRUE)
  
}

```


Use the following value as depth bottom- this value is the median of the acreage weighted combined component high/low depth to restriction: 
```{r}
if (all(is.na(EDIT_updates_filtered$restrict_depth_h_in)) &
    all(is.na(EDIT_updates_filtered$restrict_depth_l_in))) {print("No restriction depth populated in NASIS. Use 60 inches as the bottom depth.")
  
} else{
  
median(restriction_r_sum_stats, na.rm = TRUE)
  
}

```


Subsurface frag (< 3 inches) prep code
```{r, message=FALSE}
# NAs in subsurface frags generally mean zero subsurface frags

EDIT_updates_filtered$subsurf_frags_hi_2mm.3in <- ifelse(is.na(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in), 0, EDIT_updates_filtered$subsurf_frags_hi_2mm.3in)

EDIT_updates_filtered$subsurf_frags_lo_2mm.3in <- ifelse(is.na(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in), 0, EDIT_updates_filtered$subsurf_frags_lo_2mm.3in)

if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in)) &
    all(is.na(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in))) {print("No data populated in NASIS for this field.")

} else{
  
# Define bin width and round by measures

my_vector <- c(1, 2, 5, 10, 20)

my_factor <- diff(range(c(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in, 
             EDIT_updates_filtered$subsurf_frags_lo_2mm.3in), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Prepare dataframe for component histogram

subsurf_frags.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, 
                                   subsurf_frags_hi_2mm.3in, subsurf_frags_lo_2mm.3in, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "frag",
    values_to = "value",
    cols = c("subsurf_frags_hi_2mm.3in", "subsurf_frags_lo_2mm.3in")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes
  
my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in, na.rm = TRUE) &
             min(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq   
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}
  
}
```


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

if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in))) {print("No data populated in NASIS for this field.")
  
} else{

  # Create component histogram

subsurf_frag_h1 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag == "subsurf_frags_hi_2mm.3in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag high (< 3 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

subsurf_frag_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

subsurf_frag_h2 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag == "subsurf_frags_hi_2mm.3in"), ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag high (< 3 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

subsurf_frag_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(subsurf_frag_h1,
                                 subsurf_frag_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Subsurface frag high (< 3 inches), acreage weighted 95th percentile** - use as 'high (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  subsurf_frag_stats <- rep(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                accuracy = 1, f = round)) |> sort(na.last = NA)
  
  subsurf_frag_stats |> quantile(probs = c(0.95), na.rm = TRUE, type = 1)

}

```

**Subsurface frag high (< 3 inches), acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in))) {print("No data populated in NASIS for this field.")

} else{

  summary(subsurf_frag_stats)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

subsurf_frags_l1 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag == "subsurf_frags_lo_2mm.3in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag low (< 3 inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

subsurf_frags_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

subsurf_frags_l2 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag == "subsurf_frags_lo_2mm.3in"), ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag low (< 3 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Subsurface frag (%)") +
    ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

subsurf_frags_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(subsurf_frags_l1,
                                 subsurf_frags_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Subsurface frag low (< 3 inches), acreage weighted 5th percentile** - use as 'low (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in))) {print("No data populated in NASIS for this field.")
  
} else{
  
    subsurf_frag_stats <- rep(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                          accuracy = 1, f = round)) |> sort(na.last = NA)
    
    subsurf_frag_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```


**Subsurface frag low (< 3 inches), acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in))) {print("No data populated in NASIS for this field.")

} else{

  summary(subsurf_frag_stats)

}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in)) &
    all(is.na(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

subsurf_frag_r1 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag %in% c("subsurf_frags_hi_2mm.3in", "subsurf_frags_lo_2mm.3in")) |> 
  dplyr::mutate(highlow = ifelse(frag == "subsurf_frags_hi_2mm.3in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag rep. (< 3 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

subsurf_frag_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

subsurf_frag_r2 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag %in% c("subsurf_frags_hi_2mm.3in", "subsurf_frags_lo_2mm.3in")) |> 
  dplyr::mutate(highlow = ifelse(frag == "subsurf_frags_hi_2mm.3in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag rep. (< 3 inches)") + 
  ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

subsurf_frag_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(subsurf_frag_r1,
                                 subsurf_frag_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Subsurface frag representative (< 3 inches) - combined component high/low 20th and 80th percentile** - use as 'low (representative)' & 'high (representative)' 
 
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in)) &
    all(is.na(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in))) {print("No data populated in NASIS for this field.")

} else{

  subsurf_frag_stats <- rep(subsurf_frags.df$value, plyr::round_any(subsurf_frags.df$acres/comp_acreage_unit, 
                            accuracy = 1, f = round)) |> sort(na.last = NA)
  
  subsurf_frag_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE, type = 1)

}

```

**Subsurface frag representative acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_2mm.3in)) &
    all(is.na(EDIT_updates_filtered$subsurf_frags_lo_2mm.3in))) {print("No data populated in NASIS for this field.")

} else{

  summary(subsurf_frag_stats)

}
```


#### **Subsurface frag (3+ inches)**

Subsurface frag (3+ inches) prep code
```{r, message=FALSE}

EDIT_updates_filtered$subsurf_frags_hi_gt_3in <- ifelse(is.na(EDIT_updates_filtered$subsurf_frags_hi_gt_3in), 0, EDIT_updates_filtered$subsurf_frags_hi_gt_3in)

EDIT_updates_filtered$subsurf_frags_lo_gt_3in <- ifelse(is.na(EDIT_updates_filtered$subsurf_frags_lo_gt_3in), 0, EDIT_updates_filtered$subsurf_frags_lo_gt_3in)

# Is variable NA?

if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_gt_3in)) &
    all(is.na(EDIT_updates_filtered$subsurf_frags_lo_gt_3in))) {print("No data populated in NASIS for this field.")
  
} else{
  
# Define bin width and round by measures

my_vector <- c(1, 2, 5, 10, 20)

my_factor <- diff(range(c(EDIT_updates_filtered$subsurf_frags_hi_gt_3in, 
             EDIT_updates_filtered$subsurf_frags_lo_gt_3in), na.rm = TRUE))/my_vector

reducer <- my_factor < 25

bin_width <- my_vector[reducer][which.min(abs(my_factor[reducer] - 15))]

round_by <- bin_width

# Default dataframe

subsurf_frags.df <-
  EDIT_updates_filtered |>  dplyr::select(coiid, subsurf_frags_lo_gt_3in,
                                   subsurf_frags_hi_gt_3in, comp_acres) |> 
  tidyr::pivot_longer(
    names_to = "frag",
    values_to = "value",
    cols = c("subsurf_frags_lo_gt_3in", "subsurf_frags_hi_gt_3in")) |> 
  dplyr::rename(acres = comp_acres)

# Sequence for axes
  
my_seq <- seq(plyr::round_any(min(EDIT_updates_filtered$subsurf_frags_lo_gt_3in, na.rm = TRUE), round_by, f = floor),
  plyr::round_any(max(EDIT_updates_filtered$subsurf_frags_hi_gt_3in, na.rm = TRUE), round_by, f = ceiling), bin_width)

my_seq <- if(min(my_seq) == min(EDIT_updates_filtered$subsurf_frags_lo_gt_3in, na.rm = TRUE) &
             min(EDIT_updates_filtered$subsurf_frags_lo_gt_3in, na.rm = TRUE) - bin_width >= 0){
  seq(plyr::round_any(min(EDIT_updates_filtered$subsurf_frags_lo_gt_3in, na.rm = TRUE), round_by, f = floor) - bin_width,
      plyr::round_any(max(EDIT_updates_filtered$subsurf_frags_hi_gt_3in, na.rm = TRUE), round_by, f = ceiling), bin_width)
} else{
  
  my_seq    
  
}

my_seq <- if(length(my_seq) == 1) {
  my_seq <- c(my_seq[1], my_seq[1] + bin_width)
} else {
  my_seq
}
  
}

```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_gt_3in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

subsurf_frag_h1 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag == "subsurf_frags_hi_gt_3in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag high (3 inches +)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

subsurf_frag_h1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

subsurf_frag_h2 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag == "subsurf_frags_hi_gt_3in"), ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag high (3 inches +)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

subsurf_frag_h2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(subsurf_frag_h1,
                                 subsurf_frag_h2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Subsurface frag high (3+ inches), 0 - 60 inches, acreage weighted 95th percentile** - use as 'high (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_gt_3in))) {print("No data populated in NASIS for this field.")
  
} else{
  
  subsurf_frag_stats <- rep(EDIT_updates_filtered$subsurf_frags_hi_gt_3in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                accuracy = 1, f = round)) |> sort(na.last = NA)
  
  subsurf_frag_stats |> quantile(probs = c(0.95), na.rm = TRUE, type = 1)

}

```

**Subsurface frag high (3+ inches), 0 - 60 inches, acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_gt_3in)) &
    all(is.na(EDIT_updates_filtered$subsurf_frags_lo_gt_3in))) {print("No data populated in NASIS for this field.")

} else{

  summary(subsurf_frag_stats)

}
```


```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$subsurf_frags_lo_gt_3in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

subsurf_frag_l1 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag == "subsurf_frags_lo_gt_3in"), ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag low (3+ inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

subsurf_frag_l1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

subsurf_frag_l2 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag == "subsurf_frags_lo_gt_3in"), ggplot2::aes(x = value)) + 
  ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag low (3+ inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

subsurf_frag_l2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(subsurf_frag_l1,
                                 subsurf_frag_l2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Subsurface frag low (3+ inches), 0 - 60 inches, acreage weighted 5th percentile** - use as 'low (actual)'
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$subsurf_frags_lo_gt_3in))) {print("No data populated in NASIS for this field.")
  
} else{
  
    subsurf_frag_stats <- rep(EDIT_updates_filtered$subsurf_frags_lo_gt_3in, plyr::round_any(EDIT_updates_filtered$comp_acres/comp_acreage_unit, 
                                          accuracy = 1, f = round)) |> sort(na.last = NA)
    
    subsurf_frag_stats |>  quantile(probs = c(0.05), na.rm = TRUE, type = 1)

}

```


**Subsurface frag low (3+ inches), 0 - 60 inches, acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$subsurf_frags_lo_gt_3in))) {print("No data populated in NASIS for this field.")

} else{

  summary(subsurf_frag_stats)

}
```

```{r message=FALSE, warning=FALSE}
# Is variable NA?

if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_gt_3in)) &
    all(is.na(EDIT_updates_filtered$subsurf_frags_lo_gt_3in))) {print("No data populated in NASIS for this field.")
  
} else{

# Create component histogram

subsurf_frag_r1 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag %in% c("subsurf_frags_hi_gt_3in", "subsurf_frags_lo_gt_3in")) |> 
  dplyr::mutate(highlow = ifelse(frag == "subsurf_frags_hi_gt_3in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2)), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag rep. (3+ inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert histogram to plotly

subsurf_frag_r1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage weighted histogram

subsurf_frag_r2 <- ggplot2::ggplot(subsurf_frags.df |>  dplyr::filter(frag %in% c("subsurf_frags_hi_gt_3in", "subsurf_frags_lo_gt_3in")) |> 
  dplyr::mutate(highlow = ifelse(frag == "subsurf_frags_hi_gt_3in", "high", "low")) |> 
  dplyr::mutate(coiid2 = paste0(coiid, " (", highlow, ")")), 
                ggplot2::aes(x = value)) + ggplot2::geom_histogram(ggplot2::aes(text = paste("coiid(s):", coiid2), weight = acres), 
  binwidth = bin_width,
  center = min(my_seq) + bin_width/2,
  color = "black",
  fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Subsurface frag rep. (3+ inches)") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Subsurface frag (%)") +
  ggplot2::scale_x_continuous(limits = c(min(my_seq), max(my_seq)),
  breaks = my_seq[seq(1, length(my_seq), by = 2)])

# Convert to plotly

subsurf_frag_r2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures

manipulateWidget::combineWidgets(subsurf_frag_r1,
                                 subsurf_frag_r2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))

}

```

**Subsurface frag representative (3+ inches), 0 - 60 inches, acreage weighted 20th and 80th percentile** -  use as 'low (representative)' and 'high (representative)' 
```{r, fig.show="hold", out.width="50%"}
if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_gt_3in)) &
    all(is.na(EDIT_updates_filtered$subsurf_frags_lo_gt_3in))) {print("No data populated in NASIS for this field.")

} else{

  subsurf_frag_stats <- rep(subsurf_frags.df$value, plyr::round_any(subsurf_frags.df$acres/comp_acreage_unit, 
                            accuracy = 1, f = round)) |> sort(na.last = NA)
  
  subsurf_frag_stats |>  quantile(c(0.2, 0.8), na.rm = TRUE, type = 1)

}

```

**Subsurface frag representative (3+ inches), 0 - 60 inches, acreage weighted summary statistics**
```{r}

if (all(is.na(EDIT_updates_filtered$subsurf_frags_hi_gt_3in)) &
    all(is.na(EDIT_updates_filtered$subsurf_frags_lo_gt_3in))) {print("No data populated in NASIS for this field.")

} else{

  summary(subsurf_frag_stats)

}
```

#### **Component name**

If there are 15 or fewer components, a histogram and table will be shown. Else, only a table will be shown. 
```{r, message=FALSE, warning=FALSE}

if (all(is.na(EDIT_updates_filtered$compname))) {print("No data populated in NASIS for this field.")
  
} else{
  
  if (length(unique(EDIT_updates_filtered$compname)) <= 15) {
    
    comp.df <- EDIT_updates_filtered |>  dplyr::select(coiid, compname, comp_acres)  
    
    #|>  dplyr::group_by(compname) |>  dplyr::summarise(count = dplyr::n(), acres = sum(comp_acres), coiid = toString(coiid)) |>  dplyr::ungroup() 

    
# Component histogram 
comp1 <- ggplot2::ggplot(comp.df, ggplot2::aes(compname)) + 
  ggplot2::geom_bar(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Component name") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 45)) +
  ggplot2::ylab("# of components") + ggplot2::xlab("Component name")

# Convert to plotly
comp1 <- plotly::ggplotly(tooltip = c("text", "value"))

# Acreage histogram
comp2 <- ggplot2::ggplot(comp.df, ggplot2::aes(compname, comp_acres)) +
  ggplot2::geom_col(ggplot2::aes(text = paste("coiid(s):", coiid)),
                    color = "black", fill = "grey", closed = "left") + 
  ggplot2::ggtitle("Component name") + ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5),
                                                      axis.text.x = ggplot2::element_text(angle = 45)) +
  ggplot2::ylab("# of acres") + ggplot2::xlab("Component name")

#Convert to plotly
comp2 <- plotly::ggplotly(tooltip = c("text", "value"))

# Join both figures
manipulateWidget::combineWidgets(comp1,
                                 comp2,
                                 nrow = 1,
                                 ncol = 2,
                                 colsize = c(1, 1))


  } else {
    
    EDIT_updates_filtered |> dplyr::select(compname, comp_acres, coiid) |> dplyr::group_by(compname) |> 
    dplyr::summarise(acres = sum(comp_acres), `# components` = dplyr::n(), coiids = toString(coiid)) |> 
    dplyr::arrange(dplyr::desc(acres)) |> knitr::kable()
    
  }

}

```