R/s2random.R

In spatstat/spatstat.sphere: Point patterns on the Sphere

s2runifbox <- function(n, region = s2lonlatbox(),
                       giveup = 1000L, warn = TRUE, ...,
                       nsim = 1L, drop = TRUE, ex = NULL){
  if(length(n) != nsim){
    if(missing(nsim)){
      nsim <- length(n)
    } else if(length(n) == 1){
      n <- rep(n, nsim)
    } else{
      stop("Length of n must be 1 or nsim.")
    }
  }
  stopifnot(all(n >= 0))
  box <- s2lonlatbox(region)
  lat <- box$lat
  lon <- box$lon
  # Check for crossning 180 degree longitude
  if(lon[1] > lon[2]){
    # Make two copies of current region
    region_pos <- region_neg <- region
    lon_pos <- c(lon[1], 180)
    lon_neg <- c(-180, lon[2])
    region_pos$lon <- lon_pos
    region_neg$lon <- lon_neg
    n_pos <- rbinom(nsim, size = n, prob = diff(lon_pos)/(diff(lon_pos) + diff(lon_neg)))
    n_neg <- n - n_pos
    X_pos <- s2runifbox(n_pos, region_pos, giveup = giveup, warn = warn, nsim = nsim, drop = FALSE)
    X_neg <- s2runifbox(n_neg, region_neg, giveup = giveup, warn = warn, nsim = nsim, drop = FALSE)
    result <- mapply(superimpose.s2pp,
                     X_pos, X_neg,
                     MoreArgs = list(region = region, check = FALSE),
                     SIMPLIFY = FALSE)
  } else{
    # Now we have a single rectangle not crossing 180 degree line
    result <- vector(mode = "list", length = nsim)
    for(isim in 1:nsim) {
      lon_i <- runif(n[isim], min = lon[1], max = lon[2])
      lat_i <- 90 - acos(runif(n[isim], sin(lat[1]*pi/180), sin(lat[2]*pi/180)))*180/pi
      result[[isim]] <- s2pp(data.frame(lon = lon_i, lat = lat_i), region = region, check = FALSE)
    }
  }
  names(result) <- paste("Simulation", 1:nsim)
  if(nsim == 1L && drop)
    result <- result[[1]]
  return(result)
}

#' Generate uniform random points on sphere
#'
#' Generate uniform random points on sphere.
#'
#' @param n Number of points. Can be a vector in which case `nsim` should be
#'   equal to `length(n)` (happens automatically if `nsim` is missing).
#' @param region Region (of class `"s2region"`) to simulate the pattern in.
#'   Defaults to the unit sphere.
#' @param giveup Number of attempts in the rejection method after which the
#'   algorithm should stop trying to generate new points.
#' @param warn Logical. Whether to issue a warning if `n`` is very large. See Details.
#' @inheritDotParams s2
#' @param nsim Number of simulated realisations to be generated.
#' @param drop Logical. If `nsim = 1`` and `drop = TRUE` (the default), the
#'   result will be a point pattern, rather than a list containing a point pattern.
#' @param ex Optional. A point pattern to use as the example. If `ex` is given and
#'   `n` and `region` are missing, then `n` and `region` will be calculated from
#'   the point pattern `ex`.
#'
#' @details This function generates `n` independent random points, uniformly
#'   distributed in the region `region` on the sphere.
#'   The algorithm depends on the type of region, as follows:
#'
#'   - If `region` is a box in longitude and latitude coordinates (`"s2lonlatbox"`)
#'   then `n` independent random points, uniformly distributed in the box, are
#'   generated by assigning uniform random values to their lon,lat coordinates
#'   using the appropriate `acos` transform for the latitude coordinates.
#'   - If `region` is a polygon or cap, the algorithm uses the rejection method.
#'   It finds a `s2lonlatbox` enclosing the region, generates points in this box,
#'   and tests whether they fall in the desired region. It gives up when
#'   `giveup * n` tests have been performed without yielding `n` successes.
#'
#'   If `warn = TRUE`, then a warning will be issued if `n` is very large. The
#'   threshold is `[spatstat.options("huge.npoints")](spatstat::spatstat.options)`.
#'   This warning has no consequences, but it helps to trap a number of common errors.
#' @return A point pattern on the sphere (an object of class "s2pp") if
#'   `nsim = 1` and `drop = TRUE`, or a list of point patterns otherwise.
#' @export
#'
#' @examples
#' X <- s2runif(10)
#' loop <- cbind(lon = c(0, 60, 60, 0), lat = c(-40, -40, 40, 40))
#' poly <- s2polygon(loop)
#' s2runif(10, region = poly)
#' s2runif(c(10, 20), region = poly)
s2runif <- function(n, region = s2(...),
                    giveup = 1000, warn = TRUE, ...,
                    nsim = 1, drop = TRUE, ex = NULL)
{
  if(missing(n) && missing(region) && !is.null(ex)) {
    stopifnot(inherits(ex, "s2pp"))
    n <- npoints(ex)
    region <- s2region(ex)
  } else {
    if(missing(nsim) && length(n) > 1)
      nsim <- length(n)
    if(nsim > 1 && length(n) == 1)
      n <- rep(n, nsim)
    if(length(n) != nsim){
        stop("Length of n must be 1 or nsim.")
    }
    stopifnot(all(n >= 0))
  }

  # if(n == 0) {
  #   emp <- s2pp(region=region)
  #   if(nsim == 1) return(emp)
  #   result <- rep(list(emp), nsim)
  #   names(result) <- paste("Simulation", 1:nsim)
  #   return(result)
  # }

  if(warn) {
    nhuge <- spatstat.options("huge.npoints")
    if(sum(n) > nhuge) {
      whinge <- paste("Attempting to generate", sum(n), "random points")
      message(whinge)
      warning(whinge, call.=FALSE)
    }
  }

  switch(class(region)[1],
         s2 = {
           result <- s2runifbox(n, region, nsim=nsim, drop = FALSE)
           result <- lapply(result, function(x){x$domain <- region; x})
         },
         s2lonlatbox = {
           result <- s2runifbox(n, region, nsim=nsim, drop = FALSE)
         },
         ## Cap and polygon treated the same way: simulate in bounding box and check containment.
         s2cap =,
         s2polygon = {
           ## make a list of nsim point patterns
           result <- vector(mode="list", length=nsim)
           for(isim in 1:nsim) {
             ## rejection method
             ## initialise empty pattern
             X <- s2pp(region=region)
             ##
             ## rectangle in which trial points will be generated
             box <- s2lonlatbox(region)
             ##
             ntries <- 0
             repeat {
               ntries <- ntries + 1
               ## generate trial points in batches of n
               qq <- s2runifbox(n[isim], box)
               ## retain those which are inside 'win'
               qq <- qq[region]
               ## add them to result
               X <- superimpose.s2pp(X, qq, region = region, check = FALSE)
               ## if we have enough points, exit
               if(npoints(X) >= n[isim]) {
                 result[[isim]] <- X[1:n[isim]]
                 break
               } else if(ntries >= giveup) {
                 ## otherwise get bored eventually
                 stop(paste("Gave up after", giveup * n[isim], "trials,",
                            npoints(X), "points accepted"))
               }
             }
           }
         },
         stop("Unrecognised s2region type")
  )

  ## list of point patterns produced.
  if(nsim == 1 && drop)
    return(result[[1]])
  names(result) <- paste("Simulation", 1:nsim)
  return(result)
}