R: Select N non-NA random points from a mulitlayer Raster* object

sample_layers.R

library(sp)
library(raster)

#' Extract valuers from a multilayer Raster*
#'
#' @param R a multilayer Raster* object
#' @param pts location info for points [cell, layer] or [row, col, layer]
#' @return a vector of values 
layers_extractPoints <- function(R, pts){
    
    if (!inherits(R, 'BasicRaster')) stop("Input R must be a Raster* class")
    nc <- as.integer(ncell(R))
    nl <- as.integer(nlayers(R))
    ny <- nrow(R)
    nx <- ncol(R)


    if (!(is.data.frame(pts) || is.matrix(pts))) 
        stop("pts must be data.frame or matrix")
    
    iz <- which(names(pts) %in% c("layer", "z"))[1]
    if (length(iz) == 0) stop("pts must have 'layer' or 'z' column")
    layer <- pts[,iz]
    if (is.character(layer)) layer <- match(layer, names(R))
    
    ic <- which(names(pts) == 'cell')[1]
    if (length(ic) == 0){
        icol <- which(names(pts) == 'col')[1]
        if (length(icol) == 0)
            stop("pts must have 'cell' or 'row' and 'col' columns")
        irow <- which(names(pts) == 'col')[1]
        if (length(irow) == 0)
            stop("pts must have 'cell' or 'row' and 'col' columns")
        cell <- irow + (icol-1) * nx           
    } else {
        cell <- pts[,'cell']
    }
    
    index <- cell + (layer-1) * nc
    
    as.vector(R)[index]
}

#' Select N non-NA random points from a mulitlayer Raster* object
#' 
#' @param R a multilayer Raster* object
#' @param pts xyz values for presence points [lon, lat, layer] - these 
#'  locations are avoided when sampling. Ignored if NULL.  
#' @param N the number of points to select
#' @return data.frame of sampled points with the following
#'  \itemize{
#'      \item{lon}
#'      \item{lat}
#'      \item{col}
#'      \item{row}
#'      \item{cell}
#'      \item{layer}
#'      \item{value}
#'  }
layers_randomPoints <- function(R, pts = NULL, N = 1000){
    
    if (!inherits(R, 'BasicRaster')) stop("Input R must be a Raster* class")
    if (!is.null(pts)){
        if (!(is.data.frame(pts) || is.matrix(pts))) 
            stop("pts must be data.frame or matrix")
        if (ncol(pts) < 3) stop("pts must have at least three columns")
        ix <- which(names(pts) %in% c("lon", "x"))[1]
        if (length(ix) == 0) stop("pts must have 'lon' or 'x' column")
        iy <- which(names(pts) %in% c("lat", "y"))[1]
        if (length(iy) == 0) stop("pts must have 'lat' or 'y' column")
        iz <- which(names(pts) %in% c("layer", "z"))[1]   
        if (length(iz) == 0) stop("pts must have 'layer' or 'z' column")
    }
    # get all of the data as a vector - storage is 1,2,3 across rows
    # starting form upper left
    v <- as.vector(R)
    nc <- as.integer(ncell(R))
    nl <- as.integer(nlayers(R))
    ny <- nrow(R)
    nx <- ncol(R)
   
    
    # if pts are present, then we flag those in vector 'v'
    if (!is.null(pts)){    
        if (is.character(pts[,iz])){
            z <- match(pts[,iz], names(R))
        } else {
            z <- match(pts[,iz], 1:nl)
        }
        col <- colFromX(R, pts[, ix])
        row <- rowFromY(R, pts[, iy])
        cell <- cellFromRowCol(R, row, col)
        # this is the number of layer steps 'extra'
        bump <- 0:(nl-1L) * nc
        index <- cell + bump[z]
        v[index] <- NA
    }
    # create a dummy index
    aix <- seq.int(from = 1L, to = (nc * nl))
    # determine where NAs occur
    nna <- !is.na(v)
    if (sum(nna) < N) stop("N exceeds number of available values to select")
    # sample the *indices*
    s <- sort(sample(aix[nna], N))
    # convert to cells, layers, rows, columns, and xy
    cell <- ((s-1L) %% nc) + 1L
    layer <- ((s-1L) %/% nc) + 1L
    rc <- rowColFromCell(R, cell)
    xy <- xyFromCell(R, cell)

    x <- data.frame(
        lon = xFromCol(R, rc[,'col']), 
        lat = yFromRow(R, rc[,'row']),
        row = rc[,'row'],
        col = rc[,'col'],
        cell = cell, 
        layer = layer, 
        value = v[cell + (layer-1)*nc],
        stringsAsFactors = TRUE)
    x <- x[order(x[,'layer']),]
    invisible(x)   
}


#' Contructs a multilayer raster (real georeference) and example points
#' 
#' @param nl the number of layers to construct in the raster
#' @param increment logical, if TRUE then increment the values in each successive layer
#' @return list of 
#'  \itemize{
#'      \item{R an N-layer raster}
#'      \item{pts a data frame of points where we might have data}
#'  }
get_example <- function(nl = 10, increment = FALSE){
    npts <- 5   # number of points per layer
    nx <- 61     # ncols
    ny <- 64    # nrows
    nc <- nx*ny # ncells
    ext <- new("Extent"
        , xmin = -77.1809315429675
        , xmax = -56.6239315429675
        , ymin = 34.5233097286473
        , ymax = 50.5233097286473
        )
    m  <- matrix(1:nc, ncol = nx, nrow = ny, byrow = TRUE)
    r <- raster::raster(m, ext@xmin, ext@xmax, ext@ymin, ext@ymax,
        crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
        
    P <- generate_polygon()
    r <- raster::mask(r,P)
    
    R <- raster::stack(r)
    if (increment){
        for (i in 2:nl) R <- stack(R,r + (nc*i))
    } else {
        for (i in 2:nl) R <- stack(R,r)
    }
    names(R) <- paste0("layer.", 1:nl)
    # make a collection of points - npts per layer
    pts <- as.data.frame(raster::sampleRandom(r, npts * nlayers(R), xy = TRUE))
    pts[,'layer'] <- sample(names(R), nrow(pts), replace = TRUE)
    
    list(R = R, pts = pts)
}
    
    
#' Constructs a 513 vertices polygon for masking 
#' @return SpatialPolygons object
generate_polygon <- function(){
    p <- structure(c(-75.3274315429675, -75.6307315429675, -75.6307315429675, 
    -75.3274315429675, -75.2937315429675, -75.2937315429675, -74.9904315429675, 
    -74.9567315429675, -74.9567315429675, -74.6534315429675, -74.6197315429675, 
    -74.6197315429675, -74.3164315429675, -74.2827315429675, -74.2827315429675, 
    -74.2827315429675, -73.9794315429675, -73.9457315429675, -73.9457315429675, 
    -73.6424315429675, -73.6087315429675, -73.6087315429675, -73.3054315429675, 
    -73.2717315429675, -73.2717315429675, -72.9684315429675, -72.9347315429675, 
    -72.9347315429675, -72.6314315429675, -72.5977315429675, -72.2944315429675, 
    -72.2607315429675, -72.2607315429675, -71.9574315429675, -71.9237315429675, 
    -71.9237315429675, -71.6204315429675, -71.5867315429675, -71.5867315429675, 
    -71.2834315429675, -71.2497315429675, -71.2497315429675, -70.9464315429675, 
    -70.9127315429675, -70.9127315429675, -70.6094315429675, -70.5757315429675, 
    -70.2724315429675, -70.2387315429675, -70.2387315429675, -69.9354315429675, 
    -69.9017315429675, -69.5984315429675, -69.2951315429675, -69.2614315429675, 
    -68.9244315429675, -68.6211315429675, -68.5874315429675, -68.2504315429675, 
    -67.9471315429675, -67.9134315429675, -67.5764315429675, -67.2731315429675, 
    -67.2394315429675, -66.9024315429675, -66.5991315429675, -66.5654315429675, 
    -66.2284315429675, -65.9251315429675, -65.8914315429675, -65.5544315429675, 
    -65.2511315429675, -65.2174315429675, -64.9141315429675, -64.8804315429675, 
    -64.5434315429675, -64.2401315429675, -64.2064315429675, -63.8694315429675, 
    -63.5661315429675, -63.5324315429675, -63.2291315429675, -63.1954315429675, 
    -62.8584315429675, -62.5551315429675, -62.5214315429675, -62.2181315429675, 
    -62.1844315429675, -61.8474315429675, -61.5441315429675, -61.5441315429675, 
    -61.8474315429675, -61.8811315429675, -62.1844315429675, -62.5214315429675, 
    -62.8584315429675, -62.8921315429675, -63.1954315429675, -63.5324315429675, 
    -63.8694315429675, -63.9031315429675, -63.9031315429675, -64.2064315429675, 
    -64.5434315429675, -64.5771315429675, -64.8804315429675, -65.2174315429675, 
    -65.2511315429675, -65.5544315429675, -65.8914315429675, -66.1947315429675, 
    -66.2284315429675, -66.5317315429675, -66.5317315429675, -66.5317315429675, 
    -66.2284315429675, -66.1947315429675, -65.8914315429675, -65.5544315429675, 
    -65.5207315429675, -65.5207315429675, -65.5544315429675, -65.8914315429675, 
    -66.2284315429675, -66.5654315429675, -66.5991315429675, -66.5991315429675, 
    -66.9024315429675, -66.9361315429675, -67.2394315429675, -67.5764315429675, 
    -67.9134315429675, -67.9471315429675, -68.2504315429675, -68.5874315429675, 
    -68.6211315429675, -68.9244315429675, -69.2614315429675, -69.5984315429675, 
    -69.6321315429675, -69.9354315429675, -69.9691315429675, -70.2724315429675, 
    -70.3061315429675, -70.3061315429675, -70.6094315429675, -70.6431315429675, 
    -70.6431315429675, -70.6431315429675, -70.6094315429675, -70.3061315429675, 
    -70.3061315429675, -70.3061315429675, -70.6094315429675, -70.9464315429675, 
    -71.2834315429675, -71.3171315429675, -71.6204315429675, -71.9574315429675, 
    -72.2944315429675, -72.6314315429675, -72.6651315429675, -72.6651315429675, 
    -72.9684315429675, -73.0021315429675, -73.3054315429675, -73.6424315429675, 
    -73.9794315429675, -74.0131315429675, -74.0131315429675, -74.3164315429675, 
    -74.6534315429675, -74.9567315429675, -74.9904315429675, -75.3274315429675, 
    40.1733097286473, 40.3983097286473, 40.6483097286473, 40.8733097286473, 
    40.8983097286473, 41.1483097286473, 41.3733097286473, 41.3983097286473, 
    41.6483097286473, 41.8733097286473, 41.8983097286473, 42.1483097286473, 
    42.3733097286473, 42.3983097286473, 42.6483097286473, 42.8983097286473, 
    43.1233097286473, 43.1483097286473, 43.3983097286473, 43.6233097286473, 
    43.6483097286473, 43.8983097286473, 44.1233097286473, 44.1483097286473, 
    44.3983097286473, 44.6233097286473, 44.6483097286473, 44.8983097286473, 
    45.1233097286473, 45.1483097286473, 45.3733097286473, 45.3983097286473, 
    45.6483097286473, 45.8733097286473, 45.8983097286473, 46.1483097286473, 
    46.3733097286473, 46.3983097286473, 46.6483097286473, 46.8733097286473, 
    46.8983097286473, 47.1483097286473, 47.3733097286473, 47.3983097286473, 
    47.6483097286473, 47.8733097286473, 47.8983097286473, 48.1233097286473, 
    48.1483097286473, 48.3983097286473, 48.6233097286473, 48.6483097286473, 
    48.8733097286473, 48.6483097286473, 48.6233097286473, 48.6233097286473, 
    48.3983097286473, 48.3733097286473, 48.3733097286473, 48.1483097286473, 
    48.1233097286473, 48.1233097286473, 47.8983097286473, 47.8733097286473, 
    47.8733097286473, 47.6483097286473, 47.6233097286473, 47.6233097286473, 
    47.3983097286473, 47.3733097286473, 47.3733097286473, 47.1483097286473, 
    47.1233097286473, 46.8983097286473, 46.8733097286473, 46.8733097286473, 
    46.6483097286473, 46.6233097286473, 46.6233097286473, 46.3983097286473, 
    46.3733097286473, 46.1483097286473, 46.1233097286473, 46.1233097286473, 
    45.8983097286473, 45.8733097286473, 45.6483097286473, 45.6233097286473, 
    45.6233097286473, 45.3983097286473, 45.1483097286473, 44.9233097286473, 
    44.8983097286473, 44.6733097286473, 44.6733097286473, 44.6733097286473, 
    44.6483097286473, 44.4233097286473, 44.4233097286473, 44.4233097286473, 
    44.3983097286473, 44.1483097286473, 43.9233097286473, 43.9233097286473, 
    43.8983097286473, 43.6733097286473, 43.6733097286473, 43.6483097286473, 
    43.4233097286473, 43.4233097286473, 43.6483097286473, 43.6733097286473, 
    43.8983097286473, 44.1483097286473, 44.3983097286473, 44.6233097286473, 
    44.6483097286473, 44.8733097286473, 44.8733097286473, 44.8983097286473, 
    45.1483097286473, 45.1733097286473, 45.1733097286473, 45.1733097286473, 
    45.1733097286473, 45.1483097286473, 44.8983097286473, 44.6733097286473, 
    44.6483097286473, 44.4233097286473, 44.4233097286473, 44.4233097286473, 
    44.3983097286473, 44.1733097286473, 44.1733097286473, 44.1483097286473, 
    43.9233097286473, 43.9233097286473, 43.9233097286473, 43.8983097286473, 
    43.6733097286473, 43.6483097286473, 43.4233097286473, 43.3983097286473, 
    43.1483097286473, 42.9233097286473, 42.8983097286473, 42.6483097286473, 
    42.3983097286473, 42.3733097286473, 42.1483097286473, 41.8983097286473, 
    41.6483097286473, 41.4233097286473, 41.4233097286473, 41.4233097286473, 
    41.3983097286473, 41.1733097286473, 41.1733097286473, 41.1733097286473, 
    41.1733097286473, 41.1483097286473, 40.8983097286473, 40.6733097286473, 
    40.6483097286473, 40.4233097286473, 40.4233097286473, 40.4233097286473, 
    40.3983097286473, 40.1483097286473, 39.9233097286473, 39.9233097286473, 
    40.1483097286473, 40.1733097286473, 40.1733097286473), .Dim = c(175L, 
    2L))
    sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(p)), '1')))
}



#' Generate an example input dataset and sample it
#'
#' @param nl the number of layers to use in the input raster
#' @param N the number of random samples to retrieve
#' @return list of 
#'  \itemize{
#'      \item{R an N-layer raster}
#'      \item{pts a data frame of points where sampling avoids}
#'      \item{x a data frame of sampled points}
#'  }
run_example <- function(nl = 10, N = 100, ...){
    X <- get_example(nl=nl, ...)
    x <- layers_randomPoints(X$R, X$pts, N = N)
    X[['x']] <- x
    invisible(X)
}

#' Contructs a multilayer raster (simple georeference) and example points
#' 
#' @return list of 
#'  \itemize{
#'      \item{R an N-layer raster}
#'      \item{pts a data frame of points where we might have data}
#'  }
get_toy_example <- function(){
    npts <- 5   # number of points per layer
    nl <- 10     # number of layers
    nx <- 8     # ncols
    ny <- 12    # nrows
    nc <- nx*ny # ncells
    
    # make a single layer raster
    m  <- matrix(1:nc, ncol = nx, nrow = ny, byrow = TRUE)
    index <- seq.int(from= 1, to = nc*nl)
    r <- raster::raster(m, 0.5,nx+0.5, 0.5, ny+0.5)
    
    # make a polygon for masking the raster
    p <- cbind(
        c( 2, 1,  3,  7, 3, 6, 2),
        c( 3, 9, 10, 10, 6, 1, 3))
    P <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(p)), '1')))
    r <- raster::mask(r, P)
    
    # make a stack
    R <- raster::stack(r)
    for (i in 2:nl) R <- stack(R,r)
    names(R) <- paste0("layer.", 1:nl)
    # make a collection of points - npts per layer
    pts <- as.data.frame(raster::sampleRandom(r, npts * nlayers(R), xy = TRUE))
    pts[,'layer'] <- sample(names(R), nrow(pts), replace = TRUE)
    
    list(R = R, pts = pts)
}