Nothing
#### Satellite data ----
#' Auxiliary unexported and \sQuote{hidden} functions
#' Retrieve cloud polygons from the zipped Sentinel 2 file
#' @name .pa_get_cloud_polygon
#' @description This function is used to retrieve the cloud mask from different sentinel products.
#' This is used in several steps to check for clouds over an area of interest.
#' @param sentinel_product the file path to a zipper sentinel file
#' @return a polygon of the clouds in the sentinel product
#' @noRd
.pa_get_cloud_polygon <- function(sentinel_product) {
imgList <- utils::unzip(sentinel_product, list = TRUE)
cloud_mask <- grep('B00\\.jp2|B00\\.gml', imgList[[1]], ignore.case = TRUE, value = TRUE)
temporary_dir <- tempdir(check = TRUE)
utils::unzip(sentinel_product, overwrite = TRUE, exdir = temporary_dir)
## The cloud mask can be in the gml format or jp2
## Here we check for the format because they have to be processed differently
if (any(grepl('\\.gml', cloud_mask))) {
doc <- XML::xmlTreeParse(file.path(temporary_dir, cloud_mask),
asTree = T)
tolist <- rapply(doc, as.list)
coords <- grep('^[0-9]{4,} .* [0-9]{4,}$', tolist, perl = TRUE, value = TRUE)
if (length(coords) < 1){
folders_to_delete <- dir(temporary_dir, pattern = '\\.SAFE')
unlink(paste0(normalizePath(temporary_dir), "/", folders_to_delete), recursive = TRUE)
return(NULL)}
img_crs <- grep('crs:', tolist, perl = TRUE, value = TRUE)
img_crs <- unlist(strsplit(img_crs, ':'))
img_crs <- as.numeric(img_crs[length(img_crs)])
res <- c()
for (i in 1:length(coords)){
numerical_coordinates <- as.numeric(unlist(strsplit(coords[i], ' ')))
coordinate_list <- split(numerical_coordinates, rep(c(1,2)/(length(numerical_coordinates) / 2 ) ))
names(coordinate_list) <- c('lat','long')
cloud_poly <- sf::st_polygon(list(matrix(c(coordinate_list$lat, coordinate_list$long), 2)),
dim = 'XY')
cloud_poly <- sf::st_sfc(cloud_poly)
cloud_poly <- sf::st_set_crs(cloud_poly, img_crs)
res <- c(res, cloud_poly)
}
res <- sf::st_sfc(res)
res <- sf::st_set_crs(res, img_crs)
folders_to_delete <- dir(temporary_dir, pattern = '\\.SAFE')
}
if (any(grepl('\\.jp2', cloud_mask))) {
clouds <- stars::read_stars(file.path(temporary_dir, cloud_mask))
clouds <- split(clouds, 'band')
clouds <- clouds[1] + clouds[2]
clouds[clouds > 0] <- 1
res <- clouds
}
return(res)
}
#' Retrieve the Data Space token
#' @name .pa_get_dataspace_token
#' @description
#' This function retrieves the security token from Data Space. This is used in the
#' functions that download data from Data Space.
#' @param username username used to authenticate the request in Data Space
#' @param password password used to authenticate the request in Data Space
#' @noRd
.pa_get_dataspace_token <- function(username, password,
oauth = FALSE
) {
url <- "https://identity.dataspace.copernicus.eu/auth/realms/CDSE/protocol/openid-connect/token"
token_data <- list( client_id = "cdse-public", username = username,
password = password, grant_type = "password")
if (oauth) {
token_data <- list( client_id = username, client_secret = password,
grant_type = "client_credentials")
}
resp <- httr::POST( url,
body = token_data,
encode = 'form')
token <- httr::content(resp)
token$timestamp <- Sys.time()
if(is.null(token$access_token)) {
stop('There was a problem when authenticating your credentials.
Please make sure your user name and password are correct.')
}
return(token)
}
#' Crop a Sentinel 2 product to the area of interest
#' @name .pa_crop_s2_to_aoi
#' @description
#' This function crops all the raster files within the Sentinel 2 zipped file to
#' an area of interest
#' @param satellite_images list of file paths to be cropped to aoi
#' @param aoi sf object to which the file raster files will be cropped
#' @param overwrite whether to overwrite the original file
#' @noRd
.pa_crop_s2_to_aoi <- function(satellite.images,
aoi,
overwrite = TRUE){
for (sat.img in satellite.images) {
image.indices <- .pa_select_s2_files(sat.img, which = 'images')
temporary.dir <- tempdir(check = TRUE)
bname <- basename(sat.img)
bname <- strsplit(bname, '\\.')[[1]][1]
temporary.dir <- file.path(temporary.dir, bname)
reduced.dir <- file.path(temporary.dir, 'reduced/IMG_DATA')
dir.create(temporary.dir, showWarnings = FALSE, recursive = TRUE)
imgList <- utils::unzip(sat.img[[1]],
list = TRUE)
files.tgt <- .pa_select_s2_files(sat.img[[1]])
files.out <- sapply(files.tgt, function(x) grep(x, imgList[[1]], value = TRUE))
utils::unzip(sat.img[[1]],
files = files.out,
overwrite = TRUE,
exdir = temporary.dir,
junkpaths = TRUE)
cmi <- grep('B00', image.indices)
cloud.mask <- image.indices[cmi]
image.indices <- image.indices[-cmi]
dir.create(reduced.dir, recursive = TRUE, showWarnings = FALSE)
for (band in image.indices) {
band.img <- stars::read_stars(file.path(temporary.dir, band))
boundary <- sf::st_geometry(sf::st_transform(aoi, sf::st_crs(band.img)))
band.img <- sf::st_crop(band.img, boundary, mask = TRUE)
band <- gsub('\\.jp2', '\\.tif', band)
img.out <- file.path(reduced.dir, band)
stars::write_stars(band.img,
img.out)
}
if (length(cloud.mask) > 0){
file.copy(file.path(temporary.dir, cloud.mask),
file.path(reduced.dir, cloud.mask))
}
metadata <- .pa_select_s2_files(fpath = sat.img, which = 'metadata')
if (length(metadata) > 0){
file.copy(file.path(temporary.dir, metadata),
file.path(reduced.dir, metadata))
}
files.to.zip <- list.files(reduced.dir, recursive = TRUE)
if(overwrite){
wrt.path <- sat.img[[1]]
unlink(wrt.path)
}else{
bname <- basename(sat.img[[1]])
f <- strsplit(bname, '\\.')[[1]][1]
extension <- strsplit(bname, '\\.')[[1]][2]
wrt.path <- file.path(dirname(sat.img[[1]]),
paste0(f, '-reduced'))
}
utils::zip(wrt.path,
files = file.path(reduced.dir,
files.to.zip),
flags = '-q') ## might need to adjust this to different OS
## Deleting temporary folders
folders.to.delete <- dir(temporary.dir, pattern = '[.]SAFE')
unlink(paste0(normalizePath(temporary.dir), "/", folders.to.delete),
recursive = TRUE)
unlink(reduced.dir,
recursive = TRUE)
}
}
#' Retrieves the eval script
#' @name .pa_get_eval_script
#' @description
#' This function is used when building the json body to retrieve the javascript
#' file to be included in the json body
#' @param vegetation.index the vegetation index corresponding to the script to be used
#' @noRd
.pa_get_eval_script <- function(collection, vegetation.index) {
js.dir <- system.file(paste0("js/", collection), package = "pacu")
vegetation.index <- paste0(tolower(vegetation.index), '.js')
index.src <- file.path(js.dir, vegetation.index)
es <- paste(readLines(index.src), collapse = '\n')
return(es)
}
#' Builds the request body for the statistics API
#' @name .build_request_body
#' @description
#' this function is used to build the json request body of the functions that interact
#' with the data space api
#' @param aoi area of interest for which the vegetation index will be requested
#' @param start_date beginning of the time window to request data
#' @param end_date end of the time window to request data
#' @param vegetation.index vegetation index to be calculated
#' @param agg.time time interval to request images
#' @noRd
.pa_build_request_body <- function(aoi,
start.date,
end.date,
vegetation.index,
agg.time,
collection) {
if(!inherits(aoi, 'sf')) {
stop('aoi must be of class sf')
}
if(nrow(aoi) == 1){
type <- 'Polygon'
}else {
type <- 'MultiPolygon'
}
geo_json_geometry <- list(type = jsonlite::unbox(type),
coordinates = list())
for(i in 1:nrow(aoi)){
poli <- aoi[i, ]
poli <- sf::st_geometry(poli)
outside <- sf::st_multipolygon(lapply(poli, function(x) x[1]))
inside <- try(sf::st_multipolygon(rev(lapply(poli, function(x) x[2]))),
silent = TRUE)
if (inherits(inside, 'try-error'))
inside <- NULL
poli <- c(outside) ## right now we only use the outside polygon
## this is not ideal. I will fix this later.
poli <- sf::st_geometry(poli)
st_crs(poli) <- sf::st_crs(aoi)
bbox <- sf::st_cast(sf::st_geometry(poli), 'POINT')
c1 <- sf::st_coordinates(bbox[1])
epsg.code <- .pa_coord2utm(unname(c1[1]), unname(c1[2]))
bbox <- sf::st_transform(bbox, epsg.code)
if(type == 'Polygon'){
geo_json_geometry[[2]][[i]] <- vector('list', length(bbox))
for (j in 1:length(bbox)) {
geo_json_geometry[[2]][[i]][[j]] <- c(sf::st_coordinates(bbox[j]))
}
}else{
geo_json_geometry[[2]][[i]] <- vector('list', 1)
geo_json_geometry[[2]][[i]][[1]] <- vector('list', length(bbox))
for (j in 1:length(bbox)) {
geo_json_geometry[[2]][[i]][[1]][[j]] <- c(sf::st_coordinates(bbox[j]))
}
}
}
geometry_filter <- list(
bounds = list(geometry = c(geo_json_geometry),
properties =
list(crs = jsonlite::unbox(paste0("http://www.opengis.net/def/crs/EPSG/0/", epsg.code)))))
vi.script <- jsonlite::unbox(.pa_get_eval_script(collection, vegetation.index))
aggregation_filter <- list(
timeRange = list(
from= jsonlite::unbox(paste0(start.date,"T00:00:00Z" )),
to= jsonlite::unbox(paste0(end.date,"T23:59:59Z" ))),
aggregationInterval = list(of = jsonlite::unbox(agg.time)),
evalscript = vi.script,
resx = jsonlite::unbox(10),
resy = jsonlite::unbox(10))
data_filter <- list(
data = list(
list(
type = jsonlite::unbox(collection),
dataFilter = c()
)
)
)
calculations_filter <- list(default = list(
statistics = list(default = list(percentiles = list(k = c(50)),
interpolation = jsonlite::unbox('higher')))))
filter_configs <- list(
input = c(geometry_filter, data_filter),
aggregation = aggregation_filter,
calculations = calculations_filter)
body_json <- jsonlite::toJSON(filter_configs,
pretty = TRUE,
digits = 12)
return(body_json)
}
#'
#' @title Select files in the S2A and S2B products to crop the file
#' @description Select files in the S2A and S2B products to crop the file
#' @name .pa_select_s2_files
#' @rdname .pa_select_s2_files
#' @param fpath a file path from which the S2A/S2B files will be select
#' @return a vector with relevant file names
#' @noRd
.pa_select_s2_files <- function(fpath,
which = c('all', 'images', 'metadata')) {
s.wrns <- get("suppress.warnings", envir = pacu.options)
s.msgs <- get("suppress.messages", envir = pacu.options)
which <- match.arg(which)
imgList <- utils::unzip(fpath, list = TRUE)
bname <- basename(fpath)
if (!grepl('^S2A|^S2B', x = bname))
stop('Only S2A and S2B functions are supported for now.')
metadata <- .pa_read_s2_metadata(fpath, to.raw.file = TRUE)
graticule <- basename(unlist(metadata$General_Info$Product_Info$Product_Organisation))
img.indices <- sapply(graticule, function(x) grep(x, imgList[[1]], value = TRUE))
cld.str <- 'B00\\.jp2|B00\\.gml'
cloud_mask <- grep(cld.str, imgList[[1]], ignore.case = TRUE, value = TRUE)
mtd <- grep('MTD_MSI', imgList[[1]], value = TRUE)
if(which == 'all')
rel.files <- basename(unlist(c(img.indices, cloud_mask, mtd)))
if(which == 'images')
rel.files <- basename(unlist(c(img.indices, cloud_mask)))
if(which == 'metadata')
rel.files <- basename(unlist(c(mtd)))
return(rel.files)
}
#'
#' @title Retrieves the file names of a band
#' @description Retrieves the file names of a band
#' @name .pa_get_band
#' @rdname .pa_get_band
#' @param band a string representing a band name
#' @param src.dir the directory in which to look for the band
#' @param pixel.res pixel resolution
#' @param extensions image extensions
#' @return filename of the selected band
#' @noRd
.pa_get_band <- function(band,
src.dir,
pixel.res = c('default', '10m', '20m', '60m'),
extensions = c('jp2', 'tif')) {
extensions <- paste(extensions, collapse = '|')
extensions <- paste0('\\.(', extensions, ')$')
pixel.res <- match.arg(pixel.res)
flist <- list.files(src.dir, recursive = TRUE, full.names = TRUE)
if (pixel.res == 'default') pixel.res <- NULL
bb <- paste0(band,'.{0,}', pixel.res, '.{0,}', extensions)
ans <- grep(bb, flist, value = TRUE)
if(length(ans) > 1)
ans <- ans[1]
if(length(ans) < 1)
stop('Could not find band ', band)
ans
}
#'
#' @title Reads the metadata from a S2 file
#' @description Reads the metadata from a S2 file
#' @name .pa_read_s2_metadata
#' @rdname .pa_read_s2_metadata
#' @param fpath a file path pointing to an S2 metadata file
#' @param to.raw.file whether the file path points to
#' a raw file. In which case the metadata needs to be extracted.
#' @return a list containing the file metadata
#' @noRd
.pa_read_s2_metadata <- function(fpath, to.raw.file = FALSE) {
if (to.raw.file){
flist <- utils::unzip(fpath, list = TRUE )
mtd <- grep('MTD_MSI', flist[[1]], value = TRUE)
if(length(mtd) < 1)
stop('No metadata found in the file. Cannot navigate file organization.')
mtd2 <- grep(pattern = mtd,
flist$Name,
value = TRUE)
utils::unzip(fpath,
files = mtd2,
junkpaths = TRUE,
exdir = tempdir())
to.xml <- file.path(tempdir(),basename(mtd))
}else{
to.xml <- fpath
}
mtd <- XML::xmlTreeParse(file = to.xml)
mtd <- XML::xmlToList(mtd)
return(mtd)
}
#'
#' @title Aligns the bounding box of stars objects in a list
#' @description Aligns the bounding box of stars objects in a list
#' @name .pa_align_bbox
#' @param object a list contaning stars objects
#' @return a list containing stars objects
#' @noRd
.pa_align_bbox <- function(object){
crt.crs <- sf::st_crs(object[[1]])
object <- lapply(object, function(x) stars::st_warp(x, crs = crt.crs))
exts <- lapply(object, function(x) sf::st_as_sf(sf::st_as_sfc(sf::st_bbox(x))))
ext <- Reduce(function(x, y) {sf::st_union(x, y)},
x = exts)
dx <- dy <- unname(st_res(object[[1]])[1])
new.bb <- stars::st_as_stars(sf::st_bbox(ext), dx = dx, dy = dy)
for ( i in 1:length(object)){
aligned <- stars::st_warp(object[[i]], new.bb)
aligned <- stars::st_set_dimensions(aligned, which = 'time',
value = stars::st_get_dimension_values(object[[i]], 'time'))
object[[i]] <- aligned
}
return(object)
}
#'
#' @title Consolidates duplicated dates into one single image
#' @description Consolidates duplicated dates into one single image
#' @name .pa_consolidate_dates
#' @rdname .pa_consolidate_date
#' @param object a list contaning stars objects
#' @param fun function to be applied to duplicated dates
#' @return a list containing stars objects
#' @noRd
.pa_consolidate_dates <- function(object,
fun = function(x) mean(x, na.rm = TRUE)){
s.wrns <- get("suppress.warnings", envir = pacu.options)
s.msgs <- get("suppress.messages", envir = pacu.options)
date.grp <- sapply(object, function(x) stars::st_get_dimension_values(x, 'time'))
res <-list()
for (gp in unique(date.grp)){
gp.i <- which(date.grp == gp)
if (length(gp.i) > 1){
gp.date <- as.Date(stars::st_get_dimension_values(object[[gp.i[1]]], 'time'))
if (!s.wrns)
warning('Date ', gp.date, ' has duplicate data. Consolidating into one layer.')
var.name <- names(object[[gp.i[1]]])
cons <- stars::st_apply(do.call(c, c(object[gp.i], along = 'time')),
c("x", "y"),
FUN = fun)
names(cons) <- var.name
cons <- stars::st_redimension(cons,
new_dims = c(dim(cons), 1))
cons <- stars::st_set_dimensions(cons, names = c('x', 'y', 'time'))
cons <- stars::st_set_dimensions(cons, 3, gp.date)
res[[length(res) + 1]] <- cons
}else{
res[[length(res) + 1]] <- object[[gp.i]]
}
}
return(res)
}
#'
#' @title Checks the integrity of the Dataspace
#' @description Checks the integrity of the zip files downloaded
#' from Dataspace
#' @name .pa_check_zip_integrity
#' @rdname .pa_check_zip_integrity
#' @param x a list of zip files to check
#' @return No return value, called for side effects
#' @noRd
.pa_check_zip_integrity <- function(x){
s.wrns <- get("suppress.warnings", envir = pacu.options)
x <- x[file.exists(x)]
x <- x[grepl('\\.zip', x)]
for (i in 1:length(x)){
val <- try(suppressWarnings(utils::unzip(x[i], list = TRUE)), silent = TRUE)
if (inherits(val, 'try-error') || is.null(val)){
if (!s.wrns)
warning('File ', x[i], ' is corrupted and will be removed.')
file.remove(x[i])
}
}
}
#'
#' @title Predicts cardinal dates
#' @description Predicts cardinal dates
#' @name .pa_predict_cardinal_dates
#' @rdname .pa_predict_cardinal_dates
#' @param df a data frame with columns: doy, rvalue
#' @param model string representing which model to use
#' to predict cardinal dates
#' @param prior.means a vector of length three with prior mean values
#' @param prior.vars a vector of length three with prior variance values
#' @return vector of length three with cardinal date predictions
#' @noRd
.pa_predict_cardinal_dates <- function(df,
model = c("card3", "scard3", "agauss", "harmonic"),
prior.means,
prior.vars,
verbose = FALSE) {
### The assumption is that 'data' should have:
### 1) doy
### 2) rvalue
model <- match.arg(model)
algorithm <- "prior"
## Rescaling the predictor variable
## this helps with convergence of the nls algorithm
df$doy <- df$doy / 365
### First step is to fit a nonlinear model
### First we try using 'nls', if it fails, we try using 'minpack.lm::nlsLM'
### If we still fail, we return the prior
if (model == "card3") {
fnls <- try(stats::nls(rvalue ~ nlraa::SScard3(doy, tb, to, tm),
data = df
), silent = TRUE)
algorithm <- "nls"
if (inherits(fnls, "try-error")) {
fnls <- try(minpack.lm::nlsLM(rvalue ~ nlraa::SScard3(doy, tb, to, tm),
data = df
), silent = TRUE)
algorithm <- "LM"
if (inherits(fnls, "try-error")) {
if (verbose) warning("Model fitting failed Returning the prior")
attr(prior.means, "algorithm") <- algorithm
return(prior.means)
}
}
R2 <- nlraa::R2M(fnls)
cfs <- stats::coef(fnls)
nls.vars <- summary(fnls)$coefficients[, 2]^2
}
if (model == "scard3") {
fnls <- try(stats::nls(rvalue ~ nlraa::SSscard3(doy, tb, to, tm),
data = df
), silent = TRUE)
algorithm <- "nls"
if (inherits(fnls, "try-error")) {
fnls <- try(minpack.lm::nlsLM(rvalue ~ nlraa::SSscard3(doy, tb, to, tm),
data = df
), silent = TRUE)
algorithm <- "LM"
if (inherits(fnls, "try-error")) {
if (verbose) warning("Model fitting failed. Returning the prior")
attr(prior.means, "algorithm") <- algorithm
return(prior.means)
}
}
R2 <- nlraa::R2M(fnls)
cfs <- stats::coef(fnls)
nls.vars <- summary(fnls)$coefficients[, 2]^2
# fnls.bt <- nlraa::boot_nls(fnls, data = data, verbose = FALSE)
# nls.vars <- apply(fnls.bt$t, 2, stats::var, na.rm = TRUE)
}
if (model == "agauss") {
delta.start <- prior.means[2] / 365
sigma1.start <- (prior.means[2] - prior.means[1]) / (2 * 365)
sigma2.start <- (prior.means[3] - prior.means[2]) / (2 * 365)
fnls <- try(stats::nls(rvalue ~ nlraa::SSagauss(doy, eta = 1, beta = 0, delta, sigma1, sigma2),
start = list(delta = delta.start, sigma1 = sigma1.start, sigma2 = sigma2.start),
data = df
), silent = TRUE)
algorithm <- "nls"
if (inherits(fnls, "try-error")) {
fnls <- try(minpack.lm::nlsLM(rvalue ~ nlraa::SSagauss(doy, eta = 1, beta = 0, delta, sigma1, sigma2),
start = list(delta = delta.start, sigma1 = sigma1.start, sigma2 = sigma2.start),
data = df
), silent = TRUE)
algorithm <- "LM"
if (inherits(fnls, "try-error")) {
if (verbose) warning("Model fitting failed. Returning the prior")
attr(prior.means, "algorithm") <- algorithm
return(prior.means)
}
}
R2 <- nlraa::R2M(fnls)
delta.tb <- car::deltaMethod(fnls, "-2 * sigma1 + delta")
delta.tm <- car::deltaMethod(fnls, "2 * sigma2 + delta")
delta.to <- car::deltaMethod(fnls, "delta")
cfs <- c(delta.tb[[1]], delta.to[[1]], delta.tm[[1]])
### Need to compute the nls.vars
nls.vars <- c(delta.tb[[2]], delta.to[[2]], delta.tm[[2]])^2
}
if (model == "harmonic") {
algorithm <- "lm"
fnls <- try(stats::lm(rvalue ~ I(doy) + I(sin(2 * pi * doy)) + I(cos(2 * pi * doy)) +
I(sin(4 * pi * doy)) + I(cos(4 * pi * doy)), data = df))
harm.dates <- try(.pa_predict_harmonic_dates(stats::coef(fnls), stats::vcov(fnls), scaling.factor = 1))
if (inherits(fnls, "try-error") || any(is.na(harm.dates))) {
if (verbose) warning("Model fitting failed. Returning the prior")
attr(prior.means, "algorithm") <- algorithm
return(prior.means)
}
cfs <- c(harm.dates)
nls.vars <- c(attr(harm.dates, "se"))^2
R2 <- nlraa::R2M(fnls)
}
### Compute the posterior estimate
ans <- numeric(3)
sds <- numeric(3)
## Bringing the model estimates back to "day of the year"
cfs <- cfs * 365
nls.vars <- nls.vars * (365**2)
for (i in 1:3) {
ans[i] <- (1 / prior.vars[i] * prior.means[i] + 1 / nls.vars[i] * cfs[i]) / (1 / prior.vars[i] + 1 / nls.vars[i])
tau1 <- 1 / (1 / prior.vars[i] + 1 / nls.vars[i])
sds[i] <- sqrt(tau1)
}
attr(ans, "algorithm") <- algorithm
attr(ans, "R2") <- R2$R2
attr(ans, "stds") <- sds
return(round(ans, 8))
}
#'
#' @title Extract cardinal dates from harmonic regression coefficients
#' @description Extract cardinal dates from harmonic regression coefficients
#' @name .pa_get_dates_harmonic_regression
#' @rdname .pa_get_dates_harmonic_regression
#' @param betas a vector of length 6 with the parameter values of
#' a harmonic regression fit
#' @param scaling.factor the scaling factor used in the harmonic regression. 365 if
#' the x variable was day of the year. 1 if the x variable was fraction of the year (0-1).
#' @return vector of length three with cardinal date predictions from a harmonic fit
#' @noRd
.pa_get_dates_harmonic_regression <- function(betas, scaling.factor = 365) {
if (length(betas) != 6) {
stop("betas should be of length 6")
}
harmonic_gradient <- stats::deriv(
y ~ b0 + b1 * x / scaling.factor + b2 * sin(2 * pi * x / scaling.factor) +
b3 * cos(2 * pi * x / scaling.factor) +
b4 * sin(4 * pi * x / scaling.factor) + b5 * cos(4 * pi * x / scaling.factor),
c("x"),
function.arg = c(
"x", "scaling.factor", "b0", "b1", "b2", "b3",
"b4", "b5"
),
hessian = TRUE
)
preds <- function(x, scaling.factor, b0, b1, b2, b3, b4, b5) {
b0 + b1 * x / scaling.factor + b2 * sin(2 * pi * x / scaling.factor) +
b3 * cos(2 * pi * x / scaling.factor) +
b4 * sin(4 * pi * x / scaling.factor) + b5 * cos(4 * pi * x / scaling.factor)
}
get_gradient <- function(deriv_obj, ...) {
gr <- deriv_obj(...)
gr <- attr(gr, "gradient")
gr
}
xseq <- scaling.factor * seq(0, 1, length.out = 366)
gradient <- harmonic_gradient(xseq,
scaling.factor = scaling.factor,
b0 = betas[1], b1 = betas[2], b2 = betas[3],
b3 = betas[4], b4 = betas[5], b5 = betas[6]
)
gradient <- attr(gradient, "gradient")
gradient.sign <- sign(gradient)
turning.points <- xseq[which(diff(gradient.sign) != 0)]
roots <- c()
for (i in turning.points) {
root <- try(
stats::uniroot(
f = get_gradient, interval = c(i, i + (scaling.factor / 365)),
deriv_obj = harmonic_gradient,
scaling.factor = scaling.factor,
b0 = betas[1], b1 = betas[2], b2 = betas[3],
b3 = betas[4], b4 = betas[5], b5 = betas[6]
),
silent = TRUE
)
if (!inherits(root, "try-error")) {
roots <- c(roots, root$root)
}
}
to.index <- which.max(
preds(
x = roots,
scaling.factor = scaling.factor,
b0 = betas[1], b1 = betas[2], b2 = betas[3],
b3 = betas[4], b4 = betas[5], b5 = betas[6]
)
)
to <- roots[to.index]
tb <- roots[to.index - 1]
tm <- roots[to.index + 1]
if (length(to) < 1 || length(tb) < 1 || length(tm) < 1) {
warning("Could not find all three roots for this function")
return(c(tb = NA, to = NA, tm = NA))
}
card.dates <- c(tb, to, tm)
return(card.dates)
}
#'
#' @title Extract cardinal dates from harmonic regression coefficients
#' @description Extract cardinal dates from harmonic regression coefficients
#' @name .pa_predict_harmonic_dates
#' @rdname .pa_predict_harmonic_dates
#' @param betas a vector of length 6 with the parameter values of
#' a harmonic regression fit
#' @param sigma variance-covariance matrix of the harmonic regression fit
#' @param nsim number of simulations performed to estimate the standard deviation
#' @param scaling.factor the scaling factor used in the harmonic regression. 365 if
#' the x variable was day of the year. 1 if the x variable was fraction of the year (0-1).
#' @return vector of length three with cardinal date predictions from a harmonic fit
#' @noRd
#'
.pa_predict_harmonic_dates <- function(
betas,
sigma,
nsim = 1e3,
scaling.factor = 365) {
coef.samples <- MASS::mvrnorm(n = nsim, mu = betas, Sigma = sigma)
if (!inherits(coef.samples, "matrix")) {
coef.samples <- matrix(coef.samples, ncol = 6)
}
cdates.samples <- suppressWarnings(apply(coef.samples, 1, function(x) {
.pa_get_dates_harmonic_regression(x, scaling.factor)
}))
cdates.mc <- apply(cdates.samples, 1, function(x) {
c(mean = mean(x, na.rm = TRUE), sd = stats::sd(x, na.rm = TRUE))
})
cdates <- cdates.mc[1, ]
attr(cdates, "se") <- cdates.mc[2, ]
cdates
}
## Weather ----
#' Convert the units in a met file to standard units
#' @name .pa_convert_met_to_standard
#' @description
#' This function is used in the weather reports to convert the units
#' of the met file to standard units.
#' @param weather.data an object of met class containing the data to be summarized
#' @noRd
.pa_convert_met_to_standard <- function(weather.data){
weather.data$maxt <- with(weather.data, maxt <- maxt * 1.8 + 32)
weather.data$mint <- with(weather.data, mint <- mint * 1.8 + 32)
weather.data$rain <- with(weather.data, rain <- rain / 25.4)
weather.data
}
## Yield Monitor ----
#' Convert lat/long to UTM EPSG Code
#' @name .pa_coord2utm
#' @description
#' This function is used to convert crs from lat/long to UTM
#' @param long longitude
#' @param lat latitude
#' @noRd
.pa_coord2utm <- function(long, lat) {
utm.zone <- (floor((long + 180)/6) %% 60) + 1
epsg.code <- ifelse(lat > 0, 32600 + utm.zone, 32700 + utm.zone)
return(epsg.code)
}
#' Chops polygons to account for overlap
#' @name .pa_chop_polygons
#' @description
#' This function is used to chop overlapping polygons and evaluate the overlapped area
#' @param polygon a vector of vehicular polygons
#' @noRd
.pa_chop_polygons <- function(polygons,
use_s2 = TRUE){
if(!use_s2){
suppressMessages(sf::sf_use_s2(FALSE))
on.exit(suppressMessages(sf::sf_use_s2(TRUE)))
}
## if there are less than two polygons, we can skip the function
if (length(polygons) < 2){return(polygons)}
## process polygons in the reverse order
polygons <- rev(polygons)
pol.intersections <- sf::st_intersects(polygons, remove_self = TRUE)
int.ps <- (1:length(polygons))[lengths(pol.intersections) >= 1]
for (i in int.ps){
crt.pol <- polygons[i]
int.pols <- unlist(pol.intersections[i])
if(length(int.pols) < 1) {next}
nb.pols <- try(sf::st_union(polygons[int.pols]), silent = TRUE)
chopped.pol <- try(sf::st_difference(crt.pol, nb.pols), silent = TRUE)
if(length(chopped.pol) < 1 || inherits(chopped.pol, 'try-error')) {chopped.pol <- NULL}
polygons[i] <- chopped.pol
}
## reverse again
polygons <- rev(polygons)
polygons
}
#' Calculate a weighted mean based on area
#' @name .pa_areal_weighted_average
#' @description
#' This function will calculated weighted averages based on the area of overlap between x and y
#' based on a function that determines the intersection between polygons
#' @param x sf object containing the polygons from which the data will be extracted
#' @param y sf object containing the polygon to which the data will be averaged
#' @param fn a function that determines the overlap between polygons
#' @param cores number of cores to be used in the operation
#' @noRd
.pa_areal_weighted_average <- function(x, y, var, fn, sum = FALSE, cores = 1L){
s.wrns <- get("suppress.warnings", envir = pacu.options)
s.msgs <- get("suppress.messages", envir = pacu.options)
min.cov <- get("minimum.coverage.fraction", envir = pacu.options)
pol.intersections <- fn(y, x)
int.ps <- (1:length(y))[lengths(pol.intersections) >= 1]
y <- sf::st_geometry(y)
y <- sf::st_as_sf(y)
y[[var]] <- NA
pol.list <- lapply(pol.intersections, function(is) x[unlist(is), ])
if(cores > 1){
on.exit(parallel::stopCluster(cl))
cores.avlb <- parallel::detectCores(logical = FALSE)
ncores <- cores
if (cores > cores.avlb){
if(!s.wrns)
warning('Argument "cores" is greater than the number of available physical cores on the machine. Setting cores argument to ', cores.avlb,
immediate. = TRUE)
ncores <- cores.avlb
}
cl <- parallel::makeCluster(ncores)
parallel::clusterExport(cl, c('y', 'pol.list', 'var', 'sum', 'min.cov'), environment())
parallel::clusterEvalQ(cl, {library('sf')})
avs <- parallel::parLapply(cl,
1:length(pol.list),
function(i) {
ol.pols <- pol.list[[i]]
ol.pols <- suppressWarnings(sf::st_intersection(ol.pols, sf::st_buffer(y[i, ], 0)))
ol.pols$area <- as.numeric(sf::st_area(ol.pols))
cov.frac <- sum(ol.pols$area)/ as.numeric(sf::st_area(sf::st_buffer(y[i, ], 0)))
if (cov.frac < min.cov) { return(NULL)}
if(sum) {
wv <- as.numeric(sf::st_area(ol.pols)) * ol.pols[[var]]
wv <- sum(wv)
pol.mean <- wv / as.numeric(sf::st_area(sf::st_union(ol.pols)))
}else{
ol.pols <- as.data.frame(ol.pols)
pol.mean <- stats::weighted.mean(ol.pols[[var]],
ol.pols$area,
na.rm = TRUE)
}
data.frame(i = i, z = pol.mean)
})
}else{
avs <- lapply(1:length(pol.list),
function(i) {
ol.pols <- pol.list[[i]]
ol.pols <- suppressWarnings(sf::st_intersection(ol.pols, sf::st_buffer(y[i, ], 0)))
ol.pols$area <- as.numeric(sf::st_area(ol.pols))
cov.frac <- sum(ol.pols$area)/ as.numeric(sf::st_area(sf::st_buffer(y[i, ], 0)))
if (cov.frac < min.cov) { return(NULL)}
if(sum) {
wv <- as.numeric(sf::st_area(ol.pols)) * ol.pols[[var]]
wv <- sum(wv)
pol.mean <- wv / as.numeric(sf::st_area(sf::st_union(ol.pols)))
}else{
ol.pols <- as.data.frame(ol.pols)
pol.mean <- stats::weighted.mean(ol.pols[[var]],
ol.pols$area,
na.rm = TRUE)
}
data.frame(i = i, z = pol.mean)
})
}
avs <- do.call(rbind, avs)
y[[var]][avs$i] <- avs$z
y
}
#' Solves vehicular overlap
#' @name .solve_vehicle_overlap
#' @description
#' This function is used to solve for vehicle overlap and chop polygons
#' @param polygons a vector of vehicular polygons
#' @param cores a numerical value indicating the number of cores to be used in this operation
#' @noRd
.pa_solve_vehicle_overlap <- function(polygons, cores = 1L, verbose = FALSE) {
s.wrns <- get("suppress.warnings", envir = pacu.options)
s.msgs <- get("suppress.messages", envir = pacu.options)
buffered.polygons <- sf::st_buffer(polygons, -0.01) ## adding some tolerance
pol.intersections <- sf::st_intersects(buffered.polygons, remove_self = FALSE)
number.of.conflicts <- sum(as.numeric(lengths(pol.intersections) > 1))
if(verbose){cat('Solving polygon boundaries of', number.of.conflicts, 'overlapping polygons in', cores, 'core(s)','\n')}
pol.list <- lapply(pol.intersections, function(x) polygons[unlist(x)])
tgt.pol <- lapply(1:length(pol.intersections), function(i) which(unlist(pol.intersections[i]) == i))
if(number.of.conflicts > 0) {
if(cores != 1L){
on.exit(parallel::stopCluster(cl))
cores.avlb <- parallel::detectCores(logical = FALSE)
ncores <- cores
if (cores > cores.avlb){
if (!s.wrns)
warning('Argument "cores" is greater than the number of available physical cores on the machine. Setting cores argument to ', cores.avlb,
immediate. = TRUE)
ncores <- cores.avlb
}
cl <- parallel::makeCluster(ncores)
parallel::clusterExport(cl, c('.pa_chop_polygons', 'pol.list', 'tgt.pol'), environment())
parallel::clusterEvalQ(cl, {library('sf')})
chp.pols <- parallel::parLapply(cl,
1:length(pol.list),
function(i) {
pols <- .pa_chop_polygons(pol.list[[i]])
pols[tgt.pol[[i]]]
})
}else{
chp.pols <- lapply(1:length(pol.list),
function(i) {
pols <- .pa_chop_polygons(pol.list[[i]])
pols[tgt.pol[[i]]]
})
}
chp.pols <- do.call(c, chp.pols)
polygons <- chp.pols
}
polygons
}
#' Make vehicular polygon
#' @name .make_vehicle_polygon
#' @description
#' This function makes individual vehicular polygons based on location, swath, distance, and trajectory angle
#' @param point a point geometry
#' @param swath a swath value in meters
#' @param distance distance, in meters, between the previous point and the current
#' @param angle an angle in degrees
#' @noRd
.pa_make_vehicle_polygon <- function(point, swath, distance, angle){
## swath and distance must be in meters
## angle must be in degrees
## distance is the distance between the current point and the previous
point <- sf::st_geometry(point)
s <- swath / 2
d1 <- distance /2
a <- angle * pi / 180
x0 <- st_coordinates(point)[1]
y0 <- st_coordinates(point)[2]
topleft <- c(x0 - s, y0)
topright <- c(x0 + s, y0)
bottomleft <- c(x0 - s, y0 - 2 * d1)
bottomright <- c(x0 + s, y0 - 2 * d1)
pts <- matrix(c(bottomleft,
bottomright,
topright,
topleft,
bottomleft),
byrow = TRUE,
ncol = 2)
pol <- st_polygon(list(pts))
pol <- sf::st_sfc(sf::st_polygon(list(pts)))
rm <- matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
rpol <- (pol -point) * rm + point
sf::st_crs(rpol) <- sf::st_crs(point)
rpol
}
#'
#' @title Wrapper to fit a variogram
#' @description Wrapper to fit a variogram
#' @name .pa_fit_variogram
#' @param formula a formula object describing the trend
#' @param df an sf object containg the variables specified in the formular
#' @param robust whether to use Cressie's robust estimator
#' @param test.variogram logical, whether to test if the variogram results in valid predictions
#' @param ... additional parameters passed to the kriging function
#' @details This function will fit and return a Mattern variogram given a formula and a sf object
#' @return returns an sf object
#' @noRd
.pa_fit_variogram <- function(formula, df, robust = TRUE, fun, test.variogram = TRUE, verbose = FALSE, ...) {
if(verbose) cat('Fitting variogram \n')
if(fun == 'log') {
resp.col <- as.character(formula[[2]])
df[[resp.col]] <- log(df[[resp.col]])
}
v1 <- gstat::variogram(object = formula,
df,
cressie = robust)
variogram.list <- lapply(gstat::vgm()$short[2:9],
function(vm){
try(suppressWarnings(
suppressMessages(
gstat::fit.variogram(v1,
gstat::vgm(NA, vm, NA, NA),
fit.sills = TRUE,
fit.ranges = TRUE,
fit.kappa = TRUE))) , silent = TRUE)
})
fail.index <- which(sapply(variogram.list, function(x) inherits(x, 'try-error')))
if (length(fail.index) > 0){
variogram.list <- variogram.list[-fail.index]
}
if(length(variogram.list) < 1)
stop('Failed to fit a variogram using the provided formula and data.')
variogram.index <- which.min(sapply(variogram.list,function(x) attr(x, 'SSErr')))
variogram.list <- variogram.list[order(sapply(variogram.list,function(x) attr(x, 'SSErr')))]
f1 <- variogram.list[[1]]
if (test.variogram){
test.df.size <- nrow(df) %/% 20
test.df <- df[test.df.size, ]
for (i in 1:length(variogram.list)){
f1 <- variogram.list[[i]]
test.pred <- gstat::krige(formula,
df,
test.df,
f1,
debug.level = 0,
...)
if (!all(is.na(test.pred$var1.pred))){
break
}
}
}
list(f1, v1)
}
#'
#' @title Wrapper to paralellize the Krige function
#' @description Wrapper to paralellize the Krige function
#' @name .pa_predict
#' @param formula a formula object describing the trend
#' @param df an sf object containg the variables specified in the formula
#' @param new.df a new data frame for which predictions will be made
#' @param model a variogram model
#' @param cores number of cores to be used in the paralellization
#' @details This function will return krigged values
#' @return returns an sf object
#' @noRd
.pa_predict <- function(formula,
smooth.method,
df,
new.df,
cores = 1L,
fun,
verbose = FALSE,
...){
## for now, we are limiting the number of cores to 1, as kriging paralellization has
## shown to be fatal in windows...
if (cores > 1){
cores <- 1
}
if (is.null(new.df)) {
new.df <- df
}
# if (smooth.method == 'krige' && fun == 'log' && cores > 1){
# if(verbose)
# cat('Lognormal kriging is only supported using 1 core. Defaulting the smoothing process to 1 core.\n')
# #cores <- 1
# }
if(verbose) cat('Smoothing method: ', smooth.method,'\n',
'Number of features: ', nrow(new.df),'\n',
'Feature type: ', as.character(sf::st_geometry_type(new.df)[1]), '\n',
'Cores: ', cores, '\n', sep = '')
if (smooth.method == 'krige'){
model <- .pa_fit_variogram(formula, df,
robust = TRUE,
fun = fun,
test.variogram = TRUE,
verbose = verbose,
...)
vari <- model[[2]]
model <- model[[1]]
if(cores > 1){
if(verbose)
cat('Kriging...\n')
on.exit(parallel::stopCluster(cl))
new.df$coreTo <- cut(1:nrow(new.df), cores, labels = FALSE)
new.df <- split(new.df, new.df$coreTo)
cores <- min(cores, parallel::detectCores())
cl <- parallel::makeCluster(cores)
parallel::clusterExport(cl, c('formula', 'df', 'new.df', 'model', 'fun'), environment())
if (length(list(...))) {
parallel::clusterExport(cl, c('...'), environment())
}
parallel::clusterEvalQ(cl, {library('sf')})
preds <- parallel::parLapply(cl,
new.df,
function(pred.locs) {
if(fun == 'log') {
gstat::krigeTg(formula,
df,
pred.locs,
model,
lambda = 0,
...)
}else{
gstat::krige(formula,
df,
pred.locs,
model,
...)
}
})
preds <- do.call(rbind, preds)
if (fun == 'log'){
preds <- preds[c('var1TG.pred', 'var1TG.var')]
names(preds) <- c('var1.pred', 'var1.var', 'geometry')
}
preds <- list(preds, model, vari)
}else{
if(fun == 'log'){
preds <- gstat::krigeTg(formula,
df,
new.df,
model,
lambda = 0,
debug.level = -1 * verbose,
...)
preds <- preds[c('var1TG.pred', 'var1TG.var')]
names(preds) <- c('var1.pred', 'var1.var', 'geometry')
} else {
preds <- gstat::krige(formula,
df,
new.df,
model,
debug.level = -1 * verbose,
...)
}
preds <- list(preds, model, vari)
}
}else if (smooth.method == 'idw'){
preds <- gstat::idw(formula,
df,
new.df,
debug.level = -1 * verbose,
...)
preds <- list(preds)
}
preds
}
#'
#' @title Retrieve a vector given the data frame and the variable
#' @description Retrieve a vector given the data frame and the variable
#' @name .pa_get_variable
#' @param df an sf object containing the columns which to search for a variable
#' @param var a string indicating which variable to search for
#' @param verbose whether to print information the steps in this function
#' @details This function will return a vector containing the variable values and converted to the units needed for the ritas steps
#' @return a vector
#' @noRd
.pa_get_variable <- function(df, var, inf.unit, inf.col, verbose){
if (is.na(inf.col)){
inf.col <- .pa_get_variable_columns(df, var, verbose)
}
if (is.null(inf.col)){ return(NULL)}
tgt <- df[[inf.col]]
if (!is.na(inf.unit)) {
units(tgt) <- units::as_units(inf.unit)
}else {
comparison.vector <- .pa_get_comparison_vector(df, var)
tgt <- .pa_guess_units(tgt, var, comparison.vector,verbose)
}
tgt <- .pa_enforce_units(tgt, var)
tgt
}
#'
#' @title retrieve a vector with known units from the georeferenced information in the
#' data to help with unit guessing
#' @description retrieve a vector with known units from the georeferenced information in the
#' data to help with unit guessing
#' @name .pa_get_comparison_vector
#' @param df an sf object containing the columns which to search for a variable
#' @param var a string indicating which variable to search for
#' @return a vector
#' @noRd
.pa_get_comparison_vector <- function(df, var){
comparison.vector <- NULL
if (var == 'distance') {
## setting this comparison vector to 10% of the
## length of the data. We can possibly set this
## to a global variable
sample.size <- nrow(df) %/% 10
comp.df <- sf::st_geometry(df)
comp.df <- utils::head(comp.df, sample.size)
comparison.vector <- sf::st_distance(comp.df[1:(sample.size-1), ],
comp.df[2:sample.size, ],
by_element = TRUE)
}
return(comparison.vector)
}
#'
#' @title Retrieve a column name given the data frame and the variable
#' @description Retrieve a column name given the data frame and the variable
#' @name .pa_get_variable_columns
#' @param df an sf object containing the columns which to search for a variable
#' @param var a string indicating which variable to search for
#' @param verbose whether to print information the steps in this function
#' @details This function will return a vector containing the variable values and converted to the units needed for the ritas steps
#' @return a vector
#' @noRd
.pa_get_variable_columns <- function(df, var, verbose){
s.wrns <- get("suppress.warnings", envir = pacu.options)
s.msgs <- get("suppress.messages", envir = pacu.options)
verbose <- min(!s.wrns, verbose)
tgt.col <- .pa_get_variable_names(var)
tgt.col <- factor(tgt.col, levels = tgt.col, ordered = TRUE)
tgt.index <- which(tolower(names(df)) %in% tgt.col)
if(length(tgt.index) < 1) {
if (verbose)
warning('Unable to find the column for ', var, '\n',
immediate. = TRUE)
return(NULL)
}
tgt.names <- names(df)[tgt.index]
if(length(tgt.names) > 1){
tgt.name <- tgt.names[tolower(tgt.names) == as.character(min(tgt.col[tgt.col %in% tolower(tgt.names)]))]
if(verbose){
warning('Columns ', paste(tgt.names, collapse = ', '), ' matched the variable ', var,
'. Returning column ', tgt.name, '.\n', sep = '')
}
} else {
tgt.name <- tgt.names
}
return(tgt.name)
}
#'
#' @title Retrieve the interval between measurements from the time column
#' @description Retrieve the interval between measurements from the time column
#' @name .pa_time2interval
#' @param x a vector with the time data
#' @param fmts.to.try a vector containing strings of formats to try.
#' @details This function will convert datetime to interval
#' @return a vector
#' @noRd
.pa_time2interval <- function(x,
fmts.to.try =c("%Y-%m-%dT%H:%M:%OSZ",
"%Y-%m-%d %H:%M:%OS",
"%Y/%m/%d %H:%M:%OS",
"%m-%d-%Y %H:%M:%OS",
"%m/%d/%Y %H:%M:%OS",
"%d-%m-%Y %H:%M:%OS",
"%d/%m/%Y %H:%M:%OS",
"%Y-%m-%d %H:%M",
"%Y/%m/%d %H:%M",
"%m-%d-%Y %H:%M",
"%m/%d/%Y %H:%M",
"%d-%m-%Y %H:%M",
"%d/%m/%Y %H:%M",
"%Y-%m-%d",
"%Y/%m/%d",
"%m-%d-%Y",
"%m/%d/%Y",
"%d-%m-%Y",
"%d/%m/%Y",
"%H:%M:%OSZ",
"%H:%M:%OS",
"%I:%M:%S %p",
"%I:%M:%S%p"
)){
s.wrns <- get("suppress.warnings", envir = pacu.options)
s.msgs <- get("suppress.messages", envir = pacu.options)
if (!inherits(x, c('character', 'POSIXlt', 'POSIXt'))) {
stop('x must be either a character or a datetime object.')
}
if (is.character(x)){
x <- as.POSIXlt(x, tryFormats = fmts.to.try)
}
if(inherits(x, c('POSIXlt', 'POSIXt'))){
interval <- diff(x, units = 'secs')
interval <- c(stats::median(interval), interval)
}
interval <- as.numeric(interval)
if (length(interval[interval < 0]) > 0) {
if (!s.wrns)
warning('Negative interval values detected when estimating interval from time. Please check the time column.')
}
return(interval)
}
#'
#' @title Enforce the units necessary for the RITAS steps
#' @description Enforce the units necessary for the RITAS steps
#' @name .pa_enforce_units
#' @param vec a vector of class "units"
#' @param var a string indicating which variable the vector represents
#' @details This function will enforce the units necessary for the ritas steps, for example meters and g m-2
#' @return a vector
#' @noRd
.pa_enforce_units <- function(vec, var){
g <- s <- l <- count <- degreeN <- m <- NULL
if (var == 'flow'){
if(units::ud_are_convertible(units(vec), 'g/s')){
## flow mass
units(vec) <- units::make_units(g/s)
}else if(units::ud_are_convertible(units(vec), 'l/s')){
## flow is volume
units(vec) <- units::make_units(l/s)
}else {
## flow is count
units(vec) <- units::make_units(count/s)
}
}
if (var == 'moisture'){
units(vec) <- units::as_units(quote('%'))
}
if (var == 'angle') {
units(vec) <- units::make_units(degreeN)
}
if (var == 'interval') {
units(vec) <- units::make_units(s)
}
if (var == 'width') {
units(vec) <- units::make_units(m)
}
if (var == 'distance') {
units(vec) <- units::make_units(m)
}
vec
}
#'
#' @title Guess the units of a variable based on the values
#' @description Guess the units of a variable based on the values
#' @name .pa_guess_units
#' @param vec a vector of class \sQuote{units}
#' @param var a string indicating which variable the vector represents
#' @param comparison.vector a vector of values in metric units that can be extracted from the
#' data to help guess the units of the dataset
#' @param verbose whether to print information on this functions steps
#' @details This function will attempt to guess the units of a variable based on normal ranges.
#' I intend to improve this function in the future.
#' @return a \sQuote{units} vector
#' @noRd
.pa_guess_units <- function(vec, var, comparison.vector = NULL, verbose = FALSE) {
s.wrns <- get("suppress.warnings", envir = pacu.options)
s.msgs <- get("suppress.messages", envir = pacu.options)
vs <- summary(vec)
mu <- 'unknown'
if(var == 'angle') {
mu <- 'degreeN'
}
if(var == 'width') {
if(vs[3] > 65) mu <- 'in'
if(vs[3] < 65) mu <- 'ft' ## few combines are above 65ft
if (vs[3] < 19) mu <- 'm' ## 65 ft is close to 19m also
}
if (var == 'interval') {
mu <- 's'
}
if (var == 'flow') {
if (vs[3] > 0 && vs[3] < 100) mu <- 'lb/s'
}
if (var == 'moisture') {
mu <- '%'
}
if (var == 'mass') {
if(vs[3] >= 1000) mu <- 'g'
if(vs[3] < 2) mu <- 'kg'
if(vs[3] > 2 && vs[3] < 1000) mu <- 'lb'
}
if (var == 'distance') {
if (is.null(comparison.vector)){
if (vs[3] < 5) mu <- 'm'
if (vs[3] >= 5 && vs[3] <= 50) mu <- 'ft'
if(vs[3] > 50) mu <- 'in'
}else{
comparison.vector <- as.numeric(comparison.vector)
possible.units <- c('m', 'in', 'cm', 'ft')
different.units <- sapply(possible.units, function(unt) {
gu <- as.numeric(vs[3])
units(gu) <- units::as_units(unt)
units(gu) <- units::as_units('m')
gu
})
unit.ratios <- abs((stats::median(comparison.vector)/ different.units) - 1)
guessed.unit <- which.min(unit.ratios)
mu <- names(guessed.unit)
}
}
if(var == 'yield') {
if (vs[3] < 500 && vs[3] > 20) mu <- 'bushel/acre'
if (vs[3] <= 20) mu <- 't/ha'
if (vs[3] >= 500) mu <- 'kg/ha'
}
if(mu == 'unknown') {stop('unable to dermine the units of ', var, call. = FALSE)}
if (!s.msgs)
message('Guessing units of ', var, ' to be ', mu)
units::set_units(vec, mu, mode = 'standard')
}
#'
#' @title Search for a variable name in a dictionary containing possible names common in yield monitors
#' @description Search for a variable name in a dictionary containing possible names common in yield monitors
#' @name .pa_get_variable_names
#' @param vartype a variable type to retrieve names from the dictionary
#' @details this function will return a vector of common names for a given variable
#' @return a vector of possible names
#' @noRd
.pa_get_variable_names <- function(vartype = c('yield', 'angle', 'width',
'interval', 'flow', 'distance',
'moisture', 'time', 'mass')) {
vartype <- match.arg(vartype)
extd.dir <- system.file("extdata", package = "pacu")
var.dict <- jsonlite::read_json(file.path(extd.dir, 'variable-names.json'))
var.dict <- sapply(var.dict, function(x) unlist(x))
var.indices <- sapply(var.dict, function(x) x['column'] == vartype )
var.dict2 <- var.dict[var.indices]
var.names <- lapply(var.dict2, function(x) {x[grepl('names', names(x))]})
var.names <- unlist(unname(var.names))
var.names <- unique(var.names)
unlist(var.names)
}
#'
#' @title Estimate the angle between successive points
#' @description Estimate the angle between successive points
#' @name .pa_estimate_angle
#' @param points a vector of points
#' @details This function will estimate the angle of the trajectory between successive points
#' @return returns an sf object
#' @noRd
.pa_estimate_angle <- function(points){
if(!inherits(points, c("sf", "sfc", "sfg")))
stop("Object 'points' should be of class 'sf', 'sfc' or 'sfg'", call. = FALSE)
angles <- c()
for(i in 2:length(points)){
# pt1 <- sf::st_coordinates(points[i - 1])
# pt2 <- sf::st_coordinates(points[i])
# vec1 <- c(0, pt2[2] - pt1[2])
# vec2 <- c(pt2[1] - pt1[1], pt2[2] - pt1[2])
# dot.prod <- vec2%*%vec1
# m.det <- as.numeric(det(matrix(c(vec2, vec1), nrow = 2)))
# a <- atan2( m.det, dot.prod)
# a <- ifelse(a < 0, a + 2 * pi, a)
# a <- a * 180 / pi
# a <- as.numeric(a)
#
pt1 <- sf::st_coordinates(points[i - 1])
pt2 <- sf::st_coordinates(points[i])
a <- atan2(pt2[1] - pt1[1] , pt2[2] - pt1[2]) * 180 / pi
a
angles <- c(angles, a)
}
angles <- c(angles[1], angles)
units(angles) <- units::as_units('degreeN')
angles
}
#'
#' @title Estimate field boundary from points
#' @description Estimate field boundary from points
#' @name .pa_field_boundary
#' @param points a point geometry vector
#' @param arc arc of a circle to scan when building the concave hull around the field. Defaults to 5.
#' @param nsamples the number of polygons to be sampled for distances between consecutive and perpendicular points.
#' This is only relevant when method is 'polygons'.
#' @param method either 'concaveman' or 'polygons'.
#' @details This function will estimate the field boundary by calculating the distance between the
#' convex hull and each of the point. The closest points to the convex hull for every arc-circle are
#' selected to define the boundary.
#' @return returns an sf object
#' @noRd
.pa_field_boundary <- function(points,
arc = 5,
nsamples = 500,
buffer.multiplier = 1,
method = c('concaveman', 'polygons'),
verbose = FALSE) {
method <- match.arg(method)
if (inherits(points, 'data.frame'))
points <- sf::st_geometry(points)
crt.crs <- sf::st_crs(points)
if(is.na(crt.crs)) {
if(verbose)
cat('No CRS found. Defaulting to EPSG:4326\n')
sf::st_crs(points) <- 'epsg:4326'
}
points <- pa_2utm(points, FALSE)
if (method == 'polygons'){
dists <- .pa_distances(points, nsamples = nsamples, verbose = verbose)
df <- dists[[1]]
dp <- dists[[2]] / 2
vex <- sf::st_convex_hull(sf::st_combine(points))
grd <- sf::st_make_grid(vex, cellsize = rep((c( buffer.multiplier * dp, df)), 2))
cp <- sf::st_intersects(grd, points)
cp <- sapply(cp, function(x) length(x) > 0)
grd <- grd[cp]
conc <- sf::st_union(sf::st_buffer(grd, max(c(buffer.multiplier * dp, df)), endCapStyle = 'SQUARE', joinStyle = 'BEVEL'))
conc <- sf::st_intersection(vex, conc)
}
if (method == 'concaveman') {
if(!requireNamespace('concaveman', quietly = TRUE)){
stop('The concaveman package is required to generate boundaries automatically.',
' You can install it with "install.packages("concaveman")".')
}
pts <- sf::st_as_sf(sf::st_combine(points))
conc <- concaveman::concaveman(pts)
conc <- sf::st_zm(conc)
conc <- sf::st_geometry(conc)
}
return(conc)
}
#'
#' @title Get the perpendicular distance between combine passes
#' @description Get the perpendicular distance between combine passes
#' @name .pa_perpendicular_distances
#' @rdname .pa_perpendicular_distances
#' @param points an sf object
#' @param nsamples an sf object
#' @details This function will clean the yield monitor data based on proximity to the edge of the field and global outliers.
#' @return returns an sf object
#'
#' @noRd
#'
.pa_distances <- function(points,
nsamples = 500,
verbose = TRUE){
if(nsamples > length(points)) {nsamples <- length(points) - 1 }
smp.points <- sample(1:(length(points) - 1), nsamples)
dists.front <- sf::st_distance(points[smp.points], points[smp.points + 1], by_element = TRUE)
close.points <- sf::st_intersects(sf::st_buffer(points[smp.points], 5 * stats::median(dists.front)), points)
close.points <- lapply(close.points, unlist)
perp.distances <- c()
for (i in seq_along(smp.points)){
if (length(close.points[[i]]) < 1) {next}
crt.point <- points[smp.points[i]]
angles <- sapply(close.points[[i]], function(j) .pa_estimate_angle(c(crt.point,
points[j]))[1])
perp.point <- which(abs(abs(angles) - 90) < 5)
if(length(perp.point) > 1){
pd <- sf::st_distance(points[smp.points[i]], points[close.points[[i]][perp.point]])
perp.distances <- c(perp.distances, pd)
}
}
df <- as.numeric(stats::median(dists.front))
dp <- as.numeric(stats::median(perp.distances))
return(list(distance = df, width = dp))
}
#'
#' @title Clean yield monitor data
#' @description Clean yield monitor data
#' @name .pa_clean_yield_monitor
#' @rdname .pa_clean_yield_monitor
#' @param df an sf object
#' @param tgt.col the target column to use when cleaning the data based on standard deviations
#' @param boundary an sf object representing the field edges
#' @param edge.threshold a numerical value. Observations closer to the edge of the field than this value will be removed.
#' @param sd.filter a numerical value used to remove yield observations that are farther from the mean than this value.
#' @param verbose whether to print information on this functions steps
#' @details This function will clean the yield monitor data based on proximity to the edge of the field and global outliers.
#' @return returns an sf object
#'
#' @noRd
#'
.pa_clean_yield_monitor <- function(df,
tgt.col,
boundary,
edge.threshold = 30,
sd.filter = 3,
verbose = FALSE){
edge.distance <- NULL
initial.obs <- nrow(df)
crt.crs <- sf::st_crs(df)
if(is.na(crt.crs)) {
if(verbose)
cat('No CRS found. Defaulting to EPSG:4326\n')
sf::st_crs(df) <- 'epsg:4326'
}
df <- pa_2utm(df, FALSE)
crt.crs <- sf::st_crs(boundary)
if(is.na(crt.crs)) {
if(verbose )
cat('No CRS found. Defaulting to EPSG:4326\n')
sf::st_crs(boundary) <- 'epsg:4326'
}
boundary <- sf::st_transform(boundary, sf::st_crs(df))
if(sf::st_geometry_type(boundary) %in% c('POLYGON', 'MULTIPOLYGON')) {
boundary <- sf::st_cast(sf::st_geometry(boundary), 'MULTILINESTRING')
}
points <- suppressWarnings(sf::st_centroid(df))
df$edge.distance <- as.numeric(st_distance(points, boundary))
tgt <- df[[tgt.col]]
df <- subset(df, edge.distance > edge.threshold)
mean.tgt <- mean(tgt)
lower.limit <- 0
sd.tgt <- stats::sd(tgt)
if(inherits(mean.tgt, 'units')){
units(sd.tgt) <- units(mean.tgt)
units(lower.limit) <- units(mean.tgt)
}
upper.b <- mean.tgt + sd.filter * sd.tgt
lower.b <- mean.tgt - sd.filter * sd.tgt
df <- subset(df, df[[tgt.col]] > lower.b &
df[[tgt.col]] < upper.b &
df[[tgt.col]] >= lower.limit)
final.obs <- nrow(df)
if(verbose ) cat('A total of', initial.obs - final.obs, ' observations were removed during the cleaning process\n')
df <- df[tgt.col]
df
}
#'
#' @title Adjust for the time lag between the crop being cut by the combine and the sensor recording a data point
#' @description Adjust for the time lag between the crop being cut by the combine and the sensor recording a data point
#' @name .pa_adjust_lag
#' @param input an sf object containing the input data from a yield monitor
#' @param time.adj time, in seconds, to used to adjust the observations
#' @param time.col time columns in the input
#' @details This function will adjust for the time lag between the crop being cut by the combine and the sensor recording
#' a data point.
#' @return returns a sf object
#' @author Caio dos Santos and Fernando Miguez
#' @noRd
#'
.pa_adjust_lag <- function(input, time.adj, time.col) {
time <- NULL
if (time.adj > 0){
input$time <- cumsum(input[[time.col]])
pts <- st_geometry(input)
idf <- as.data.frame(input)
idf <- idf[-grep('geometry', names(idf))]
idf <- subset(idf, time > time.adj)
obs.shift <- dim(input)[1] - dim(idf)[1]
shifted <- pts[1:(length(pts) - obs.shift) ]
shifted <- cbind(idf, shifted)
shifted <- st_as_sf(shifted)
return(shifted)
}else{
return(input)
}
}
#'
#' @title Convert crop flow to mass
#' @description Convert crop flow to mass
#' @name .pa_flow2mass
#' @param crop.flow a vector of crop flow in lb/s
#' @param interval a vector the interval between measurements in seconds
#' @param moisture vector of the recorded grain moisture in \%
#' @details This function will convert crop flow in lb/s to dry mass in g
#' @return returns a vector containing mass, in grams
#' @author Caio dos Santos and Fernando Miguez
#' @noRd
#'
.pa_flow2mass <- function(crop.flow, interval, moisture){
## crop.flow must be in g/s
## interval in seconds
## moisture in %
mass <- crop.flow * interval
attributes(mass) <- NULL
attributes(moisture) <- NULL
mass <- .pa_moisture(mass, moisture, 0)
mass
}
#'
#' @title Calculate the length and width of rectangular experimental units
#' @description Calculate the length and width of rectangular experimental units
#' @name .pa_eu_dimentions
#' @param x an experimental unit geometry
#' @details This function will calculate the plot length and with of a rectangular geometry
#' @return returns a vector containing width and length in meters
#' @author Caio dos Santos and Fernando Miguez
#' @noRd
#'
.pa_eu_dimensions <- function(x){
area <- sf::st_area(x)
perimeter <- sf::st_perimeter(x)
l = (perimeter + sqrt((perimeter ^ 2) - (16 * area))) / 4
w = (perimeter - (2*l))/2
res <- c(w,l)
names(res) <- c('width', 'length')
return(res)
}
#'
#' @title Retrieve the kriging weights from the kriging process
#' @description Retrieve the kriging weights from the kriging process
#' @name .pa_kriging_weights
#' @rdname .pa_kriging_weights
#' @param x an sf object
#' @param formula a formula object describing the trend
#' @param new.data a new data frame for which predictions will be made
#' @param model a variogram model
#' @details This function will retrieve the kriging weights from an object. It was not written by me. I found it in one of the gstat vignettes and thought it
#' could be useful to diagnose problems.
#' @return returns an sf object
#'
#' @noRd
#'
.pa_kriging_weights <- function(x, formula, newdata, model, ...) {
weighti <- function(x, i, formula,...) {
ret <- rep(0,nrow(x))
ret[i]<-1
x[[1]]<-ret
gstat::krige(formula = formula,locations = x,...)
}
ret <- lapply(1:nrow(x), weighti, x=x, newdata=newdata[1,], model=model,formula=formula, ...)
ret <- do.call(rbind, ret)
unlist(ret[[1]])
}
#'
#' @title Retrieve the kriging weights from the kriging process
#' @description Retrieve the kriging weights from the kriging process
#' @name .available_memory
#' @details This function will retrieve the number of ram in byte.
#' @return available RAM in bytes
#'
#' @noRd
#'
.available_memory <- function()
{
## I did not write this function, I found it here:
## https://stackoverflow.com/questions/6457290/how-to-check-the-amount-of-ram
## I believe it can be useful for diagnosing code
# Get operating system
OS <- tolower(Sys.info()["sysname"])
# Branch based on OS
if(OS == "windows"){ # Windows
# System information
system_info <- system("systeminfo", intern = TRUE)
# Get available memory
value <- system_info[
grep("Available Physical Memory", system_info)
]
# Remove extraneous information
value <- gsub("Available Physical Memory: ", "", value)
value <- gsub("\\,", "", value)
# Convert to bytes
value_split <- unlist(strsplit(value, split = " "))
# Check for second value
bytes <- as.numeric(value_split[1]) * switch(
value_split[2],
"KB" = 1e03,
"MB" = 1e06,
"GB" = 1e09
)
}else if(OS == "linux"){ # Linux
# Split system information
info_split <- strsplit(system("free", intern = TRUE), split = " ")
# Remove "Mem:" and "Swap:"
info_split <- lapply(info_split, function(x){gsub("Mem:", "", x)})
info_split <- lapply(info_split, function(x){gsub("Swap:", "", x)})
# Get actual values
info_split <- lapply(info_split, function(x){x[x != ""]})
# Bind values
info_split <- do.call(rbind, info_split[1:2])
# Get free values
bytes <- as.numeric(info_split[2, info_split[1,] == "free"])
}else{ # Mac
# System information
system_info <- system("top -l 1 -s 0 | grep PhysMem", intern = TRUE)
# Get everything after comma
unused <- gsub(" .*,", "", system_info)
# Get values only
value <- gsub("PhysMem: ", "", unused)
value <- gsub(" unused.", "", value)
# Check for bytes
if(grepl("M", value)){
bytes <- as.numeric(gsub("M", "", value)) * 1e06
}else if(grepl("G", value)){
bytes <- as.numeric(gsub("G", "", value)) * 1e09
}else if(grepl("K", value)){
bytes <- as.numeric(gsub("K", "", value)) * 1e03
}
}
# Return bytes
return(bytes)
}
## Units ----
.pa_unit_system <- function(x, us, conv.fac, p = 1){
if (us == 'standard') {
if (is.null(conv.fac)){
stop('"conv.fac" argument is necessary to convert units to US standard unit system.')
}
x <- .pa_bushel_metric(x, conv.fac, to = 'standard', p)
units(x) <- (units::as_units(('bushel/acre')) ^ p)
}
if (us == 'metric') {
x <- x * (1e4 ^ p) / (1e6 ^ p)
units(x) <- (units::as_units('t/ha') ^ p)
}
x
}
.pa_bushel_metric <- function(x,
conv.fac,
to = c('standard', 'metric'),
p = 1 ) {
##conv.fac is the conversion factor in lb/bu
to <- match.arg(to)
## if we want to go to bushel, I will assume we are always coming from
## grams/m2
if (to == 'standard') {
x <- 4046^p * x / (454^p * conv.fac^p) ## x will go from g/m2 to bu/ac
}
## if we want to go to metric, I will assume we are always going to grams/m2
if(to == 'metric'){
x <- (x * conv.fac * 454) / (4046) ## x will go from bushel/ac to gram/m2
}
x
}
.pa_moisture <- function(x, crt, to,
verbose = FALSE) {
if (length(crt) == 1) {crt <- rep(crt, length(x))}
if (length(to) == 1) {to <- rep(to, length(x))}
if (verbose ) cat('Adjusting moisture to', round(to[1], 2), '%\n')
x <- x * (100 - crt) / (100 - to)
x
}
.pa_check_polygon_overlap <- function(input,
n = 500){
input <- utils::head(input, n)
crt.crs <- sf::st_crs(input)
if(is.na(crt.crs)) {
sf::st_crs(input) <- 'epsg:4326'
}
input <- pa_2utm(input, FALSE)
distance <- suppressMessages(.pa_get_variable(input, 'distance', NA, NA, FALSE))
swath <- suppressMessages(.pa_get_variable(input, 'width', NA, NA, FALSE))
angle <- suppressMessages(.pa_get_variable(input, 'angle', NA, NA, FALSE))
if(is.null(angle)) {
angle <- .pa_estimate_angle(sf::st_geometry(input))
}
swath <- units::drop_units(swath)
distance <- units::drop_units(distance)
angle <- units::drop_units(angle)
pols <- pa_make_vehicle_polygons(points = sf::st_geometry(input),
swath = swath,
distance = distance,
angle = angle)
c.pols <- pa_adjust_obs_effective_area(pols,
1:length(pols))
stats::median(c.pols$area.ratio, na.rm = TRUE)
}
## Plotting ----
.pa_plot_ts <- function(x,
time = 'time',
plot.var,
by = 'year',
xlab = 'Day of the year',
ylab = '',
main = '',
pch = 16,
legend.title = NULL,
legend.outside = FALSE,
palette = 'Set 2'){
if(inherits(x, 'stars'))
x <- as.data.frame(x)
## converting dates to doy
dates <- x[[time]]
if(any(is.na(as.Date(dates)))){
## assume doy
if (any(dates < 1) || any(dates > 365))
stop('Time is not in a recognizable date format or day of the year')
doy <- dates
year <- 1
}else{
doy <- as.numeric(strftime(dates, '%j'))
year <- as.numeric(strftime(dates, '%Y'))
}
x$year <- year
grps <- interaction(x[by])
## setting up colors and groups
number.of.groups <- length(unique(grps))
col.map <- hcl.colors(number.of.groups, palette)
## graph parms
xlims <- range(doy, na.rm = TRUE)
xlims <- xlims + c(-10, 10)
ylims <- range(x[[plot.var]], na.rm = TRUE)
ylims <- ylims * c(0.9, 1.1)
old.par <- graphics::par('mar', 'xpd')
on.exit(graphics::par(old.par))
inset <- c(0, 0)
if (legend.outside){
graphics::par(mar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)
inset <- c(-0.2, 0)
}
plot(x = 1,
xlab = xlab,
ylab = ylab,
xlim = xlims,
ylim = ylims,
main = main,
type = "n")
for (i in 1:number.of.groups){
y <- x[, plot.var]
graphics::points(doy[grps == levels(grps)[i]],
y[grps == levels(grps)[i]],
col = col.map[i],
pch = pch)
}
graphics::legend("topright",
inset = inset,
col = col.map,
legend= levels(grps),
title = legend.title,
pch= pch)
}
## Output ----
.pa_print_table <- function(x,
headers = TRUE,
frame = TRUE,
sep = ' ',
width = NULL,
digits = 3){
x1 <- x
x1.num <- suppressWarnings(apply(x1, 1:2, function(y){try(!is.na(as.numeric(y)), silent = TRUE)} ))
x1[x1.num] <- floor(as.numeric(x[x1.num]))
max.char <- apply(x1, 2, function(y) max(nchar(y)))
max.char <- ifelse(nchar(names(x)) > max.char, nchar(names(x)), max.char)
max.char <- unlist(unname(max.char))
if (!is.null(width)){
if(length(width) == 1) { width <- rep(width, length(max.char))}
if(length(width) != length(max.char)){stop('Width must be length one or same as number of columns')}
width <- ifelse(!is.na(width), width, max.char)
max.char <- width
}
frame.width <- sum(max.char) + floor(1.5 * length(max.char))
if(frame)
cat(rep('-', frame.width), sep = '', fill = TRUE)
if(headers){
hh <- names(x)
ho <- data.frame()
ho[1, 1:length(hh)] <- hh
ho[2, 1:length(hh)] <- rep('---', length(hh))
names(ho) <- names(x)
x <- rbind(ho, x)
}
for(i in 1:nrow(x)){
row.out <- x[i, ]
row.out <- unlist(row.out)
to.print <- list()
for (j in 1:length(max.char)){
entry <- row.out[j]
can.be.numeric <- suppressWarnings(try(as.numeric(entry), silent = TRUE))
if (!inherits(can.be.numeric, 'try-error') && !is.na(can.be.numeric)){
entry <- as.numeric(entry)
if (entry %% 1 == 0) {ff <- 'd'}else {ff <- 'f'}
dd <- max(max.char[j] - nchar(floor(entry) - 1), 0)
}
to.print[[j]] <- formatC(entry,
width = max.char[j] + 1,
format = ff,
digits = dd,
flag = '-')
}
cat(unlist(to.print), fill = TRUE, sep = sep)
}
if(frame)
cat(rep('-', frame.width), sep = '', fill = TRUE)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.