# R/s2random.R In spatstat/spatstat.sphere: Point patterns on the Sphere

#### Documented in s2runif

```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]
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)
}
```
spatstat/spatstat.sphere documentation built on Feb. 11, 2018, 11:33 p.m.