diff --git a/DESCRIPTION b/DESCRIPTION
index 60c2107..5ff584c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -10,6 +10,10 @@ Authors@R: c(person(given = "Marc",
                     family = "Hill",
                     role = "ctb",
                     email = "hill.ryan@epa.gov"),
+             person(given = "Selia",
+                    family = "Markley",
+                    role = "ctb",
+                    email = "markley.selia@epa.gov"),
              person(given = "Travis",
                     family = "Hudson",
                     role = "ctb",
@@ -39,10 +43,14 @@ Imports:
     nhdplusTools,
     jsonlite,
     httr2,
-    curl (>= 6.0.0)
+    curl (>= 6.0.0),
+    ggpattern,
+    patchwork,
+    cowplot,
+    tigris,
+    ggplot2
 Suggests: 
   dplyr,
-  ggplot2,
   mapview,
   testthat, 
   knitr, 
diff --git a/NAMESPACE b/NAMESPACE
index 4ac1597..bf661fb 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,11 +4,16 @@ export(lc_fullname)
 export(lc_get_comid)
 export(lc_get_data)
 export(lc_get_metric_names)
+export(lc_get_nni)
 export(lc_get_params)
 export(lc_nlcd)
+export(lc_plotnni)
 export(sc_fullname)
 export(sc_get_comid)
 export(sc_get_data)
 export(sc_get_metric_names)
 export(sc_get_params)
 export(sc_nlcd)
+export(sc_plotnni)
+import(ggpattern)
+import(ggplot2)
diff --git a/R/lc_get_data.R b/R/lc_get_data.R
index 4d3200f..94947b7 100644
--- a/R/lc_get_data.R
+++ b/R/lc_get_data.R
@@ -171,6 +171,8 @@ lc_get_data <- function(comid = NULL,
 #' @title Get NLCD Data
 #'
 #' @description
+#' `r lifecycle::badge("deprecated")`
+#' `lc_nlcd()` was renamed to `lc_get_nlcd()` to create a more consistent API.
 #' Function to specifically retrieve all NLCD metrics for a given year using the StreamCat API.
 #'
 #' @author
@@ -205,21 +207,23 @@ lc_get_data <- function(comid = NULL,
 #' @examples
 #' \dontrun{
 #'
-#' df <- lc_nlcd(comid='23783629', year='2019', aoi='ws')
-#'
-#' df <- lc_nlcd(year='2016', aoi='cat',
+#' df <- lc_nlcd(comid='23783629', year='2019', aoi='ws')  # Will show a deprecation warning
+#' 
+#' df <- lc_get_nlcd(comid='23783629', year='2019', aoi='ws')
+#' 
+#' df <- lc_get_nlcd(year='2016', aoi='cat',
 #' comid='23783629,23794487,23812618', showAreaSqKm=FALSE, showPctFull=TRUE)
 #'
-#' df <- lc_nlcd(year='2016', aoi='cat',
+#' df <- lc_get_nlcd(year='2016', aoi='cat',
 #' comid='23783629,23794487,23812618', countOnly=TRUE)
 #'
-#' df <- lc_nlcd(year='2016, 2019', aoi='cat,ws',
+#' df <- lc_get_nlcd(year='2016, 2019', aoi='cat,ws',
 #' comid='23783629,23794487,23812618')
 #' }
 #' @export
 
 
-lc_nlcd <- function(year = '2019', aoi = NULL, comid = NULL,
+lc_get_nlcd <- function(year = '2019', aoi = NULL, comid = NULL,
                     showAreaSqKm = NULL, showPctFull = NULL,
                     countOnly = NULL) {
   # year must be a character string.
@@ -282,3 +286,184 @@ lc_nlcd <- function(year = '2019', aoi = NULL, comid = NULL,
   # End of function. Return a data frame.
   return(final_df)
 }
+
+#' @rdname lc_get_nlcd
+#' @export
+#' @keywords internal
+lc_nlcd <- function(year = '2019',
+                    comid = NULL,
+                    aoi = NULL,
+                    showAreaSqKm = NULL,
+                    showPctFull = NULL,
+                    state = NULL,
+                    county = NULL,
+                    region = NULL,
+                    conus = NULL,
+                    countOnly = NULL) {
+  lifecycle::deprecate_warn("0.10.0", "lc_nlcd()", "lc_get_nlcd()")
+  lc_get_nlcd(year = '2019',
+              comid = NULL,
+              aoi = NULL,
+              showAreaSqKm = NULL,
+              showPctFull = NULL,
+              state = NULL,
+              county = NULL,
+              region = NULL,
+              conus = NULL,
+              countOnly = NULL)
+}
+
+#' @title Get NNI
+#' 
+#' @description 
+#' Function to get all NNI data available for a given year.
+#' 
+#' @author
+#' Selia Markley
+#'
+#' @param year Years(s) of NNI metrics to query.
+#' Only valid NNI years are accepted (1987:2017)
+#' Syntax: year=<year1>,<year2>
+#' 
+#' @param aoi Specify the area of interest described by a metric. By default, all available areas of interest
+#' for a given metric are returned.
+#' Syntax: areaOfInterest=<value1>,<value2>
+#' Values: catchment|watershed
+#' 
+#' @param comid Return metric information for specific COMIDs
+#' Syntax: comid=<comid1>,<comid2>
+#' 
+#' @param showAreaSqKm Return the area in square kilometers of a given area of interest.
+#' The default value is true.
+#' Values: true|false
+#' 
+#' @param showPctFull Return the pctfull for each dataset. The default value is false.
+#' Values: true|false
+#' 
+#' @param countOnly Return a CSV containing only the row count (ROWCOUNT) and the column
+#' count (COLUMNCOUNT) that the server expects to return in a request. The default value is false.
+#' Values: true|false
+#'
+#' @return A tibble of desired StreamCat metrics
+#' @export
+#'
+#' @examples\donttest{
+#' df <- lc_get_nni(year='1987, 1990, 2005, 2017', aoi='cat,ws',
+#' comid='23783629,23794487,23812618')
+#' 
+#' df <- lc_get_nni(year='2015', aoi='cat',
+#' comid='23783629', countOnly=TRUE)
+#' 
+#' df <- lc_get_nni(comid='23783629', year='2011, 2012', aoi='ws')
+#' }
+
+lc_get_nni <- function(year, aoi = NULL, comid = NULL,
+                    showAreaSqKm = TRUE, showPctFull = NULL,
+                    countOnly = NULL) {
+  # year must be a character string.
+  year_chr <-  as.character(year)
+  # split multiple years supplied as a single string into
+  # a vector of years.
+  year_vec <- unlist(strsplit(x = year_chr,
+                              split = ",|, "))
+  # Vector of valid NNI years to check inputs against.
+  valid_years <- c('1987',
+                   '1988',
+                   '1989',
+                   '1990',
+                   '1991',
+                   '1992',
+                   '1993',
+                   '1994',
+                   '1995',
+                   '1996',
+                   '1997',
+                   '1998',
+                   '1999',
+                   '2000',
+                   '2001',
+                   '2002',
+                   '2003',
+                   '2004',
+                   '2005',
+                   '2006',
+                   '2007',
+                   '2008',
+                   '2009',
+                   '2010',
+                   '2011',
+                   '2012',
+                   '2013',
+                   '2014',
+                   '2015',
+                   '2016',
+                   '2017')
+  # Stop early if any of the year(s) supplied are not found in the valid
+  # years vec.
+  stopifnot(
+    "year must be a valid NNI year" = any(year_vec %in% valid_years)
+  )
+  # Vector of NNI metric names.
+  nni <- c(
+      'n_leg_',
+      'n_ags_',
+      'n_ff_',
+      'n_uf_',
+      'n_cf_',
+      'n_cr_',
+      'n_hw_',
+      'n_lw_',
+      'p_leg_',
+      'p_ags_',
+      'p_ff_',
+      'p_uf_',
+      'p_cr_',
+      'p_hw_',
+      'p_lw_'
+    )
+  # Add n_dep for available years
+  ndep_year_vec <- year_vec[!year_vec %in% c('1987', '1988', '1989')]
+  ndep_comb <- expand.grid('n_dep_', ndep_year_vec)
+  ndep_mets <- paste0(ndep_comb$Var1,
+                     ndep_comb$Var2,
+                     collapse = ",",
+                     recycle0 = TRUE)
+  # Add p_dep for available years
+  pdep_year_vec <- year_vec[!year_vec %in% c('1987', '1988', '1989', '1990', '1991', '1992', '1993', '1994', '1995', '1996',
+                                             '2014', '2015', '2016', '2017')]
+  pdep_comb <- expand.grid('p_dep_', pdep_year_vec)
+  pdep_mets <- paste0(pdep_comb$Var1,
+                      pdep_comb$Var2,
+                      collapse = ",",
+                      recycle0 = TRUE)
+  # Add n_usgsww and p_usgsww for available years
+  ww_year_vec <- year_vec[year_vec %in% c('1988', '1990', '1992', '1996', '2000', '2004', '2008', '2012')]
+  ww_comb <- expand.grid(c('p_usgsww_', 'n_usgsww_'), ww_year_vec)
+  ww_mets <- paste0(ww_comb$Var1,
+                      ww_comb$Var2,
+                      collapse = ",",
+                      recycle0 = TRUE)
+  # Create a data frame of all NNI Metric and year combinations.
+  all_comb <- expand.grid(nni, year_vec)
+  # Concatenate the NLCD metric name with the supplied year(s) to create
+  # valid metric names to submit to the API.
+  nni_mets <- paste0(all_comb$Var1,
+                      all_comb$Var2,
+                      collapse = ",",
+                      recycle0 = TRUE)
+  # Combine all NNI metrics
+  nni_mets_all <- paste0(nni_mets, ",", ndep_mets, ",", pdep_mets, ",", ww_mets)
+  
+  # Query the API.
+  final_df <- lc_get_data(
+    metric = nni_mets_all,
+    aoi = aoi,
+    comid = comid,
+    showAreaSqKm = showAreaSqKm,
+    showPctFull = showPctFull,
+    countOnly = countOnly
+  )
+  # End of function. Return a data frame.
+  return(final_df)
+}
+
diff --git a/R/lc_plot.R b/R/lc_plot.R
new file mode 100644
index 0000000..6ca1c33
--- /dev/null
+++ b/R/lc_plot.R
@@ -0,0 +1,326 @@
+#' Plot National Nutrient Inventory data for lakes
+#'
+#' @description
+#' Function to plot time series of nitrogen and phosphorus budgets for a given lake
+#' COMID. This function allows a user to return a time series of major inputs, 
+#' outputs, and derived metrics of nitrogen and phosphorus. Plot is returned as an 
+#' object
+#' 
+#' @author 
+#' Selia Markley
+#'
+#' @param comid Identifier of lake COMID user wants to plot NNI data for. Must be a character string
+#' with the COMID digit.
+#' Syntax: com=<COMID>
+#' 
+#' @param include.nue Include time series of nitrogen use efficiency in the returned plot. 
+#' The default value is false. 
+#' Values: true|false 
+#'
+#' @return
+#' Return plot as an object.
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' p <- lc_plotnni(comid='23794487')
+#' p <- lc_plotnni(comid='23794487', include.nue=TRUE)
+#' }
+
+lc_plotnni <- function(comid, include.nue = FALSE){
+  
+  # Get StreamCat data
+  nni <- lc_get_data(metric = 'n_dep_1990,n_ff_1990,n_uf_1990,n_lw_1990,n_hw_1990,n_ags_1990,n_cf_1990,n_cr_1990,p_cr_1990,p_lw_1990,p_hw_1990,p_uf_1990,p_ff_1990,p_ags_1990,n_dep_1991,n_ff_1991,n_uf_1991,n_lw_1991,n_hw_1991,n_ags_1991,n_cf_1991,n_cr_1991,p_cr_1991,p_lw_1991,p_hw_1991,p_uf_1991,p_ff_1991,p_ags_1991,n_dep_1992,n_ff_1992,n_uf_1992,n_lw_1992,n_hw_1992,n_ags_1992,n_cf_1992,n_cr_1992,p_cr_1992,p_lw_1992,p_hw_1992,p_uf_1992,p_ff_1992,p_ags_1992,n_dep_1993,n_ff_1993,n_uf_1993,n_lw_1993,n_hw_1993,n_ags_1993,n_cf_1993,n_cr_1993,p_cr_1993,p_lw_1993,p_hw_1993,p_uf_1993,p_ff_1993,p_ags_1993,n_dep_1994,n_ff_1994,n_uf_1994,n_lw_1994,n_hw_1994,n_ags_1994,n_cf_1994,n_cr_1994,p_cr_1994,p_lw_1994,p_hw_1994,p_uf_1994,p_ff_1994,p_ags_1994,n_dep_1995,n_ff_1995,n_uf_1995,n_lw_1995,n_hw_1995,n_ags_1995,n_cf_1995,n_cr_1995,p_cr_1995,p_lw_1995,p_hw_1995,p_uf_1995,p_ff_1995,p_ags_1995,n_dep_1996,n_ff_1996,n_uf_1996,n_lw_1996,n_hw_1996,n_ags_1996,n_cf_1996,n_cr_1996,p_cr_1996,p_lw_1996,p_hw_1996,p_uf_1996,p_ff_1996,p_ags_1996,n_dep_1997,n_ff_1997,n_uf_1997,n_lw_1997,n_hw_1997,n_ags_1997,n_cf_1997,n_cr_1997,p_cr_1997,p_lw_1997,p_hw_1997,p_uf_1997,p_ff_1997,p_ags_1997,n_dep_1998,n_ff_1998,n_uf_1998,n_lw_1998,n_hw_1998,n_ags_1998,n_cf_1998,n_cr_1998,p_cr_1998,p_lw_1998,p_hw_1998,p_uf_1998,p_ff_1998,p_ags_1998,n_dep_1999,n_ff_1999,n_uf_1999,n_lw_1999,n_hw_1999,n_ags_1999,n_cf_1999,n_cr_1999,p_cr_1999,p_lw_1999,p_hw_1999,p_uf_1999,p_ff_1999,p_ags_1999,n_dep_2000,n_ff_2000,n_uf_2000,n_lw_2000,n_hw_2000,n_ags_2000,n_cf_2000,n_cr_2000,p_cr_2000,p_lw_2000,p_hw_2000,p_uf_2000,p_ff_2000,p_ags_2000,n_dep_2001,n_ff_2001,n_uf_2001,n_lw_2001,n_hw_2001,n_ags_2001,n_cf_2001,n_cr_2001,p_cr_2001,p_lw_2001,p_hw_2001,p_uf_2001,p_ff_2001,p_ags_2001,n_dep_2002,n_ff_2002,n_uf_2002,n_lw_2002,n_hw_2002,n_ags_2002,n_cf_2002,n_cr_2002,p_cr_2002,p_lw_2002,p_hw_2002,p_uf_2002,p_ff_2002,p_ags_2002,n_dep_2003,n_ff_2003,n_uf_2003,n_lw_2003,n_hw_2003,n_ags_2003,n_cf_2003,n_cr_2003,p_cr_2003,p_lw_2003,p_hw_2003,p_uf_2003,p_ff_2003,p_ags_2003,n_dep_2004,n_ff_2004,n_uf_2004,n_lw_2004,n_hw_2004,n_ags_2004,n_cf_2004,n_cr_2004,p_cr_2004,p_lw_2004,p_hw_2004,p_uf_2004,p_ff_2004,p_ags_2004,n_dep_2005,n_ff_2005,n_uf_2005,n_lw_2005,n_hw_2005,n_ags_2005,n_cf_2005,n_cr_2005,p_cr_2005,p_lw_2005,p_hw_2005,p_uf_2005,p_ff_2005,p_ags_2005,n_dep_2006,n_ff_2006,n_uf_2006,n_lw_2006,n_hw_2006,n_ags_2006,n_cf_2006,n_cr_2006,p_cr_2006,p_lw_2006,p_hw_2006,p_uf_2006,p_ff_2006,p_ags_2006,n_dep_2007,n_ff_2007,n_uf_2007,n_lw_2007,n_hw_2007,n_ags_2007,n_cf_2007,n_cr_2007,p_cr_2007,p_lw_2007,p_hw_2007,p_uf_2007,p_ff_2007,p_ags_2007,n_dep_2008,n_ff_2008,n_uf_2008,n_lw_2008,n_hw_2008,n_ags_2008,n_cf_2008,n_cr_2008,p_cr_2008,p_lw_2008,p_hw_2008,p_uf_2008,p_ff_2008,p_ags_2008,n_dep_2009,n_ff_2009,n_uf_2009,n_lw_2009,n_hw_2009,n_ags_2009,n_cf_2009,n_cr_2009,p_cr_2009,p_lw_2009,p_hw_2009,p_uf_2009,p_ff_2009,p_ags_2009,n_dep_2010,n_ff_2010,n_uf_2010,n_lw_2010,n_hw_2010,n_ags_2010,n_cf_2010,n_cr_2010,p_cr_2010,p_lw_2010,p_hw_2010,p_uf_2010,p_ff_2010,p_ags_2010,n_dep_2011,n_ff_2011,n_uf_2011,n_lw_2011,n_hw_2011,n_ags_2011,n_cf_2011,n_cr_2011,p_cr_2011,p_lw_2011,p_hw_2011,p_uf_2011,p_ff_2011,p_ags_2011,n_dep_2012,n_ff_2012,n_uf_2012,n_lw_2012,n_hw_2012,n_ags_2012,n_cf_2012,n_cr_2012,p_cr_2012,p_lw_2012,p_hw_2012,p_uf_2012,p_ff_2012,p_ags_2012,n_dep_2013,n_ff_2013,n_uf_2013,n_lw_2013,n_hw_2013,n_ags_2013,n_cf_2013,n_cr_2013,p_cr_2013,p_lw_2013,p_hw_2013,p_uf_2013,p_ff_2013,p_ags_2013,n_dep_2014,n_ff_2014,n_uf_2014,n_lw_2014,n_hw_2014,n_ags_2014,n_cf_2014,n_cr_2014,p_cr_2014,p_lw_2014,p_hw_2014,p_uf_2014,p_ff_2014,p_ags_2014,n_dep_2015,n_ff_2015,n_uf_2015,n_lw_2015,n_hw_2015,n_ags_2015,n_cf_2015,n_cr_2015,p_cr_2015,p_lw_2015,p_hw_2015,p_uf_2015,p_ff_2015,p_ags_2015,n_dep_2016,n_ff_2016,n_uf_2016,n_lw_2016,n_hw_2016,n_ags_2016,n_cf_2016,n_cr_2016,p_cr_2016,p_lw_2016,p_hw_2016,p_uf_2016,p_ff_2016,p_ags_2016,n_dep_2017,n_ff_2017,n_uf_2017,n_lw_2017,n_hw_2017,n_ags_2017,n_cf_2017,n_cr_2017,p_cr_2017,p_lw_2017,p_hw_2017,p_uf_2017,p_ff_2017,p_ags_2017',
+                     aoi='ws',
+                     comid = comid,
+                     showAreaSqKm = FALSE,
+                     showPctFull = FALSE)
+  
+  #Create N inputs df
+  
+  nin <- nni[, grepl("^(n)", names(nni)) & !grepl("(cr)", names(nni)) & !grepl("(ags)", names(nni))]
+  
+  names(nin) <- sapply(names(nin), function(col){
+    substr(col, 3, nchar(col) -2)
+  })
+  
+  nin <- nin |>
+    tidyr::pivot_longer(
+      cols = everything(),
+      names_to = c("metric", "year"),
+      names_sep = "_",
+      values_to = "value"
+    ) |>
+    dplyr::mutate(year = as.integer(year)) |>
+    dplyr::mutate(value = value / 1000000)
+  
+  #Create P inputs df
+  
+  pin <- nni[, grepl("^(p)", names(nni)) & !grepl("(cr)", names(nni)) & !grepl("(ags)", names(nni))]
+  
+  names(pin) <- sapply(names(pin), function(col){
+    substr(col, 3, nchar(col) -2)
+  })
+  
+  pin <- pin |>
+    tidyr::pivot_longer(
+      cols = everything(),
+      names_to = c("metric", "year"),
+      names_sep = "_",
+      values_to = "value"
+    ) |>
+    dplyr::mutate(year = as.integer(year)) |>
+    dplyr::mutate(value = value / 1000000)
+  
+  #Create N dfs for lines (cr, agsur, nue)
+  
+  nlines <- nni[, grepl("^n_cr|^n_ags", names(nni))]
+  
+  names(nlines) <- sapply(names(nlines), function(col){
+    substr(col, 3, nchar(col) -2)
+  })
+  
+  nlines <- nlines |>
+    tidyr::pivot_longer(
+      cols = everything(),
+      names_to = c("metric","year"),
+      names_sep = "_",
+      values_to = "value"
+    ) |>
+    dplyr::mutate(year = as.integer(year)) |>
+    dplyr::mutate(value = value / 1000000) |>
+    tidyr::pivot_wider(
+      names_from = 'metric', 
+      values_from = 'value'
+    ) |>
+    dplyr::mutate(totag = ags + cr) |>
+    dplyr::mutate(nue = (cr / totag) * 100) |>
+    tidyr::pivot_longer(
+      cols = !year, 
+      names_to="metric", 
+      values_to="value")
+  
+  ncrag <- nlines |>
+    dplyr::filter(metric %in% c('ags', 'cr'))
+  
+  nue <- nlines |>
+    dplyr::filter(metric == 'nue') |>
+    tidyr::pivot_wider(names_from = 'metric',
+                values_from = 'value')
+  
+  #Create P dfs for lines (cr, agsur, pue)
+  
+  plines <- nni[, grepl("^p_cr|^p_ags", names(nni))]
+  
+  names(plines) <- sapply(names(plines), function(col){
+    substr(col, 3, nchar(col) -2)
+  })
+  
+  
+  plines <- plines |>
+    tidyr::pivot_longer(
+      cols = everything(),
+      names_to = c("metric","year"),
+      names_sep = "_",
+      values_to = "value"
+    ) |>
+    dplyr::mutate(year = as.integer(year)) |>
+    dplyr::mutate(value = value / 1000000) |>
+    tidyr::pivot_wider(
+      names_from = 'metric', 
+      values_from = 'value'
+    ) |>
+    dplyr::mutate(totag = ags + cr) |>
+    dplyr::mutate(nue = (cr / totag) * 100) |>
+    tidyr::pivot_longer(
+      cols = !year, 
+      names_to="metric", 
+      values_to="value")
+  
+  pcrag <- plines |>
+    dplyr::filter(metric %in% c('ags', 'cr'))
+  
+  pue <- plines |>
+    dplyr::filter(metric == 'nue') |>
+    tidyr::pivot_wider(names_from = 'metric',
+                values_from = 'value')
+  
+  pdf <- dplyr::bind_rows(plines, pin)
+  
+  ndf <- dplyr::bind_rows(nlines, nin)
+  
+  #create estimate column
+  knownfertyrs <- c(1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,
+                    2005,2006,2007,2008,2009,2010,2011,2012,2017)
+  nwsin <- nin |>
+    dplyr::mutate(estimated=dplyr::case_when(
+      metric == "dep" ~ FALSE,
+      metric == "hw" ~ FALSE,
+      metric == "cf" & year %in% c(1987,1992,1997,2002,2007,2012, 2017) ~ FALSE,
+      metric == "ff" & year %in% knownfertyrs ~ FALSE,
+      metric == "uf" & year %in% knownfertyrs ~ FALSE,
+      metric == "lw" & year %in% c(1987,1992,1997,2002,2007,2012,2017) ~ FALSE,
+      TRUE ~ TRUE
+    )) 
+  
+  pwsin <- pin |>
+    dplyr::filter(metric != 'cr') |>
+    dplyr::mutate(estimated=dplyr::case_when(
+      metric == "hw" ~ FALSE,
+      metric == "lw" & year %in% c(1987,1992,1997,2002,2007,2012,2017) ~ FALSE,
+      metric == "ff" & year %in% knownfertyrs ~ FALSE,
+      metric == "uf" & year %in% knownfertyrs ~ FALSE,
+      TRUE ~ TRUE
+    ))
+  
+  #get ready for plot
+  colorsn <- c('ff' = '#A3CC51', 'lw'='#B26F2C','hw'='#E51932','uf'='black', 'dep'='#6db6ff', 'cf'='#FFD700')
+  colorsp <- c('ff' = '#A3CC51', 'lw'='#B26F2C','hw'='#E51932','uf'='black')
+  
+  nwsin$metric <- factor(nwsin$metric, levels = c('uf','hw','dep','lw','cf','ff'))
+  pwsin$metric <- factor(pwsin$metric, levels = c('uf','hw','lw','ff'))
+  
+  
+  #create titles with higher level include.nue param
+  if (include.nue == TRUE){
+    nbartitle <- 'b)'
+    pbartitle <- 'd)'
+    nuetitle <- 'a)'
+    puetitle <- 'c)'
+  } else{
+    nbartitle <- 'a)'
+    pbartitle <- 'b)'
+    nuetitle <- ' '
+    puetitle <- ' '
+  }
+  
+  #create N bar plot
+  nbar <- ggplot() + 
+    ggpattern::geom_bar_pattern(data = nwsin, 
+                     aes(x=year,y=value, fill=metric, 
+                         pattern=factor(estimated, levels=c(TRUE,FALSE),
+                                        labels=c('Estimated','Non-Estimated'))),
+                     pattern = ifelse(nwsin$estimated, 'stripe','none'),
+                     pattern_color='white',
+                     pattern_density=0.05, 
+                     pattern_fill = 'white', 
+                     pattern_alpha = 0.5, 
+                     pattern_spacing=0.025,
+                     stat='identity', position='stack', 
+                     pattern_size=0.05) +
+    labs(title = 'Nitrogen (million kg)',
+         y = "Budget",
+         x = " ") +
+    scale_fill_manual(values=colorsn,
+                      labels = c('ff' = 'Farm Fertilizer',
+                                 'uf' = 'Urban Fertilizer',
+                                 'cf' = 'Crop N-Fixation',
+                                 'lw' = 'Livestock Manure',
+                                 'hw' = 'Human Waste',
+                                 'dep' = 'Total Deposition')) + 
+    scale_pattern_manual(name='Estimate Status',
+                         values=c('Estimated'='stripe','Non-estimated'='none')) + 
+    geom_line(data=ncrag, 
+              aes(x=year,y=value, linetype=metric), 
+              linewidth=1.25, color="black") + 
+    scale_linetype_manual(values = 
+                            c("ags"="solid", "cr"="dotted"),
+                          labels = c('ags' = 'Agricultural Surplus',
+                                     'cr' = 'Crop Removal')) +
+    guides(fill=
+             guide_legend(order=1, override.aes = list(pattern='none')), 
+           pattern=
+             guide_legend(order=2, override.aes = list(fill='grey')), 
+           linetype=
+             guide_legend(title=NULL)) + 
+    scale_x_continuous(breaks=seq(1987,2017,by=5)) + 
+    scale_color_manual(values=c("Agricultural Surplus"="black"), 
+                       guide = 'none') + 
+    theme_bw() + 
+    theme(plot.title = element_text(size=9, face="bold"), 
+          axis.title.y = element_text(size=9), 
+          legend.background = element_rect(fill="white", colour = "black"), 
+          legend.title = element_blank())
+  
+  #create p bar plot
+  pbar <- ggplot() + 
+    ggpattern::geom_bar_pattern(data = pwsin, 
+                     aes(x=year,y=value, fill=metric, 
+                         pattern=factor(estimated, levels=c(TRUE,FALSE),
+                                        labels=c('Estimated','Non-Estimated'))), 
+                     pattern = ifelse(pwsin$estimated, 'stripe','none'),
+                     pattern_color='white',
+                     pattern_density=0.05, 
+                     pattern_fill = 'white', 
+                     pattern_alpha = 0.5, 
+                     pattern_spacing=0.025,
+                     stat='identity', position='stack', pattern_size=0.05) +
+    labs(title = 'Phosphorus (million kg)',
+         y = "Budget",
+         x = " ") +
+    scale_fill_manual(values=colorsp) + 
+    scale_pattern_manual(name='Estimate Status',
+                         values=c('Estimated'='stripe',
+                                  'Non-estimated'='none')) + 
+    geom_line(data=pcrag, 
+              aes(x=year,y=value, 
+                  linetype=metric), 
+              linewidth=1.25, color="black") + 
+    scale_linetype_manual(values = 
+                            c("ags"="solid", "cr"="dotted")) +
+    guides(fill=
+             guide_legend(order=1, override.aes = list(pattern='none')), 
+           pattern=
+             guide_legend(order=2, override.aes = list(fill='grey')), 
+           linetype=
+             guide_legend(title=NULL)) + 
+    guides(fill="none", pattern = "none", linetype="none") +
+    scale_x_continuous(breaks=seq(1987,2017,by=5)) + 
+    scale_color_manual(values=c("Agricultural Surplus"="black"), 
+                       guide = 'none') + 
+    theme_bw() + 
+    theme(plot.title = element_text(size=9, face="bold"), 
+          axis.title.y = element_text(size=9), 
+          legend.background = element_rect(fill="white", colour = "black"), 
+          legend.title = element_blank())
+  
+  #create nue line plots
+  nue <- ggplot() + 
+    geom_line(data=nue, aes(x=year,y=nue), linewidth=1.25, color='seagreen')+
+    theme_bw() + 
+    scale_x_continuous(breaks=seq(1987,2017,by=5)) + 
+    labs(title = 'Nitrogen Use Efficiency',
+         y = "%", 
+         x=" ") + 
+    theme(plot.title = element_text(size=9, face="bold"), 
+          axis.title.x = element_text(size=9), 
+          axis.title.y = element_text(size=9)) 
+  
+  pue <- ggplot() + 
+    geom_line(data=pue, aes(x=year,y=nue, lty='Nutrient Use Efficiency'), linewidth=1.25, color="seagreen") + 
+    theme_bw() + 
+    scale_x_continuous(breaks=seq(1987,2017,by=5)) + 
+    labs(title = 'Phosphorus Use Efficiency',
+         y = "%", 
+         x="Year") + 
+    theme(plot.title = element_text(size=9, face="bold"), 
+          axis.title.x = element_text(size=9), 
+          axis.title.y = element_text(size=9, hjust=0.5),
+          legend.background = element_rect(fill="white", colour = "black"), 
+          legend.title = element_blank())  +
+    guides(fill="none", pattern = "none", linetype="none")
+  
+  #export final figure
+  inputs <- patchwork::wrap_plots(nbar, pbar, ncol=1, guides="collect")
+  nue <- patchwork::wrap_plots(nue, pue, ncol=1, guides="collect")
+  
+  if (include.nue == TRUE){
+    timenni <- patchwork::wrap_plots(nue, inputs, ncol=2) 
+  }
+  else {
+    timenni <- inputs 
+  }
+  
+  return(timenni)
+  
+}
diff --git a/R/sc_get_data.R b/R/sc_get_data.R
index 906be8f..3edea51 100644
--- a/R/sc_get_data.R
+++ b/R/sc_get_data.R
@@ -188,8 +188,11 @@ sc_get_data <- function(comid = NULL,
 #' @title Get NLCD Data
 #'
 #' @description
+#' `r lifecycle::badge("deprecated")`
+#' `sc_nlcd()` was renamed to `sc_get_nlcd()` to create a more consistent API.
 #' Function to retrieve all NLCD metrics for a given year using the StreamCat API.
 #'
+#' 
 #' @author
 #' Marc Weber
 #'
@@ -236,22 +239,25 @@ sc_get_data <- function(comid = NULL,
 #'
 #' @examples
 #' \dontrun{
-#' df <- sc_nlcd(year='2001', aoi='cat',comid='179,1337,1337420')
+#' 
+#' df <- sc_nlcd(year='2001', aoi='cat',comid='179') # Will show a deprecation warning
+#' 
+#' df <- sc_get_nlcd(year='2001', aoi='cat',comid='179,1337,1337420')
 #'
-#' df <- sc_nlcd(year='2001', aoi='ws', region='Region01')
+#' df <- sc_get_nlcd(year='2001', aoi='ws', region='Region01')
 #'
-#' df <- sc_nlcd(year='2001', aoi='ws', region='Region01',
+#' df <- sc_get_nlcd(year='2001', aoi='ws', region='Region01',
 #' countOnly=TRUE)
 #'
 #' df <- sc_nlcd(year='2001', aoi='ws', region='Region01',
 #' showAreaSqKm=FALSE, showPctFull=TRUE)
 #'
-#' df <- sc_nlcd(year='2001, 2006', aoi='cat,ws',
+#' df <- sc_get_nlcd(year='2001, 2006', aoi='cat,ws',
 #' comid='179,1337,1337420')
 #' }
 #' @export
 
-sc_nlcd <- function(year = '2019',
+sc_get_nlcd <- function(year = '2019',
                     comid = NULL,
                     aoi = NULL,
                     showAreaSqKm = NULL,
@@ -326,6 +332,189 @@ sc_nlcd <- function(year = '2019',
   return(final_df)
 }
 
-ignore_unused_imports <- function() {
-  curl::curl_parse_url()
+#' @rdname sc_get_nlcd
+#' @export
+#' @keywords internal
+sc_nlcd <- function(year = '2019',
+                    comid = NULL,
+                    aoi = NULL,
+                    showAreaSqKm = NULL,
+                    showPctFull = NULL,
+                    state = NULL,
+                    county = NULL,
+                    region = NULL,
+                    conus = NULL,
+                    countOnly = NULL) {
+  lifecycle::deprecate_warn("0.10.0", "sc_nlcd()", "sc_get_nlcd()")
+  sc_get_nlcd(year = '2019',
+              comid = NULL,
+              aoi = NULL,
+              showAreaSqKm = NULL,
+              showPctFull = NULL,
+              state = NULL,
+              county = NULL,
+              region = NULL,
+              conus = NULL,
+              countOnly = NULL)
+}
+
+#' @title Get NNI
+#' 
+#' @description 
+#' Function to get all NNI data available for a given year.
+#' 
+#' @author
+#' Selia Markley
+#'
+#' @param year Years(s) of NNI metrics to query.
+#' Only valid NNI years are accepted (1987:2017)
+#' Syntax: year=<year1>,<year2>
+#' 
+#' @param aoi Specify the area of interest described by a metric. By default, all available areas of interest
+#' for a given metric are returned.
+#' Syntax: areaOfInterest=<value1>,<value2>
+#' Values: catchment|watershed
+#' 
+#' @param comid Return metric information for specific COMIDs
+#' Syntax: comid=<comid1>,<comid2>
+#' 
+#' @param showAreaSqKm Return the area in square kilometers of a given area of interest.
+#' The default value is true.
+#' Values: true|false
+#' 
+#' @param showPctFull Return the pctfull for each dataset. The default value is false.
+#' Values: true|false
+#' 
+#' @param countOnly Return a CSV containing only the row count (ROWCOUNT) and the column
+#' count (COLUMNCOUNT) that the server expects to return in a request. The default value is false.
+#' Values: true|false
+#'
+#' @return A tibble of desired StreamCat metrics
+#'
+#' @examples
+#' \dontrun{
+#' df <- sc_get_nni(year='1987, 1990, 2005, 2017', aoi='cat,ws',
+#' comid=179,1337,1337420')
+#' 
+#' df <- sc_get_nni(year='2015', aoi='cat',
+#' comid='179', countOnly=TRUE)
+#' 
+#' df <- sc_get_nni(comid='179', year='2011, 2012', aoi='ws')
+#' 
+#' df <- sc_get_nni(year='2015, 2016, 2017', county='41003', aoi='ws')
+#' }
+#' #' @export
+
+sc_get_nni <- function(year, aoi = NULL, comid = NULL,
+                      showAreaSqKm = TRUE, state = NULL,
+                      county = NULL, region = NULL,conus = NULL, 
+                      showPctFull = NULL,countOnly = NULL) {
+  # year must be a character string.
+  year_chr <-  as.character(year)
+  # split multiple years supplied as a single string into
+  # a vector of years.
+  year_vec <- unlist(strsplit(x = year_chr,
+                              split = ",|, "))
+  # Vector of valid NNI years to check inputs against.
+  valid_years <- c('1987',
+                   '1988',
+                   '1989',
+                   '1990',
+                   '1991',
+                   '1992',
+                   '1993',
+                   '1994',
+                   '1995',
+                   '1996',
+                   '1997',
+                   '1998',
+                   '1999',
+                   '2000',
+                   '2001',
+                   '2002',
+                   '2003',
+                   '2004',
+                   '2005',
+                   '2006',
+                   '2007',
+                   '2008',
+                   '2009',
+                   '2010',
+                   '2011',
+                   '2012',
+                   '2013',
+                   '2014',
+                   '2015',
+                   '2016',
+                   '2017')
+  # Stop early if any of the year(s) supplied are not found in the valid
+  # years vec.
+  stopifnot(
+    "year must be a valid NNI year" = any(year_vec %in% valid_years)
+  )
+  # Vector of NNI metric names.
+  nni <- c(
+    'n_leg_',
+    'n_ags_',
+    'n_ff_',
+    'n_uf_',
+    'n_cf_',
+    'n_cr_',
+    'n_hw_',
+    'n_lw_',
+    'p_leg_',
+    'p_ags_',
+    'p_ff_',
+    'p_uf_',
+    'p_cr_',
+    'p_hw_',
+    'p_lw_'
+  )
+  # Add n_dep for available years
+  ndep_year_vec <- year_vec[!year_vec %in% c('1987', '1988', '1989')]
+  ndep_comb <- expand.grid('n_dep_', ndep_year_vec)
+  ndep_mets <- paste0(ndep_comb$Var1,
+                      ndep_comb$Var2,
+                      collapse = ",",
+                      recycle0 = TRUE)
+  # Add p_dep for available years
+  pdep_year_vec <- year_vec[!year_vec %in% c('1987', '1988', '1989', '1990', '1991', '1992', '1993', '1994', '1995', '1996',
+                                             '2014', '2015', '2016', '2017')]
+  pdep_comb <- expand.grid('p_dep_', pdep_year_vec)
+  pdep_mets <- paste0(pdep_comb$Var1,
+                      pdep_comb$Var2,
+                      collapse = ",",
+                      recycle0 = TRUE)
+  # Add n_usgsww and p_usgsww for available years
+  ww_year_vec <- year_vec[year_vec %in% c('1988', '1990', '1992', '1996', '2000', '2004', '2008', '2012')]
+  ww_comb <- expand.grid(c('p_usgsww_', 'n_usgsww_'), ww_year_vec)
+  ww_mets <- paste0(ww_comb$Var1,
+                    ww_comb$Var2,
+                    collapse = ",",
+                    recycle0 = TRUE)
+  # Create a data frame of all NNI Metric and year combinations.
+  all_comb <- expand.grid(nni, year_vec)
+  # Concatenate the NLCD metric name with the supplied year(s) to create
+  # valid metric names to submit to the API.
+  nni_mets <- paste0(all_comb$Var1,
+                     all_comb$Var2,
+                     collapse = ",",
+                     recycle0 = TRUE)
+  # Combine all NNI metrics
+  nni_mets_all <- paste0(nni_mets, ",", ndep_mets, ",", pdep_mets, ",", ww_mets)
+  
+  # Query the API.
+  final_df <- sc_get_data(
+    metric = nni_mets_all,
+    aoi = aoi,
+    comid = comid,
+    state = state,
+    county = county,
+    showAreaSqKm = showAreaSqKm,
+    showPctFull = showPctFull,
+    conus = conus,
+    countOnly = countOnly
+  )
+  # End of function. Return a data frame
+  return(final_df)
 }
diff --git a/R/sc_plot.R b/R/sc_plot.R
new file mode 100644
index 0000000..86ea89a
--- /dev/null
+++ b/R/sc_plot.R
@@ -0,0 +1,348 @@
+#' @title Plot National Nutrient Inventory data for streams
+#'
+#' @description 
+#' Function to plot time series of nitrogen and phosphorus budgets for a given stream
+#' COMID. This function allows a user to return a time series of major inputs, 
+#' outputs, and derived metrics of nitrogen and phosphorus. Plot is returned as an 
+#' object
+#' 
+#' @author 
+#' Selia Markley
+#'
+#' @param comid Identifier of stream COMID user wants to plot NNI data for. Must be a character string
+#' with the COMID digit.
+#' Syntax: com=<COMID>
+#' 
+#' @param include.nue Include time series of nitrogen use efficiency in the returned plot. 
+#' The default value is false. 
+#' Values: true|false 
+#' 
+#' @param include.inset Include inset map that shows the location of the COMID and its basin. 
+#' The default value is true. 
+#' Values: true|false 
+#'
+#' @return
+#' Return plot as an object.
+#' @export
+#'
+#' @import ggplot2
+#' @import ggpattern
+#' 
+#' @examples
+#' \dontrun{
+#' p <- sc_plotnni(comid='1337420')
+#' p <- sc_plotnni(comid='1337420', include.nue=TRUE)
+#' p <- sc_plotnni(comid='1337420', include.inset=FALSE)
+#' }
+
+
+sc_plotnni <- function(comid, include.nue = FALSE, include.inset = TRUE){
+  message("If the plot does not render to the plot window when calling the function either save the plot or resize the plot window")
+  # Get StreamCat data
+  nni <- sc_get_data(metric = 'n_dep_1990,n_ff_1990,n_uf_1990,n_lw_1990,n_hw_1990,n_ags_1990,n_cf_1990,n_cr_1990,p_cr_1990,p_lw_1990,p_hw_1990,p_uf_1990,p_ff_1990,p_ags_1990,n_dep_1991,n_ff_1991,n_uf_1991,n_lw_1991,n_hw_1991,n_ags_1991,n_cf_1991,n_cr_1991,p_cr_1991,p_lw_1991,p_hw_1991,p_uf_1991,p_ff_1991,p_ags_1991,n_dep_1992,n_ff_1992,n_uf_1992,n_lw_1992,n_hw_1992,n_ags_1992,n_cf_1992,n_cr_1992,p_cr_1992,p_lw_1992,p_hw_1992,p_uf_1992,p_ff_1992,p_ags_1992,n_dep_1993,n_ff_1993,n_uf_1993,n_lw_1993,n_hw_1993,n_ags_1993,n_cf_1993,n_cr_1993,p_cr_1993,p_lw_1993,p_hw_1993,p_uf_1993,p_ff_1993,p_ags_1993,n_dep_1994,n_ff_1994,n_uf_1994,n_lw_1994,n_hw_1994,n_ags_1994,n_cf_1994,n_cr_1994,p_cr_1994,p_lw_1994,p_hw_1994,p_uf_1994,p_ff_1994,p_ags_1994,n_dep_1995,n_ff_1995,n_uf_1995,n_lw_1995,n_hw_1995,n_ags_1995,n_cf_1995,n_cr_1995,p_cr_1995,p_lw_1995,p_hw_1995,p_uf_1995,p_ff_1995,p_ags_1995,n_dep_1996,n_ff_1996,n_uf_1996,n_lw_1996,n_hw_1996,n_ags_1996,n_cf_1996,n_cr_1996,p_cr_1996,p_lw_1996,p_hw_1996,p_uf_1996,p_ff_1996,p_ags_1996,n_dep_1997,n_ff_1997,n_uf_1997,n_lw_1997,n_hw_1997,n_ags_1997,n_cf_1997,n_cr_1997,p_cr_1997,p_lw_1997,p_hw_1997,p_uf_1997,p_ff_1997,p_ags_1997,n_dep_1998,n_ff_1998,n_uf_1998,n_lw_1998,n_hw_1998,n_ags_1998,n_cf_1998,n_cr_1998,p_cr_1998,p_lw_1998,p_hw_1998,p_uf_1998,p_ff_1998,p_ags_1998,n_dep_1999,n_ff_1999,n_uf_1999,n_lw_1999,n_hw_1999,n_ags_1999,n_cf_1999,n_cr_1999,p_cr_1999,p_lw_1999,p_hw_1999,p_uf_1999,p_ff_1999,p_ags_1999,n_dep_2000,n_ff_2000,n_uf_2000,n_lw_2000,n_hw_2000,n_ags_2000,n_cf_2000,n_cr_2000,p_cr_2000,p_lw_2000,p_hw_2000,p_uf_2000,p_ff_2000,p_ags_2000,n_dep_2001,n_ff_2001,n_uf_2001,n_lw_2001,n_hw_2001,n_ags_2001,n_cf_2001,n_cr_2001,p_cr_2001,p_lw_2001,p_hw_2001,p_uf_2001,p_ff_2001,p_ags_2001,n_dep_2002,n_ff_2002,n_uf_2002,n_lw_2002,n_hw_2002,n_ags_2002,n_cf_2002,n_cr_2002,p_cr_2002,p_lw_2002,p_hw_2002,p_uf_2002,p_ff_2002,p_ags_2002,n_dep_2003,n_ff_2003,n_uf_2003,n_lw_2003,n_hw_2003,n_ags_2003,n_cf_2003,n_cr_2003,p_cr_2003,p_lw_2003,p_hw_2003,p_uf_2003,p_ff_2003,p_ags_2003,n_dep_2004,n_ff_2004,n_uf_2004,n_lw_2004,n_hw_2004,n_ags_2004,n_cf_2004,n_cr_2004,p_cr_2004,p_lw_2004,p_hw_2004,p_uf_2004,p_ff_2004,p_ags_2004,n_dep_2005,n_ff_2005,n_uf_2005,n_lw_2005,n_hw_2005,n_ags_2005,n_cf_2005,n_cr_2005,p_cr_2005,p_lw_2005,p_hw_2005,p_uf_2005,p_ff_2005,p_ags_2005,n_dep_2006,n_ff_2006,n_uf_2006,n_lw_2006,n_hw_2006,n_ags_2006,n_cf_2006,n_cr_2006,p_cr_2006,p_lw_2006,p_hw_2006,p_uf_2006,p_ff_2006,p_ags_2006,n_dep_2007,n_ff_2007,n_uf_2007,n_lw_2007,n_hw_2007,n_ags_2007,n_cf_2007,n_cr_2007,p_cr_2007,p_lw_2007,p_hw_2007,p_uf_2007,p_ff_2007,p_ags_2007,n_dep_2008,n_ff_2008,n_uf_2008,n_lw_2008,n_hw_2008,n_ags_2008,n_cf_2008,n_cr_2008,p_cr_2008,p_lw_2008,p_hw_2008,p_uf_2008,p_ff_2008,p_ags_2008,n_dep_2009,n_ff_2009,n_uf_2009,n_lw_2009,n_hw_2009,n_ags_2009,n_cf_2009,n_cr_2009,p_cr_2009,p_lw_2009,p_hw_2009,p_uf_2009,p_ff_2009,p_ags_2009,n_dep_2010,n_ff_2010,n_uf_2010,n_lw_2010,n_hw_2010,n_ags_2010,n_cf_2010,n_cr_2010,p_cr_2010,p_lw_2010,p_hw_2010,p_uf_2010,p_ff_2010,p_ags_2010,n_dep_2011,n_ff_2011,n_uf_2011,n_lw_2011,n_hw_2011,n_ags_2011,n_cf_2011,n_cr_2011,p_cr_2011,p_lw_2011,p_hw_2011,p_uf_2011,p_ff_2011,p_ags_2011,n_dep_2012,n_ff_2012,n_uf_2012,n_lw_2012,n_hw_2012,n_ags_2012,n_cf_2012,n_cr_2012,p_cr_2012,p_lw_2012,p_hw_2012,p_uf_2012,p_ff_2012,p_ags_2012,n_dep_2013,n_ff_2013,n_uf_2013,n_lw_2013,n_hw_2013,n_ags_2013,n_cf_2013,n_cr_2013,p_cr_2013,p_lw_2013,p_hw_2013,p_uf_2013,p_ff_2013,p_ags_2013,n_dep_2014,n_ff_2014,n_uf_2014,n_lw_2014,n_hw_2014,n_ags_2014,n_cf_2014,n_cr_2014,p_cr_2014,p_lw_2014,p_hw_2014,p_uf_2014,p_ff_2014,p_ags_2014,n_dep_2015,n_ff_2015,n_uf_2015,n_lw_2015,n_hw_2015,n_ags_2015,n_cf_2015,n_cr_2015,p_cr_2015,p_lw_2015,p_hw_2015,p_uf_2015,p_ff_2015,p_ags_2015,n_dep_2016,n_ff_2016,n_uf_2016,n_lw_2016,n_hw_2016,n_ags_2016,n_cf_2016,n_cr_2016,p_cr_2016,p_lw_2016,p_hw_2016,p_uf_2016,p_ff_2016,p_ags_2016,n_dep_2017,n_ff_2017,n_uf_2017,n_lw_2017,n_hw_2017,n_ags_2017,n_cf_2017,n_cr_2017,p_cr_2017,p_lw_2017,p_hw_2017,p_uf_2017,p_ff_2017,p_ags_2017',
+                     aoi='ws',
+                     comid = comid,
+                     showAreaSqKm = FALSE,
+                     showPctFull = FALSE)
+  
+  #Create N inputs df
+  
+  nin <- nni[, grepl("^(n)", names(nni)) & !grepl("(cr)", names(nni)) & !grepl("(ags)", names(nni))]
+  
+  names(nin) <- sapply(names(nin), function(col){
+    substr(col, 3, nchar(col) -2)
+  })
+  
+  nin <- nin |>
+    tidyr::pivot_longer(
+      cols = everything(),
+      names_to = c("metric", "year"),
+      names_sep = "_",
+      values_to = "value"
+    ) |>
+    dplyr::mutate(year = as.integer(year)) |>
+    dplyr::mutate(value = value / 1000000)
+  
+  #Create P inputs df
+  
+  pin <- nni[, grepl("^(p)", names(nni)) & !grepl("(cr)", names(nni)) & !grepl("(ags)", names(nni))]
+  
+  names(pin) <- sapply(names(pin), function(col){
+    substr(col, 3, nchar(col) -2)
+  })
+  
+  pin <- pin |>
+    tidyr::pivot_longer(
+      cols = everything(),
+      names_to = c("metric", "year"),
+      names_sep = "_",
+      values_to = "value"
+    ) |>
+    dplyr::mutate(year = as.integer(year)) |>
+    dplyr::mutate(value = value / 1000000)
+  
+  #Create N dfs for lines (cr, agsur, nue)
+  
+  nlines <- nni[, grepl("^n_cr|^n_ags", names(nni))]
+  
+  names(nlines) <- sapply(names(nlines), function(col){
+    substr(col, 3, nchar(col) -2)
+  })
+  
+  nlines <- nlines |>
+    tidyr::pivot_longer(
+      cols = everything(),
+      names_to = c("metric","year"),
+      names_sep = "_",
+      values_to = "value"
+    ) |>
+    dplyr::mutate(year = as.integer(year)) |>
+    dplyr::mutate(value = value / 1000000) |>
+    tidyr::pivot_wider(
+      names_from = 'metric', 
+      values_from = 'value'
+    ) |>
+    dplyr::mutate(totag = ags + cr)  |>
+    dplyr::mutate(nue = (cr / totag) * 100) |>
+    tidyr::pivot_longer(
+      cols = !year, 
+      names_to="metric", 
+      values_to="value")
+  
+  ncrag <- nlines |>
+    dplyr::filter(metric %in% c('ags', 'cr'))
+  
+  nue <- nlines |>
+    dplyr::filter(metric == 'nue') |>
+    tidyr::pivot_wider(names_from = 'metric',
+                values_from = 'value')
+  
+  #Create P dfs for lines (cr, agsur, pue)
+  
+  plines <- nni[, grepl("^p_cr|^p_ags", names(nni))]
+  
+  names(plines) <- sapply(names(plines), function(col){
+    substr(col, 3, nchar(col) -2)
+  })
+  
+  
+  plines <- plines |>
+    tidyr::pivot_longer(
+      cols = everything(),
+      names_to = c("metric","year"),
+      names_sep = "_",
+      values_to = "value"
+    ) |>
+    dplyr::mutate(year = as.integer(year)) |>
+    dplyr::mutate(value = value / 1000000) |>
+    tidyr::pivot_wider(
+      names_from = 'metric', 
+      values_from = 'value'
+    ) |>
+    dplyr::mutate(totag = ags + cr)  |>
+    dplyr::mutate(nue = (cr / totag) * 100) |>
+    tidyr::pivot_longer(
+      cols = !year, 
+      names_to="metric", 
+      values_to="value")
+  
+  pcrag <- plines |>
+    dplyr::filter(metric %in% c('ags', 'cr'))
+  
+  pue <- plines |>
+    dplyr::filter(metric == 'nue') |>
+    tidyr::pivot_wider(names_from = 'metric',
+                values_from = 'value')
+  
+  pdf <- dplyr::bind_rows(plines, pin)
+  
+  ndf <- dplyr::bind_rows(nlines, nin)
+  
+  #create estimate column
+  knownfertyrs <- c(1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,
+                    2005,2006,2007,2008,2009,2010,2011,2012,2017)
+  nwsin <- nin |>
+    dplyr::mutate(estimated=dplyr::case_when(
+      metric == "dep" ~ FALSE,
+      metric == "hw" ~ FALSE,
+      metric == "cf" & year %in% c(1987,1992,1997,2002,2007,2012, 2017) ~ FALSE,
+      metric == "ff" & year %in% knownfertyrs ~ FALSE,
+      metric == "uf" & year %in% knownfertyrs ~ FALSE,
+      metric == "lw" & year %in% c(1987,1992,1997,2002,2007,2012,2017) ~ FALSE,
+      TRUE ~ TRUE
+    )) 
+  
+  pwsin <- pin |>
+    dplyr::filter(metric != 'cr') |>
+    dplyr::mutate(estimated=dplyr::case_when(
+      metric == "hw" ~ FALSE,
+      metric == "lw" & year %in% c(1987,1992,1997,2002,2007,2012,2017) ~ FALSE,
+      metric == "ff" & year %in% knownfertyrs ~ FALSE,
+      metric == "uf" & year %in% knownfertyrs ~ FALSE,
+      TRUE ~ TRUE
+    ))
+  
+  #get ready for plot
+  colorsn <- c('ff' = '#A3CC51', 'lw'='#B26F2C','hw'='#E51932','uf'='black', 'dep'='#6db6ff', 'cf'='#FFD700')
+  colorsp <- c('ff' = '#A3CC51', 'lw'='#B26F2C','hw'='#E51932','uf'='black')
+  
+  nwsin$metric <- factor(nwsin$metric, levels = c('uf','hw','dep','lw','cf','ff'))
+  pwsin$metric <- factor(pwsin$metric, levels = c('uf','hw','lw','ff'))
+  
+  #Get COMID location and states for inset map
+  #states
+  states <- 
+    tigris::states(cb = TRUE, progress_bar = FALSE) |>
+    dplyr::filter(!STUSPS %in% c('HI', 'PR', 'AK', 'MP', 'GU', 'AS', 'VI')) |> 
+    sf::st_transform(crs = 5070)
+  
+  #comid
+  comidint <- as.integer(comid)
+  flowline <- nhdplusTools::get_nhdplus(comid = comidint, realization = "flowline")
+  point <- sf::st_centroid(flowline)
+  
+  #create N bar plot
+  nbar <- ggplot() + 
+    ggpattern::geom_bar_pattern(data = nwsin, 
+                     aes(x=year,y=value, fill=metric, 
+                         pattern=factor(estimated, levels=c(TRUE,FALSE),
+                                        labels=c('Estimated','Non-Estimated'))),
+                     pattern = ifelse(nwsin$estimated, 'stripe','none'),
+                     pattern_color='white',
+                     pattern_density=0.05, 
+                     pattern_fill = 'white', 
+                     pattern_alpha = 0.5, 
+                     pattern_spacing=0.025,
+                     stat='identity', position='stack', 
+                     pattern_size=0.05) +
+    labs(title = 'Nitrogen (million kg/year)',
+         y = "Budget",
+         x = " ") +
+    scale_fill_manual(values=colorsn,
+                      labels = c('ff' = 'Farm Fertilizer',
+                                 'uf' = 'Urban Fertilizer',
+                                 'cf' = 'Crop N-Fixation',
+                                 'lw' = 'Livestock Manure',
+                                 'hw' = 'Human Waste',
+                                 'dep' = 'Total Deposition')) + 
+    scale_pattern_manual(name='Estimate Status',
+                         values=c('Estimated'='stripe','Non-estimated'='none')) + 
+    geom_line(data=ncrag, 
+              aes(x=year,y=value, linetype=metric), 
+              linewidth=1.25, color="black") + 
+    scale_linetype_manual(values = 
+                            c("ags"="solid", "cr"="dotted"),
+                          labels = c('ags' = 'Agricultural Surplus',
+                                     'cr' = 'Crop Removal')) +
+    guides(fill=
+             guide_legend(order=1, override.aes = list(pattern='none')), 
+           pattern=
+             guide_legend(order=2, override.aes = list(fill='grey')), 
+           linetype=
+             guide_legend(title=NULL)) + 
+    scale_x_continuous(breaks=seq(1987,2017,by=5)) + 
+    scale_color_manual(values=c("Agricultural Surplus"="black"), 
+                       guide = 'none') + 
+    theme_bw() + 
+    theme(plot.title = element_text(size=9, face="bold"), 
+          axis.title.y = element_text(size=9), 
+          legend.background = element_rect(fill="white", colour = "black"), 
+          legend.title = element_blank())
+  
+  #create p bar plot
+  pbar <- ggplot() + 
+    ggpattern::geom_bar_pattern(data = pwsin, 
+                     aes(x=year,y=value, fill=metric, 
+                         pattern=factor(estimated, levels=c(TRUE,FALSE),
+                                        labels=c('Estimated','Non-Estimated'))), 
+                     pattern = ifelse(pwsin$estimated, 'stripe','none'),
+                     pattern_color='white',
+                     pattern_density=0.05, 
+                     pattern_fill = 'white', 
+                     pattern_alpha = 0.5, 
+                     pattern_spacing=0.025,
+                     stat='identity', position='stack', pattern_size=0.05) +
+    labs(title = 'Phosphorus (million kg/year)',
+         y = "Budget",
+         x = " ") +
+    scale_fill_manual(values=colorsp) + 
+    scale_pattern_manual(name='Estimate Status',
+                         values=c('Estimated'='stripe',
+                                  'Non-estimated'='none')) + 
+    geom_line(data=pcrag, 
+              aes(x=year,y=value, 
+                  linetype=metric), 
+              linewidth=1.25, color="black") + 
+    scale_linetype_manual(values = 
+                            c("ags"="solid", "cr"="dotted")) +
+    guides(fill=
+             guide_legend(order=1, override.aes = list(pattern='none')), 
+           pattern=
+             guide_legend(order=2, override.aes = list(fill='grey')), 
+           linetype=
+             guide_legend(title=NULL)) + 
+    guides(fill="none", pattern = "none", linetype="none") +
+    scale_x_continuous(breaks=seq(1987,2017,by=5)) + 
+    scale_color_manual(values=c("Agricultural Surplus"="black"), 
+                       guide = 'none') + 
+    theme_bw() + 
+    theme(plot.title = element_text(size=9, face="bold"), 
+          axis.title.y = element_text(size=9), 
+          legend.background = element_rect(fill="white", colour = "black"), 
+          legend.title = element_blank())
+  
+  #create nue line plots
+  nue <- ggplot() + 
+    geom_line(data=nue, aes(x=year,y=nue), linewidth=1.25, color='seagreen')+
+    theme_bw() + 
+    scale_x_continuous(breaks=seq(1987,2017,by=5)) + 
+    labs(title = 'Nitrogen Use Efficiency',
+         y = "%", 
+         x=" ") + 
+    theme(plot.title = element_text(size=9, face="bold"), 
+          axis.title.x = element_text(size=9), 
+          axis.title.y = element_text(size=9)) 
+  
+  pue <- ggplot() + 
+    geom_line(data=pue, aes(x=year,y=nue, lty='Nutrient Use Efficiency'), linewidth=1.25, color="seagreen") + 
+    theme_bw() + 
+    scale_x_continuous(breaks=seq(1987,2017,by=5)) + 
+    labs(title = 'Phosphorus Use Efficiency',
+         y = "%", 
+         x="Year") + 
+    theme(plot.title = element_text(size=9, face="bold"), 
+          axis.title.x = element_text(size=9), 
+          axis.title.y = element_text(size=9, hjust=0.5),
+          legend.background = element_rect(fill="white", colour = "black"), 
+          legend.title = element_blank())  +
+    guides(fill="none", pattern = "none", linetype="none")
+  
+  #Create inset map
+  inset <- ggplot() +
+    geom_sf(data = states, color = "grey", fill = 'transparent', lwd = .2) + 
+    geom_sf(data = point, size = 3.5, color = "red") + theme_void()
+  
+  #export final figure
+  inputs <- patchwork::wrap_plots(nbar, pbar, ncol=1, guides="collect")
+  nue <- patchwork::wrap_plots(nue, pue, ncol=1, guides="collect")
+  
+  if (include.nue == TRUE){
+      timenni <- patchwork::wrap_plots(nue, inputs, ncol=2)
+    }
+  else {
+      timenni <- inputs
+    }
+
+  if (include.inset == TRUE){
+      timenni <- cowplot::plot_grid(
+        timenni, inset,
+        ncol = 1,
+        rel_heights = c(3,1)
+      )
+      return(timenni)
+    }
+  else {
+      return(timenni)
+  }
+  
+}
diff --git a/README.md b/README.md
index 373894c..381901c 100644
--- a/README.md
+++ b/README.md
@@ -5,7 +5,7 @@
 [![CRAN](http://www.r-pkg.org/badges/version/StreamCatTools)](https://cran.r-project.org/package=StreamCatTools)
 [![cran checks](https://badges.cranchecks.info/worst/StreamCatTools.svg)](https://cran.r-project.org/web/checks/check_results_StreamCatTools.html)
 [![R-CMD-check](https://github.com/USEPA/StreamCatTools/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/USEPA/StreamCatTools/actions/workflows/R-CMD-check.yaml)
-[![Downloads](https://cranlogs.r-pkg.org/badges/StreamCatTools)](https://cran.r-project.org/package=StreamCatTools)
+[![Downloads](https://cranlogs.r-pkg.org/badges/grand-total/StreamCatTools)](https://cran.r-project.org/package=StreamCatTools)
 [![CodeCov](https://img.shields.io/badge/test%20coverage-86.9%25-388600.svg)](https://img.shields.io/badge/test%20coverage-86.9%25-388600.svg)
 <!-- badges: end -->
 
diff --git a/tests/testthat/test-sc_nni.R b/tests/testthat/test-sc_nni.R
new file mode 100644
index 0000000..7b0698e
--- /dev/null
+++ b/tests/testthat/test-sc_nni.R
@@ -0,0 +1,19 @@
+context("Test that sc_nlcd is pulling in StreamCat API data")
+
+
+test_that("sc_get_nni for a sample COMID returns a data frame", {
+  testthat::skip_on_cran()
+  df <- sc_get_nni(year='1987, 1990, 2005, 2017', aoi='cat',
+                comid='1337420')
+  expect_true(exists("df"))
+  expect_equal(nrow(df), 1)
+  expect_equal(ncol(df), 71)
+})
+
+test_that("sc_get_data for a county and ws metrics returns a data frame", {
+  testthat::skip_on_cran()
+  df <- sc_get_nni(year='2015, 2016, 2017', county='41003', aoi='ws')
+  expect_true(exists("df"))
+  expect_equal(nrow(df), 632)
+  expect_equal(ncol(df), 53)
+})
diff --git a/tests/testthat/test-sc_plot.R b/tests/testthat/test-sc_plot.R
new file mode 100644
index 0000000..4a4587f
--- /dev/null
+++ b/tests/testthat/test-sc_plot.R
@@ -0,0 +1,16 @@
+context("Test that sc_plot is creating plot object")
+
+
+test_that("sc_plot is producing basic NNI plot", {
+  testthat::skip_on_cran()
+  p <- sc_plotnni(comid='1337420')
+  expect_true(exists("p"))
+  expect_true(class(p)[1]=="ggplot2::ggplot")
+})
+
+test_that("sc_plot is producing basic NNI plot", {
+  testthat::skip_on_cran()
+  p <- sc_plotnni(comid='1337420', include.nue=TRUE)
+  expect_true(exists("p"))
+  expect_true(class(p)[1]=="ggplot2::ggplot")
+})
diff --git a/vignettes/Articles/Applications.Rmd b/vignettes/Articles/Applications.Rmd
index 9e497bc..4386224 100644
--- a/vignettes/Articles/Applications.Rmd
+++ b/vignettes/Articles/Applications.Rmd
@@ -96,8 +96,8 @@ mmi <- readr::read_csv("https://www.epa.gov/sites/production/files/2015-09/bentc
 
 # join mmi to NARS info data frame with StreamCat PctCrop metric
 nrsa_sf <- dplyr::left_join(nrsa_sf, mmi[,c('SITE_ID','BENT_MMI_COND')], by='SITE_ID')
-bxplt <- nrsa_sf %>% 
-  tidyr::drop_na(BENT_MMI_COND) %>%
+bxplt <- nrsa_sf |> 
+  tidyr::drop_na(BENT_MMI_COND)  |> 
   ggplot2::ggplot(aes(x=pctcrop2006ws, y=BENT_MMI_COND))+
   ggplot2::geom_boxplot()+
   ggplot2::ggtitle('NRSA Benthic MMI versus % Crop in Watershed from 2006 NLCD')
diff --git a/vignettes/Articles/NNI.Rmd b/vignettes/Articles/NNI.Rmd
new file mode 100644
index 0000000..0cd94fd
--- /dev/null
+++ b/vignettes/Articles/NNI.Rmd
@@ -0,0 +1,261 @@
+---
+title: "National Nutrient Inventory (NNI) Data in StreamCat"
+output: 
+  html_document:
+    theme: flatly
+    keep_md: yes
+    number_sections: true
+    highlighted: default 
+    toc: yes
+    toc_float:
+      collapsed: no
+      smooth_scroll: no
+    toc_depth: 2
+vignette: >
+  %\VignetteIndexEntry{Applications withs StreamCatTools}
+  %\VignetteEncoding{UTF-8}{inputenc}
+  %\VignetteEngine{knitr::rmarkdown}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>", 
+  warning = FALSE, 
+  message = FALSE
+)
+```
+
+## Use sc_plotnni and lc_plotnni
+In the following examples, we use sc_plotnni and lc_plotnni functions to plot annual time series 
+of nitrogen and phosphorus budget data for a given watershed. These data are plotted from 1990-2017. Hatching on the bars
+is reflective of data linearly interpolated between agricultural census years.
+
+Mississippi-Atchafalaya River Basin
+```{r, warning = FALSE, message = FALSE}
+library(StreamCatTools)
+library(ggplot2)
+library(ggpattern)
+com <- '22812041'
+sc_plotnni(comid = com, include.nue = TRUE)
+```
+
+Neuse River, North Carolina
+```{r, warning=FALSE, message=FALSE}
+com <- '10975909'
+sc_plotnni(comid = com)
+```
+
+Lake Kanasatka, New Hampshire
+```{r,warning=FALSE, message=FALSE}
+com = '6738112'
+lc_plotnni(comid = com)
+```
+## Use sc_get_nni and lc_get_nni to pull NNI budget data
+In these examples, we access all available NNI metrics for user-defined years. Most NNI metrics are available from 1987
+to 2017, except for nitrogen deposition, phosphorus deposition, and nitrogen and phosphorus point source loads. These metrics
+are returned for available years. 
+
+## Cuyahoga River
+Return dataframe of NNI metrics for Cuyahoga River
+```{r,warning=FALSE, message=FALSE}
+com <- 15588532
+years <- '1987, 1992, 1997, 2002, 2007, 2012, 2017'
+cuyahoga <- sc_get_nni(year=years,
+           comid=com,
+           aoi='ws',
+           showAreaSqKm=FALSE)
+colnames(cuyahoga)
+```
+Transform sc_get_nni dataframe into long dataframe
+```{r,warning=FALSE, message=FALSE}
+nin <- cuyahoga[, grepl("^(n)", names(cuyahoga)) & 
+                  !grepl("(cr)", names(cuyahoga)) & 
+                  !grepl("(ags)", names(cuyahoga)) & 
+                  !grepl("(leg)", names(cuyahoga))]
+  
+names(nin) <- sapply(names(nin), function(col){
+  substr(col, 3, nchar(col) -2)
+})
+  
+nin <- nin |>
+  tidyr::pivot_longer(
+    cols = everything(),
+    names_to = c("metric", "year"),
+    names_sep = "_",
+    values_to = "value"
+  ) |>
+  dplyr::mutate(year = as.integer(year))
+
+pin <- cuyahoga[, grepl("^(p)", names(cuyahoga)) & 
+                  !grepl("(cr)", names(cuyahoga)) & 
+                  !grepl("(ags)", names(cuyahoga)) & 
+                  !grepl("(leg)", names(cuyahoga))]
+  
+names(pin) <- sapply(names(pin), function(col){
+  substr(col, 3, nchar(col) -2)
+})
+  
+pin <- pin |>
+  tidyr::pivot_longer(
+    cols = everything(),
+    names_to = c("metric", "year"),
+    names_sep = "_",
+    values_to = "value"
+  ) |>
+  dplyr::mutate(year = as.integer(year))
+```
+
+Pie Chart of Sources of N and P to the Cuyahoga River in 2012
+```{r,warning=FALSE, message=FALSE}
+n2012 <- nin |> 
+  dplyr::filter(year == '2012')
+
+p2012 <- pin |>
+  dplyr::filter(year == '2012')
+
+colorsn <- c('ff' = '#A3CC51', 'lw'='#B26F2C','hw'='#E51932','uf'='black', 'dep'='#6db6ff', 'cf'='#FFD700', 'usgsww' = 'pink')
+colorsp <- c('ff' = '#A3CC51', 'lw'='#B26F2C','hw'='#E51932','uf'='black', 'usgsww' = 'pink')
+
+n <- ggplot(nin, aes(x = "", y = value, fill = metric)) +
+  geom_bar(stat = "identity", width = 1) +
+  coord_polar(theta = "y") +
+  scale_fill_manual(values = colorsn, labels = 
+                       c('ff' = 'Farm Fertilizer',
+                         'uf' = 'Urban Fertilizer',
+                         'hw' = 'Human Waste',
+                         'lw' = 'Livestock Waste',
+                         'dep' = 'Total Deposition',
+                         'usgsww' = 'Point Source Loads',
+                         'cf' = 'Crop N-Fixation')
+                     ) +
+  theme_void() +
+  labs(title = "Nitrogen",
+       fill = 'Source')
+
+p <- ggplot(pin, aes(x = "", y = value, fill = metric)) +
+  geom_bar(stat = "identity", width = 1) +
+  coord_polar(theta = "y") +
+  scale_fill_manual(values = colorsp) +
+  theme_void() +
+  labs(title = "Phosphorus") + 
+  guides(fill = 'none')
+
+n + p + patchwork::plot_layout(guides = "collect")
+```
+
+Percent change of N and P inputs from 1992 to 2017
+```{r,warning=FALSE, message=FALSE}
+nperc <- nin |>
+  dplyr::filter(year %in% c('1992', '2017')) |>
+  tidyr::pivot_wider(values_from = 'value',
+              names_from = 'year')
+nperc$perchange <- (nperc$`2017` - nperc$`1992`) / nperc$`1992` * 100
+nperc <- nperc |>
+  dplyr::select(metric, perchange) |>
+  dplyr::filter(metric != 'usgsww') |>
+  dplyr::rename('Source' = 'metric', 'Nitrogen Percent Change' = 'perchange')
+
+pperc <- pin |>
+  dplyr::filter(year %in% c('1992', '2017')) |>
+  tidyr::pivot_wider(values_from = 'value',
+              names_from = 'year')
+pperc$perchange <- (pperc$`2017` - pperc$`1992`) / pperc$`1992` * 100
+pperc <- pperc |>
+  dplyr::select(metric, perchange) |>
+  dplyr::filter(metric != 'usgsww') |>
+  dplyr::rename('Source' = 'metric', 'Phosphorus Percent Change' = 'perchange')
+
+perc <- dplyr::left_join(nperc, pperc, by = 'Source')
+
+perc <- perc |> dplyr::mutate(
+  Source = dplyr::case_when(
+    Source == 'ff' ~ 'Farm Fertilizer',
+    Source == 'uf' ~ 'Urban Fertilizer',
+    Source == 'lw' ~ 'Livestock Waste',
+    Source == 'cf' ~ 'Crop N-Fixation',
+    Source == 'dep' ~ 'Atmos. Deposition',
+    Source == 'hw' ~ 'Human Waste'
+  )
+)
+tibble::tibble(perc)
+```
+
+Local catchment and watershed N and P inputs to the Platte River Basin
+```{r,warning=FALSE, message=FALSE}
+com <- 17416474
+tibble::tibble(sc_get_nni(year='2017', comid=com, aoi='cat,ws'))
+```
+
+
+## Plot a single NNI metric for a given watershed
+In this example we access a single National Nutrient Inventory (NNI) metric for the Calapooia River basin using the `sc_get_data` function. We use the `nhdplusTools` library to pull in flowlines and the watershed boundary for the Calapooia River basin, plot the selected NNI metric for the Calapooia River and show the watershed.
+```{r farmn, results='hide'}
+start_comid = 23763517
+nldi_feature <- list(featureSource = "comid", featureID = start_comid)
+
+flowline_nldi <- nhdplusTools::navigate_nldi(nldi_feature, mode = "UT", data_source = "flowlines", distance=5000)
+
+# get StreamCat metrics
+comids <- paste(as.integer(flowline_nldi$UT_flowlines$nhdplus_comid), collapse=",",sep="")
+
+df <- sc_get_data(metric='n_ff_2016', aoi='cat', comid=comids, showAreaSqKm=TRUE)
+
+flowline_nldi <- flowline_nldi$UT_flowlines
+flowline_nldi$Farm_Nitrogon_2016 <- df$n_ff_2016cat[match(flowline_nldi$nhdplus_comid, df$comid)]
+
+basin <- nhdplusTools::get_nldi_basin(nldi_feature = nldi_feature)
+```
+
+## Map the Results
+```{r map1, results='hide', message=FALSE}
+library(ggplot2)
+library(ggspatial)
+flowline_nldi |> 
+  ggplot() + geom_sf(aes(colour = Farm_Nitrogon_2016)) + 
+  scale_y_continuous() +
+  scale_color_distiller(palette = "Spectral") +
+  labs(color = "Kg Nitrogen") + 
+  theme_minimal(12) + 
+  ggtitle('Farm Nitrogen at the Catchment Scale for \nthe Calapooia River Basin 2016')
+```
+
+## Look at change through time Calapooia Farm Nitrogen
+```{r change, results='hide', message=FALSE}
+df1 <- sc_get_data(metric='n_ff_1987', aoi='cat', comid=comids)
+df2 <- sc_get_data(metric='n_ff_2017', aoi='cat', comid=comids)
+df2$n_ff_1987cat <- df1$n_ff_1987cat[match(df2$comid, df1$comid)]
+df2$Farm_Nitrogen_Difference <- df2$n_ff_2017cat-df2$n_ff_1987cat
+flowline_nldi$Farm_Nitrogen_Difference <- df2$Farm_Nitrogen_Difference[match(flowline_nldi$nhdplus_comid, df2$comid)]
+```
+
+## Map the Results
+```{r map2, results='hide', message=FALSE}
+flowline_nldi |> 
+  ggplot() + geom_sf(aes(colour = Farm_Nitrogen_Difference)) + 
+  scale_y_continuous() + 
+  labs(color = "Kg Nitrogen")+ 
+  scale_color_distiller(palette = "Spectral") +
+  theme_minimal(12) + 
+  ggtitle('Farm Nitrogen Change at the Catchment Scale for \nthe Calapooia River Basin from 1987 to 2017')
+```
+
+## Crop Fixation for the Calapooia Watershed
+```{r crop, results='hide', message='false'}
+options(scipen=3)
+# get StreamCat metrics
+df <- sc_get_data(metric='n_cf_2017', aoi='ws', comid=comids)
+
+flowline_nldi$CropNFixation <- df$n_cf_2017ws[match(flowline_nldi$nhdplus_comid, df$comid)]
+```
+
+## Map the Results
+```{r map3, results='hide', message=FALSE}
+flowline_nldi |> 
+  ggplot() + geom_sf(aes(colour = CropNFixation)) + 
+  scale_y_continuous() +
+  scale_color_distiller(palette = "Spectral") +
+  labs(colour = "Kg Nitrogen") + 
+  theme_minimal(12) + 
+  ggtitle('Watershed Crop Nitrogen Fixation \nfor the Calapooia River Basin 2017')
+```
\ No newline at end of file