various updates

Author: logan-berner

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: LandsatTS
 Title: An R package to facilitate retrieval, cleaning, cross-calibration, and phenological modeling of Landsat time-series data
-Version: 1.0.2
+Version: 1.1.0
 Authors@R:
     c(
     person(given = "Logan T.",
@@ -24,7 +24,7 @@ Authors@R:
            email = "",
            comment = c(ORCID = "0000-0002-6326-4308"))
            )
-Description: This software package facilitates sample-based time series analysis of surface reflectance and spectral indices derived from sensors on the Landsat satellites. The package includes functions that enable extraction of the full Landsat record for point sample locations or small study regions using Google Earth Engine directly accessed from R. Moreover, the package includes functions for (1) rigorous data cleaning, (2) cross-sensor calibration with machine learning, (3) phenological modeling, and (4) time series analysis.
+Description: This software package facilitates sample-based time series analysis of surface reflectance and spectral indices derived from sensors on the Landsat satellites. The package includes functions that enable extraction of the full Landsat record for point sample locations or small study regions using Google Earth Engine directly accessed from R. Moreover, the package includes functions for (1) data cleaning, (2) cross-sensor calibration, (3) phenological modeling, and (4) time series analysis.
 License: Modified MIT
 Encoding: UTF-8
 LazyData: true

LICENSE

@@ -1,12 +1,11 @@
 MIT License
 
-Copyright (c) 2021 by
-    Authors:  Logan T. Berner, Jakob J. Assmann, Richard Massey, Signe Normand
-              and Scott Goetz
-    Citation: Logan T. Berner, Jakob J. Assmann, Richard Massey, Signe Normand
-              and Scott J. Goetz. 2021. lsatTS - an R package for deriving
-              vegetation greenness time series using Landsat satellite data.
-              https://github.com/logan-berner/lsatTS
+Copyright (c) 2022 by
+    Authors:  Logan T. Berner, Jakob J. Assmann, Signe Normand, and Scott Goetz
+    Citation: Logan T. Berner, Jakob J. Assmann, Signe Normand, and Scott J. Goetz. 
+    2022. LandsatTS: an R package to facilitate retrieval, cleaning, cross-calibration, 
+    and phenological modeling of Landsat time-series data. 
+    https://github.com/logan-berner/LandsatTS
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

NAMESPACE

@@ -12,6 +12,7 @@ export(lsat_fit_phenological_curves)
 export(lsat_format_data)
 export(lsat_get_pixel_centers)
 export(lsat_neighborhood_mean)
+export(lsat_plot_trend_hist)
 export(lsat_summarize_data)
 export(lsat_summarize_growing_seasons)
 import(data.table)

NEWS.md

@@ -1,11 +1,14 @@
 # LandsatTS changelog
 
-## 1.0.2 (released: 30 November 2022)
-- Changed package name from lsatTS to LandsatTS.
-- Added function lsat_calibrate_poly().
+## 1.1.0 (released: 2022-12-07)
+- Changed package name from lsatTS to LandsatTS
+- Added function lsat_calibrate_poly()
+- Added function lsat_plot_trend_hist()
+- Revised function lsat_fit_phenological_curves() with improve outlier screening
+- Revised function lsat_calc_trend() so it no longer creates figures
 
-## 1.0.1 (released: 2 June 2022)
+## 1.0.1 (released: 2022-6-2)
 - Bug fix to address error in lsat_get_pixel_centers() caused by change in the URL for the WRS2 geometries hosted by USGS. See Issue #35: https://github.com/logan-berner/lsatTS/issues/35
 
-## 1.0.0 (released: 12 April 2022)
+## 1.0.0 (released: 2022-4-12)
 - First public release of lsatTS.

R/lsat_calc_trend.R

@@ -1,7 +1,7 @@
 #' Calculate non-parametric vegetation greenness trends
 #'
 #' @description This function evaluates and summarizes interannual trends in 
-#' vegetation greenness for sample sites over a user-specificed time period. 
+#' vegetation greenness for sample sites over a user-specified time period. 
 #' Potential interannual trends in vegetation greenness are assessed using 
 #' Mann-Kendall trend tests and Theil-Sen slope indicators after prewhitening
 #' each time series. This trend assessment relies on the zyp.yuepilon() function 
@@ -15,13 +15,11 @@
 #' @param nyr.min.frac Fraction of years within the time period for which observations 
 #'     must be available if a trend is to be computed.
 #' @param sig A p-value significance cutoff used to categories trends (e.g., 0.10)
-#' @param legend.position Legend position for output plot, specified as a vector with x and y values
-#'    ranging from 0 to 1.   
-#' @param legend.direction Legend direction for output plot, either "horizontal" or "vertical"
-#' @return Data.table with summary of temporal trends by site and a multi-panel figure with:
-#'     (1) a histogram of relative changes in vegetation greenness among sample sites and
-#'     (2) a time-series plot of mean vegetation greenness for sample sites grouped
-#'         by trend category.     
+#' @return A list that includes: 
+#'     (1) a summary message about the mean relative change across sample sites;
+#'     (2) a data.table summarizing the number and percentage of sites that
+#'     fall into each trend category;
+#'     (3) a data.table with trend statistics for each sample site.
 #' 
 #' @export lsat_calc_trend
 #' @import data.table
@@ -41,9 +39,7 @@ lsat_calc_trend <- function(dt,
                             yrs, 
                             yr.tolerance = 1, 
                             nyr.min.frac = 0.66, 
-                            sig = 0.10, 
-                            legend.position=c(0.8, 0.2), 
-                            legend.direction = 'horizontal'){
+                            sig = 0.10){
 
   dt <- data.table::data.table(dt)
   data.table::setnames(dt, si, 'si')
@@ -58,7 +54,7 @@ lsat_calc_trend <- function(dt,
   site.smry[, trend.period := paste0(min(yrs),'to',max(yrs))]
 
   # note which si is being used
-  site.smry[, si.name := si]
+  site.smry[, si := si]
 
   # identify sites with observations within the specified tolerance of first and last years
   site.smry[, ':='(first.yr.abs.dif = abs(first.yr - min(yrs)),
@@ -94,57 +90,6 @@ lsat_calc_trend <- function(dt,
   trend.dt[pval <= sig & slope < 0, trend.cat := 'browning']
   trend.dt[pval > sig, trend.cat := 'no_trend']
 
-  # histogram of vegetation greenness trends
-  fig1 <- ggplot2::ggplot(trend.dt, ggplot2::aes(total.change.pcnt, fill=after_stat(x))) +
-    ggplot2::geom_histogram(bins = 50, size = 0.25, color = 'gray20') +
-    ggplot2::scale_fill_gradient2(low="darkgoldenrod4", mid='white', high="darkgreen", 
-                                  limits = c(-50,50), midpoint = 0) +
-    ggplot2::labs(y = 'Number of sample sites', 
-                  x = paste0("Relative change in Landsat ", 
-                             gsub('.MAX', 'max', toupper(si)),
-                             ' from ', min(yrs), ' to ', max(yrs), ' (%)')) +
-    ggplot2::theme_bw() + 
-    ggplot2::theme(legend.position = 'none', 
-                   axis.text=ggplot2::element_text(size=12), 
-                   axis.title=ggplot2::element_text(size=14)) + 
-    ggplot2::xlim(-50, 50)
-  
-  # Create time series figure for each trend class
-  dt <- dt[trend.dt, on = 'sample.id']
-  
-  trend.cls.yrly.dt <- dt[, .(si.avg = mean(si, na.rm = T), si.sd = sd(si, na.rm = T),
-                              n = .N), by = c('trend.cat','year')]
-  
-  trend.cls.yrly.dt[, si.se := si.sd/sqrt(n)]
-  trend.cls.yrly.dt[is.na(si.se), si.se := 0]
-  
-  trend.cls.yrly.dt[, trend.cat := factor(trend.cat, levels = c('browning','no_trend','greening'), 
-                                          labels = c('browning','no trend','greening'))]
-  
-  trend.cls.yrly.dt[, trend.cat := droplevels.factor(trend.cat)]
-  
-  trend.cols <- data.table(trend.cat = c('browning','no trend','greening'),
-                           cols = c('darkgoldenrod4','ivory3','darkgreen'))
-  trend.cols <- trend.cols[trend.cat %in% trend.cls.yrly.dt$trend.cat]
-  
-  fig2 <- ggplot2::ggplot(trend.cls.yrly.dt, 
-                          ggplot2::aes(year, si.avg, group = trend.cat, color = trend.cat)) + 
-    ggplot2::labs(y=paste0('Mean Landsat ', gsub('.MAX', 'max', toupper(si))), x='Year') + 
-    ggplot2::geom_ribbon(ggplot2::aes(ymin = si.avg-si.se, ymax = si.avg + si.se, fill=trend.cat), 
-                         alpha=0.3, linetype=0)+
-    ggplot2::geom_line(ggplot2::aes(color = trend.cat), alpha = 1, linewidth=1) + 
-    ggplot2::scale_fill_manual(values = trend.cols$cols, name = 'Trend class') + 
-    ggplot2::scale_color_manual(values = trend.cols$cols, name = 'Trend class')+
-    ggplot2::theme_bw() +
-    ggplot2::theme(legend.position=legend.position, legend.direction=legend.direction,  
-                   axis.text=ggplot2::element_text(size=12), 
-                   axis.title=ggplot2::element_text(size=14),
-                   plot.title=ggplot2::element_text(hjust = 0.5))
-  
-  # combo figure
-  fig <- ggpubr::ggarrange(fig1, fig2, ncol = 1, nrow = 2, 
-                           labels = c('(a)','(b)'), vjust=0.9, hjust = -0.1)
-  
   # create output message about average relative change 
   avg <- round(mean(trend.dt$total.change.pcnt),2)
   std <- round(sd(trend.dt$total.change.pcnt),2)
@@ -159,12 +104,14 @@ lsat_calc_trend <- function(dt,
   names(trend.freqs.dt) <- c('Trend category','Number of sites','Percent of sites')
 
   # output
-  print(fig)
-  print(print.avg.msg)
-  print(trend.freqs.dt)
+  output.message <- list()
+  output.message$message <- print.avg.msg
+  output.message$trend.frequencies <- trend.freqs.dt
+  print(output.message)
   trend.dt
 }
 
+
 calc.trends <- function(x,y){
   xx <- zyp::zyp.yuepilon(y,x) ## note the order of x and y are switched in this call!!!
   return(data.table(slope=round(xx['trend'],5), 

R/lsat_plot_trend_hist.R

@@ -0,0 +1,46 @@
+#' Create a histogram summarizing relative temporal changes in a spectral index
+#' across all sample sites.   
+#'
+#' @param dt A data.table output from lsat_calc_trend()
+#' @param xlim Numeric vector specifying the minimum and maximum values for the
+#' histogram x-axis. 
+#'
+#' @return A histogram generated by ggplot2
+#' @export lsat_plot_trend_hist
+#'
+#' @examples
+#' data(lsat.example.dt)
+#' lsat.dt <- lsat_format_data(lsat.example.dt)
+#' lsat.dt <- lsat_clean_data(lsat.dt)
+#' lsat.dt <- lsat_calc_spectral_index(lsat.dt, 'ndvi')
+#' # lsat.dt <- lsat_calibrate_rf(lsat.dt, band.or.si = 'ndvi', write.output = F)
+#' lsat.pheno.dt <- lsat_fit_phenological_curves(lsat.dt, si = 'ndvi') 
+#' lsat.gs.dt <- lsat_summarize_growing_seasons(lsat.pheno.dt, si = 'ndvi')
+#' lsat.trend.dt <- lsat_calc_trend(lsat.gs.dt, si = 'ndvi.max', yrs = 2000:2020)
+#' lsat_plot_trend_hist(lsat.trend.dt)
+
+lsat_plot_trend_hist <- function(dt, xlim = c(-30,30)){
+  
+  dt <- data.table::data.table(dt)
+  si <- first(dt$si)
+  first.yr <- first(dt$first.yr)
+  last.yr <-  first(dt$last.yr)
+  
+  # histogram of vegetation greenness trends
+  hist <- ggplot2::ggplot(dt, ggplot2::aes(total.change.pcnt, fill=after_stat(x))) +
+  ggplot2::geom_histogram(bins = 50, size = 0.25, color = 'gray20') +
+  ggplot2::scale_fill_gradient2(low="darkgoldenrod4", mid='white', high="darkgreen", 
+                                limits = c(-50,50), midpoint = 0) +
+  ggplot2::labs(y = 'Number of sample sites', 
+                x = paste0("Relative change in Landsat ", 
+                           gsub('.MAX', 'max', toupper(si)),
+                           ' from ', first.yr, ' to ', last.yr, ' (%)')) +
+  ggplot2::theme_bw() + 
+  ggplot2::theme(legend.position = 'none', 
+                 axis.text=ggplot2::element_text(size=12), 
+                 axis.title=ggplot2::element_text(size=14)) + 
+  ggplot2::xlim(xlim[1], xlim[2])
+
+  hist
+
+}

data-raw/generate_external_data_for_itex_examples.R

@@ -19,17 +19,17 @@ itex.sites.dt <- itex.dt %>% rename(sample_id = study_area) %>%
   group_by(sample_id) %>%
   top_n(n = 1) %>% as.data.table()
 
-itex.sites.dt
+# exclude the Bogong site located in Australia
+itex.sites.dt <- itex.sites.dt[sample_id != 'Bogong']
 save(itex.sites.dt, file="data/itex.sites.dt.RData")
 
-itex.sites.sf <- itex.sites.dt %>% as.data.frame() %>% st_as_sf(coords = c("long", "lat"), crs = 4326)
-save(itex.sites.sf, file="data/itex.sites.sf.RData")
-
-jpeg('man/manuscript/figures/itex_site_map.jpg', width = 7, height = 5, units = 'in', res = 300)
-plot(st_geometry(ne_countries(returnclass = "sf")))
-plot(st_geometry(itex.sites.sf), col = "black", bg = 'red', add = T, pch = 21)
-dev.off()
+# itex.sites.sf <- itex.sites.dt %>% as.data.frame() %>% st_as_sf(coords = c("long", "lat"), crs = 4326)
+# save(itex.sites.sf, file="data/itex.sites.sf.RData")
 
+# jpeg('man/manuscript/figures/itex_site_map.jpg', width = 7, height = 5, units = 'in', res = 300)
+# plot(st_geometry(ne_countries(returnclass = "sf")))
+# plot(st_geometry(itex.sites.sf), col = "black", bg = 'red', add = T, pch = 21)
+# dev.off()
 
 
 # EXPORT TIME SERIES ==============================================================
@@ -44,5 +44,6 @@ task_list <- lsat_export_ts(pixel_coords_sf = itex.sites.sf,
 
 
 # GRAB DATA FROM GOOGLE DRIVE AND MOVE TO LOCAL DIRECTORY =========================
-itex.lsat.dt <- fread('C:/Users/Logan/My Drive/gee_export/itex_sites_chunk_1.csv')
-save(itex.lsat.dt, file="data/itex.lsat.dt.RData")
+itex.lsat.dt <- fread('C:/Users/Logan/My Drive/gee_export/itex_sites__chunk_1.csv')
+itex.lsat.dt <- itex.lsat.dt[sample_id != "Bogong"]
+save(itex.lsat.dt, file="data/itex.lsat.dt.RData")