R/spm.R

In lcgodoy/smile: Spatial Misalignment: Interpolation, Linkage, and Estimation

##' Transforming a \code{sf} into a \code{spm} object (Internal use)
##'
##' @name sf_to_spm
##'
##' @title single \code{sf} to \code{spm}
##'
##' @param sf_obj a \code{sf} object s.t. its geometries are polygons.
##' @param n_pts a \code{numeric} scalar representing the number of points to
##'     create a grid in the study region on which the polygons in \code{sf_obj}
##'     is observed. Alternatively, it can be a vector of the same length as
##'     \code{nrow(sf_obj)}. In this case, it generates the given number of
##'     points for each polygon in \code{sf_obj}.
##' @param type a \code{character} indicating the type of grid to be
##'     generated. The options are \code{c("random", "regular",
##'     "hexagonal")}. For more details, see \code{st_sample} in the \code{sf}
##'     package.
##' @param by_polygon a \code{logical} indicating whether we should generate
##'     \code{n_pts} by polygon or for the \code{n_pts} for the whole study
##'     region.
##' @param poly_ids a \code{character} vector informing the name of the variable
##'     in \code{sf_obj} that represents the polygons unique identifiers. In
##'     case this is not informed, we assume the id of the polygons are given by
##'     their row numbers.
##' @param var_ids a scalar or vector of type \code{character} indicating the
##'     (numerical) variables that are going to be analyzed.
##' @param trunc_d truncation distance for grid points. Consider using half of
##'     the maximum distance between polygons
##' @return a named \code{list} of size 6 belonging to the class
##'     \code{spm}. This list stores all the objects necessary to fit models
##'     using the \code{\link{fit_spm}}.
##'
##' @examples
##' data(liv_lsoa) # loading the LSOA data
##'
##' msoa_spm <- sf_to_spm(sf_obj = liv_msoa, n_pts = 1000,
##'                       type = "regular", by_polygon = FALSE,
##'                       poly_ids = "msoa11cd",
##'                       var_ids = "leb_est")
single_sf_to_spm <- function(sf_obj,
                             n_pts, type = "regular",
                             by_polygon = FALSE,
                             poly_ids = NULL,
                             var_ids  = NULL,
                             trunc_d  = NULL) {
  stopifnot(inherits(sf_obj, "sf"))
  stopifnot(all(grepl("POLYGON", sf::st_geometry_type(sf_obj))))
  stopifnot(! is.null(var_ids))
  stopifnot(inherits(sf_obj, "sf"))
  stopifnot(length(n_pts) == 1 | length(n_pts) == nrow(sf_obj))
  stopifnot(type %in% c("random", "regular", "hexagonal"))
  if (is.null(poly_ids) || (! poly_ids %in% names(sf_obj))) {
    sf_obj <- transform(sf_obj,
                        poly_ids = seq_len(nrow(sf_obj)))
    poly_ids <- "poly_ids"
  }
  if (by_polygon && length(n_pts) == 1) {
    out_grid <-
      lapply(sf::st_cast(sf::st_geometry(sf_obj),
                         "POLYGON"),
             function(x, .sz, .tp) {
               sf::st_sample(x    = x,
                             size = .sz, 
                             type = .tp)
             },
             .sz = n_pts,
             .tp = type)
  } else if (length(n_pts) > 1) {
    out_grid <-
      Map(function(x, .sz, .tp) {
        sf::st_sample(x    = x,
                      size = .sz,
                      type = .tp)
      },
      x = sf::st_cast(sf::st_geometry(sf_obj),
                      "POLYGON"),
      .sz = n_pts,
      .tp = type)
  } else {
    out_grid <-
      ## sf::st_sample(x = sf::st_cast(sf::st_geometry(sf_obj),
      ##                               "POLYGON"),
      sf::st_sample(x = sf::st_geometry(sf_obj),
                    size = n_pts,
                    type = type,
                    by_polygon = by_polygon)
  }

  if (by_polygon) {
    out_grid <- sf::st_set_crs(do.call("c", out_grid),
                               sf::st_crs(sf_obj))
    ## out_grid <- sf::st_set_crs(do.call("rbind", out_grid),
    ##                            sf::st_crs(sf_obj))
  }
  out_grid_pt <-
    sf::st_join(x = sf::st_sf(geometry = out_grid),
                y = sf_obj[poly_ids],
                join = sf::st_within)
  if (any(! sf_obj[[poly_ids]] %in% out_grid_pt[[poly_ids]])) {
    empty_polys <- which(! sf_obj[[poly_ids]] %in% out_grid_pt[[poly_ids]]) 
    out_grid_aux <-
      sf::st_centroid(x = sf_obj[empty_polys, poly_ids])
    out_grid_pt <- rbind(out_grid_pt, sf::st_sf(out_grid_aux))
  }
  
  out_grid_pt <- out_grid_pt[order(out_grid_pt[[poly_ids]]), ]

  out_grid_pt <- transform(out_grid_pt,
                           x = sf::st_coordinates(out_grid_pt)[, 1],
                           y = sf::st_coordinates(out_grid_pt)[, 2])

  if (is.null(trunc_d)) {
    out_dists <- dist_from_grids(out_grid_pt, poly_ids)
  } else {
    aux_dist <-
      sf::st_is_within_distance(sf::st_geometry(sf_obj),
                                dist   = trunc_d,
                                sparse = FALSE)
    aux_lo <- aux_dist[lower.tri(aux_dist, diag = TRUE)]
    out_dists <- dist_from_grids_tr(out_grid_pt, poly_ids,
                                    aux_lo, trunc_d * 1.01)
  }
  if (length(var_ids) == 1) {
    out_var <- sf_obj[[var_ids]]
  } else {
    stop("var_ids must be a single variable")
  }
  npix <- as.numeric(sf::st_area(sf_obj))
  output <- list(
      var     = out_var,
      dists   = out_dists,
      ids_var = poly_ids,
      grid    = out_grid_pt[poly_ids],
      npix    = npix,
      sf_poly = sf::st_geometry(sf_obj)
  )
  class(output) <- append(class(output), "spm")
  return(output)
}

##' @name sf_to_spm
##' @export
sf_to_spm <- single_sf_to_spm

## print.spm <- function(x, ...) {
##   cat(sprintf("\n Number of variables to be analyzed: %d \n", NCOL(x$var)))
##   cat(sprintf("\n ID variable: %ds\n", x$poly_ids))
##   if (NCOL(x$var) == 1) {
##     cat("\n Variable summary: ")
##     summary(x$var)
##   } else {
##     apply(x$var, 2, summary)
##   }
##   cat(sprintf("\n Number of polygons in the partition: %d \n", nrow(x$sf_poly)))
##   cat(sprintf("\n Number of points in the grid over: %d \n", nrow(x$grid)))
## }

## ##' @export
## summary.spm <- function(object, ...) {
##     print.spm(object, ...)
## }