Skip to content

Commit

Permalink
Merge pull request #41 from USEPA/develop
Browse files Browse the repository at this point in the history
CRAN v5.5.1
  • Loading branch information
michaeldumelle authored Jan 9, 2024
2 parents 30bf38f + 3ce9190 commit a629b19
Show file tree
Hide file tree
Showing 70 changed files with 2,856 additions and 2,139 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: spsurvey
Title: Spatial Sampling Design and Analysis
Version: 5.5.0
Version: 5.5.1
Authors@R: c(
person("Michael", "Dumelle", role=c("aut","cre"),
email = "[email protected]", comment = c(ORCID = "0000-0002-3393-5529")),
Expand Down Expand Up @@ -40,5 +40,5 @@ BugReports: https://github.com/USEPA/spsurvey/issues
VignetteBuilder: knitr
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.1.2
RoxygenNote: 7.2.3
NeedsCompilation: no
10 changes: 9 additions & 1 deletion NEWS.md
Original file line number Diff line number Diff line change
@@ -1,3 +1,11 @@
# spsurvey 5.5.1

## Bug Fixes

* Fixed a bug in `revisit_dsgn()` that prevented proper printing of panels when there were multiple panels.
* Fixed a bug that prevented `grts()` and `irs()` from working properly with empty `LINESTRING` or `POLYGON` geometries.
* Fixed a bug that prevented `grts()` and `irs()` from returning coordinates when the the `geometry` column of `sframe` was not `"geometry"` and `legacy_sites` was specified ([#40](https://github.com/USEPA/spsurvey/issues/40)).

# spsurvey 5.5.0

## Minor Updates
Expand All @@ -14,7 +22,7 @@

## Bug Fixes

* Fixed a bug that caused an erorr in `grts()` and `irs()` occurred when at least
* Fixed a bug that caused an error in `grts()` and `irs()` occurred when at least
one variable name in `sframe` was named `"siteID"`, `"siteuse"`, `"replsite"`,
`"lon_WGS84"`, `"lat_WGS84"`, `"stratum"`, `"wgt"`, `"ip"`, `"caty"`, `"aux"`,
`xcoord`, `ycoord`, or `idpts` and the name of the geometry column in `sframe`
Expand Down
125 changes: 125 additions & 0 deletions R/contrast_analysis.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,125 @@
# contrast_analysis <- function(dframe, vars, siteID = NULL, weight = "weight", xcoord, ycoord, vartype = "Local",
# conf = 95, statistics = c("Mean", "Total")) {
#
#
#
# if (is.null(siteID)) {
# siteID <- "siteID"
# dframe$siteID <- paste("site", seq_len(nrow(dframe)), sep = "-")
# }
#
# # create design object
# stratum_ind <- FALSE
# stratumID <- NULL
# cluster_ind <- FALSE
# clusterID <- NULL
# weight1 <- NULL
# sizeweight <- FALSE
# sweight <- NULL
# sweight1 <- NULL
# fpcfactor_ind <- FALSE
# fpcsize <- NULL
# Ncluster <- NULL
# stage1size <- NULL
# jointprob <- NULL
#
# design <- survey_design(
# dframe, siteID, weight, stratum_ind, stratumID, cluster_ind, clusterID,
# weight1, sizeweight, sweight, sweight1, fpcfactor_ind, fpcsize, Ncluster,
# stage1size, vartype, jointprob
# )
#
# all_vars <- all.vars(vars)
# # incorporate form when all_vars has only length one (function of single variable)?
# # what about three variables? Function to assign matrix to localmean_cov?
# vars_form <- paste(all_vars, collapse = " + ")
#
# if (vartype == "Local") {
# local_weights <- localmean_weight(x = dframe[[xcoord]], y = dframe[[ycoord]], prb = 1 / dframe[[weight]])
# }
#
# if ("Mean" %in% statistics) {
# rslt <- svymean(make.formula(vars_form), design)
# rslt_con <- svycontrast(rslt, vars)
# if (vartype == "Local") {
# tw <- sum(dframe[[weight]])
# local_vars <- do.call(cbind, lapply(all_vars, function(x) get_local_vars(dframe, weight, x, rslt, "Mean")))
# cov_mx <- localmean_cov(local_vars, local_weights) / tw^2
# derivs <- lapply(all_vars, D, expr = vars)
# grad <- rbind(get_grad(all_vars, derivs, rslt))
# se_val <- sqrt(grad %*% cov_mx %*% t(grad))
# }
#
# if (vartype == "SRS") {
# se_val <- SE(rslt_con)
# }
#
# mean_df <- data.frame(Estimate = rslt_con[1], StdError = se_val, LCB95 = rslt_con[1] - 1.96 * se_val, UCB95 = rslt_con[1] + 1.96 * se_val)
# }
#
#
# if ("Total" %in% statistics) {
# rslt <- svytotal(make.formula(vars_form), design)
# rslt_con <- svycontrast(rslt, vars)
# if (vartype == "Local") {
# tw <- sum(dframe[[weight]])
# local_vars <- do.call(cbind, lapply(all_vars, function(x) get_local_vars(dframe, weight, x, rslt, "Total")))
# cov_mx <- localmean_cov(local_vars, local_weights)
# derivs <- lapply(all_vars, D, expr = vars)
# grad <- rbind(get_grad(all_vars, derivs, rslt))
# se_val <- sqrt(grad %*% cov_mx %*% t(grad))
# }
#
# if (vartype == "SRS") {
# se_val <- SE(rslt_con)
# }
#
# tot_df <- data.frame(Estimate = rslt_con[1], StdError = se_val, LCB95 = rslt_con[1] - 1.96 * se_val, UCB95 = rslt_con[1] + 1.96 * se_val)
#
# }
#
# contr_out <- list()
#
# if ("Mean" %in% statistics) {
# contr_out$Mean <- mean_df
# } else {
# contr_out$Mean <- NULL
# }
#
# if ("Total" %in% statistics) {
# contr_out$Total <- tot_df
# } else {
# contr_out$Total <- NULL
# }
#
# contr_out
#
# }
#
# get_grad <- function(all_vars, derivs, svyout) {
#
# # hoping for no name conflict with .use suffix
# grad.use.spsurvey <- rep(0, length(all_vars))
# derivs.use.spsurvey <- derivs
# svyout.use.spsurvey <- svyout
#
# for (x in all_vars) {
# assign(x, svyout.use.spsurvey[[x]])
# }
#
# for (i in seq_along(grad.use.spsurvey)) {
# grad.use.spsurvey[i] <- eval(derivs.use.spsurvey[[i]])
# }
# grad.use.spsurvey
# }
#
# get_local_vars <- function(dframe, weight, var, rslt, statistic) {
# if (statistic == "Mean") {
# val <- (dframe[[var]] - rslt[[var]]) * dframe[[weight]]
# }
#
# if (statistic == "Total") {
# val <- dframe[[var]] * dframe[[weight]]
# }
# val
# }
39 changes: 14 additions & 25 deletions R/grts.R
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@
#' geometry representing the legacy sites. spsurvey assumes that
#' the legacy sites were selected from a previous sampling design that
#' incorporated randomness into site selection and that the legacy sites
#' are elements of the current sampling frame. If \code{sframe} has a
#' are elements of the current sampling frame. If \code{sframe} has a
#' \code{POINT} or \code{MULTIPOINT} geometry, the observations in \code{legacy_sites}
#' should not also be in \code{sframe} (i.e., duplicates are not removed). Thus, \code{sframe}
#' and \code{legacy_sites} together compose the current sampling frame. If m or z values
Expand Down Expand Up @@ -153,8 +153,8 @@
#' If \code{n_over} is an unnamed, length-one vector, it's value is recycled
#' and used for each stratum. Note that if the
#' sampling design has unequal selection probabilities (\code{seltype = "unequal"}), then \code{n_over} sites
#' are given the same proportion of \code{caty_n} values as \code{n_base}.
#'
#' are given the same proportion of \code{caty_n} values as \code{n_base}.
#'
#'
#' @param n_near The number of nearest neighbor (nn) replacement sites.
#' If the sampling design is unstratified, \code{n_near} is integer from \code{1}
Expand Down Expand Up @@ -340,6 +340,7 @@ grts <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var =
legacy_caty_var = NULL, legacy_aux_var = NULL, mindis = NULL,
maxtry = 10, n_over = NULL, n_near = NULL, wgt_units = NULL,
pt_density = NULL, DesignID = "Site", SiteBegin = 1, sep = "-", projcrs_check = TRUE) {

if (inherits(sframe, c("tbl_df", "tbl"))) { # identify if tibble class elements are present
class(sframe) <- setdiff(class(sframe), c("tbl_df", "tbl"))
# remove tibble class for rownames warning
Expand All @@ -350,18 +351,17 @@ grts <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var =
# remove tibble class for rownames warning
}

geom_col_name <- attr(sframe, "sf_column")
st_geometry(sframe) <- "geometry"
if (!is.null(legacy_sites)) {
sframe_geom_name <- attr(sframe, "sf_column")
legacy_geom_name <- attr(legacy_sites, "sf_column")
names(legacy_sites)[names(legacy_sites) == legacy_geom_name] <- sframe_geom_name
st_geometry(legacy_sites) <- sframe_geom_name
st_geometry(legacy_sites) <- "geometry"
}

# save initial variable specifications for the design list later
initial_stratum_var <- stratum_var
initial_caty_var <- caty_var
initial_aux_var <- aux_var

# Create warning indicator and data frame to collect all potential issues during
# sample selection
warn_ind <- FALSE
Expand Down Expand Up @@ -436,14 +436,6 @@ grts <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var =
legacy_sites_names <- names(legacy_sites)
}

# Find geometry column name
geom_col_name <- attr(sframe, "sf_column")
if (geom_col_name != "geometry") {
# Force to geometry for other sf consistency
names(sframe)[names(sframe) == geom_col_name] <- "geometry"
st_geometry(sframe) <- "geometry"
}

## Create variables in sampling frame if needed.
# Create unique sampling frame ID values
sframe$id <- 1:nrow(sframe)
Expand Down Expand Up @@ -495,7 +487,7 @@ grts <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var =
legacy_sites$legacy <- legacy_sites[[legacy_var]]
}
}

# save initial variable specifications for the design list later
initial_legacy_stratum_var <- legacy_stratum_var
initial_legacy_caty_var <- legacy_caty_var
Expand Down Expand Up @@ -764,29 +756,26 @@ grts <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var =
if (geom_col_name != "geometry") {
# sframe prefix if necessary
if (geom_col_name %in% dsgn_names_extra) {
new_geom_col_name <- paste("sframe", geom_col_name, sep = "_")
sframe_names[sframe_names == geom_col_name] <- new_geom_col_name
geom_col_name <- new_geom_col_name
geom_col_name <- paste("sframe", geom_col_name, sep = "_")
}

# restore original column names
if (!is.null(sites_legacy)) {
names(sites_legacy)[names(sites_legacy) == "geometry"] <- geom_col_name
st_geometry(sites_legacy) <- geom_col_name
legacy_sites_names[legacy_sites_names == "geometry"] <- geom_col_name
st_geometry(legacy_sites) <- geom_col_name
}

if (!is.null(sites_base)) {
names(sites_base)[names(sites_base) == "geometry"] <- geom_col_name
sframe_names[sframe_names == "geometry"] <- geom_col_name
st_geometry(sites_base) <- geom_col_name
}

if (!is.null(sites_over)) {
names(sites_over)[names(sites_over) == "geometry"] <- geom_col_name
st_geometry(sites_over) <- geom_col_name
}

if (!is.null(sites_near)) {
names(sites_near)[names(sites_near) == "geometry"] <- geom_col_name
st_geometry(sites_near) <- geom_col_name
}
}
Expand Down
5 changes: 3 additions & 2 deletions R/grts_stratum.R
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ grts_stratum <- function(stratum, dsgn, sframe, sf_type, wgt_units = NULL, pt_de
n_size <- as.integer(ceiling(pmin(1e9, pt_density * (n_base + n_over))))
sfpts <- st_sample(sftmp, size = n_size, type = "regular", exact = TRUE)
sfpts <- st_as_sf(as.data.frame(sfpts), crs = st_crs(sftmp))
sfpts <- st_cast(sfpts, to = "POINT")
# sfpts <- st_cast(sfpts, to = "POINT")
# drop features with no points
sfpts <- sfpts[!st_is_empty(sfpts), ]
# join sites with linear features
Expand Down Expand Up @@ -145,9 +145,10 @@ grts_stratum <- function(stratum, dsgn, sframe, sf_type, wgt_units = NULL, pt_de
n_size <- as.integer(ceiling(pmin(1e9, pt_density * (n_base + n_over))))
sfpts <- st_sample(sftmp, size = n_size, type = "hexagonal", exact = TRUE)
sfpts <- st_as_sf(as.data.frame(sfpts), crs = st_crs(sftmp))
sfpts <- st_cast(sfpts, to = "POINT")
# sfpts <- st_cast(sfpts, to = "POINT")
# drop features with no points
sfpts <- sfpts[!st_is_empty(sfpts), ]
sfpts <- st_cast(sfpts, to = "POINT")
sftmp <- st_join(sfpts, sftmp)
sftmp$xcoord <- st_coordinates(sftmp)[, "X"]
sftmp$ycoord <- st_coordinates(sftmp)[, "Y"]
Expand Down
36 changes: 12 additions & 24 deletions R/irs.R
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ irs <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var = N
legacy_caty_var = NULL, legacy_aux_var = NULL, mindis = NULL,
maxtry = 10, n_over = NULL, n_near = NULL, wgt_units = NULL,
pt_density = NULL, DesignID = "Site", SiteBegin = 1, sep = "-", projcrs_check = TRUE) {

if (inherits(sframe, c("tbl_df", "tbl"))) { # identify if tibble class elements are present
class(sframe) <- setdiff(class(sframe), c("tbl_df", "tbl"))
# remove tibble class for rownames warning
Expand All @@ -48,18 +49,17 @@ irs <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var = N
# remove tibble class for rownames warning
}

geom_col_name <- attr(sframe, "sf_column")
st_geometry(sframe) <- "geometry"
if (!is.null(legacy_sites)) {
sframe_geom_name <- attr(sframe, "sf_column")
legacy_geom_name <- attr(legacy_sites, "sf_column")
names(legacy_sites)[names(legacy_sites) == legacy_geom_name] <- sframe_geom_name
st_geometry(legacy_sites) <- sframe_geom_name
st_geometry(legacy_sites) <- "geometry"
}

# save initial variable specifications for the design list later
initial_stratum_var <- stratum_var
initial_caty_var <- caty_var
initial_aux_var <- aux_var

# Create warning indicator and data frame to collect all potential issues during
# sample selection
warn_ind <- FALSE
Expand Down Expand Up @@ -132,15 +132,6 @@ irs <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var = N
legacy_sites_names <- names(legacy_sites)
}

# Find geometry column name
geom_col_name <- attr(sframe, "sf_column")
if (geom_col_name != "geometry") {
# Force to geometry for other sf consistency
names(sframe)[names(sframe) == geom_col_name] <- "geometry"
st_geometry(sframe) <- "geometry"
}


## Create variables in sample frame if needed.
# Create unique sample frame ID values
sframe$id <- 1:nrow(sframe)
Expand Down Expand Up @@ -192,7 +183,7 @@ irs <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var = N
legacy_sites$legacy <- legacy_sites[[legacy_var]]
}
}

# save initial variable specifications for the design list later
initial_legacy_stratum_var <- legacy_stratum_var
initial_legacy_caty_var <- legacy_caty_var
Expand Down Expand Up @@ -461,29 +452,26 @@ irs <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var = N
if (geom_col_name != "geometry") {
# sframe prefix if necessary
if (geom_col_name %in% dsgn_names_extra) {
new_geom_col_name <- paste("sframe", geom_col_name, sep = "_")
sframe_names[sframe_names == geom_col_name] <- new_geom_col_name
geom_col_name <- new_geom_col_name
geom_col_name <- paste("sframe", geom_col_name, sep = "_")
}

# restore original column names
if (!is.null(sites_legacy)) {
names(sites_legacy)[names(sites_legacy) == "geometry"] <- geom_col_name
st_geometry(sites_legacy) <- geom_col_name
legacy_sites_names[legacy_sites_names == "geometry"] <- geom_col_name
st_geometry(legacy_sites) <- geom_col_name
}

if (!is.null(sites_base)) {
names(sites_base)[names(sites_base) == "geometry"] <- geom_col_name
sframe_names[sframe_names == "geometry"] <- geom_col_name
st_geometry(sites_base) <- geom_col_name
}

if (!is.null(sites_over)) {
names(sites_over)[names(sites_over) == "geometry"] <- geom_col_name
st_geometry(sites_over) <- geom_col_name
}

if (!is.null(sites_near)) {
names(sites_near)[names(sites_near) == "geometry"] <- geom_col_name
st_geometry(sites_near) <- geom_col_name
}
}
Expand Down Expand Up @@ -608,7 +596,7 @@ irs <- function(sframe, n_base, stratum_var = NULL, seltype = NULL, caty_var = N
seltype = dsgn$seltype, caty_var = initial_caty_var, caty_n = dsgn$caty_n, aux_var = initial_aux_var, legacy = dsgn$legacy_option,
mindis = dsgn$mindis, n_over = dsgn$n_over, n_near = dsgn$n_near
)

if (any(dsgn$legacy)) {
dsgn <- c(dsgn, list(legacy_stratum_var = initial_legacy_stratum_var, legacy_caty_var = initial_legacy_caty_var,
legacy_aux_var = initial_legacy_aux_var))
Expand Down
Loading

0 comments on commit a629b19

Please sign in to comment.