##' Create object to store binned images.
##'
##' These functions provide a data structure tools to store binned
##' samples generated by the Metropolis samplers.
##'
##' These functions provide spatial binning of samples. A spatial
##' summary image is stored separately for each time step and may be
##' mosaiced into the entire study region. Separate summaries may be
##' combined to create a multiple-track summary.
##'
##' If \code{pimg} is supplied \code{grid} and \code{proj} are ignored
##' and binning is added to the existing \code{pimg}. If \code{pimg}
##' or is not supplied \code{grid} is used to build one with the
##' details from the fit object, and \code{proj} is ignored. If only
##' \code{proj} is supplied a grid is build using that projection and
##' the details from the fit object.
##'
##' The \code{proj} argument should be a PROJ.4 string, see
##' \code{\link[raster]{projection}} and \code{\link[sp]{CRS}}, or an
##' incomplete PROJ.4 name string. If the string consists only of the
##' projection family name then a central coordinate is calculated
##' from the samples.See \code{rgdal::projInfo("proj")$name} for
##' candidate strings, and \url{http://www.spatialreference.org} for
##' more details.
##'
##' @title Bin MCMC samples
##' @param x vector of POSIXct date-times for locations
##' @param type character, samples to bin, "primary" or "intermediate"
##' @param pimg a Pimage object to accept these samples, if not
##' supplied one will be created base on the inputs
##' @param grid object to use as a template for grid specification,
##' can be a anything accepted by \code{\link{raster}} or a Pimage
##' object (any previous binning will be reset), see Details
##' @param proj a character string of projection arguments, see
##' Details
##' @param ... arguments passed to \code{chain.bin}
##' @return \code{\link{Pimage}}
##' @importFrom raster raster xmin xmax ymin ymax res projection projectExtent
##' @importFrom sp CRS spTransform coordinates
##' @export
##' @examples
##' \dontrun{
##' ## Brownian motion tethered at each end
##' brownian.bridge <- function(n, r) {
##' x <- cumsum(rnorm(n, 0, 1))
##' x <- x - (x[1] + seq(0, 1, length=n) * (x[n] - x[1]))
##' r * x
##' }
##'
##' ## Number of days and number of obs
##' days <- 50
##' n <- 200
##'
## Make separation between obs gamma distributed
##' x <- rgamma(n, 3)
##' x <- cumsum(x)
##' x <- x/x[n]
##'
## Track is lissajous + brownian bridge
##' b.scale <- 0.6
##' r.scale <- sample(c(0.1, 2, 10.2), n, replace=TRUE,
##' prob=c(0.8, 0.18, 0.02))
##' set.seed(71)
##'
##' tms <- ISOdate(2001, 1, 1) + trunc(days * 24 * 60 * 60 *x)
##' lon <- 120 + 20 * sin(2 * pi * x) +
##' brownian.bridge(n, b.scale) + rnorm(n, 0, r.scale)
##' lat <- -40 + 10 *(sin(3 * 2 * pi * x) + cos(2 * pi * x) - 1) +
##' brownian.bridge(n, b.scale) + rnorm(n, 0, r.scale)
##' x0 <- cbind(lon, lat)
##' z0 <- trackMidpts(x0)
##' n3 <- 1500
##' fx <- list(x = list(array(NA_real_, c(length(lon), 2L, n3))),
##' z = list(array(NA_real_, c(length(lon)-1L, 2L, n3))),
##' model = list(time = tms))
##' fx$x[[1L]][,,1L] <- x0
##' fx$z[[1L]][,,1L] <- z0
##' for (i in seq(n3)[-1L]) {
##' fx$x[[1L]][,,i] <- jitter(x0, factor = 8L)
##' fx$z[[1L]][,,i] <- jitter(z0, factor = 12L)
##' }
##'
##' g <- raster(extent(x0) + 5, nrows = 350, ncols = 375, crs = "+proj=longlat +datum=WGS84")
##'px <- Pimage(fx, grid = g)
##'pz <- Pimage(fx, type = "intermediate", grid = g)
##'
##' for (i in seq(n3)[-1L]) {
##' fx$x[[1L]][,,i] <- jitter(x0, factor = 8L)
##' fx$z[[1L]][,,i] <- jitter(z0, factor = 12L)
##' }
##'
##' px2 <- Pimage(fx, pimg = px)
##' pz2 <- Pimage(fx, type = "intermediate", pimg = pz)
##'
##' px3 <- Pimage(fx, grid = g)
##' pz3 <- Pimage(fx, grid = g, type = "intermediate")
##'
##' ## first
##' px$p[[80]]
##' ## this should be the sum of first and last
##' px2$p[[80]]
##' ## last
##' px3$p[[80]]
##' }
Pimage <- function(x, ...) UseMethod("Pimage")
##' @rdname Pimage
##' @method Pimage POSIXct
##' @export
Pimage.POSIXct <- function(x, type = c("primary", "intermediate"),
pimg = NULL,
grid = NULL,
proj = NULL, ...) {
## process arguments
##if (!inherits(x, "POSIXct")) stop("input times must be POSIXct")
type <- match.arg(type)
## guaranteed to be a raster, of at least the world
if (inherits(grid, "Pimage")) grid <- grid[]
if (is.null(grid)) {
grid <- raster()
} else {
grid <- raster(grid)
}
dims <- rev(dim(grid)[1L:2L])
resolution <- res(grid)
## offset to cell centres
xmn <- xmin(grid) + resolution[1L]/2
xmx <- xmax(grid) - resolution[1L]/2
ymn <- ymin(grid) + resolution[2L]/2
ymx <- ymax(grid) - resolution[2L]/2
##
pim <- .Pimage(x, c(xmn, xmx, dims[1L]), c(ymn, ymx, dims[2L]), type = type,
itersbin = 0,
projection = projection(grid)
)
pim
}
.pimg <- function (tbnd) {
list(offset = c(1L, 1L), image = NULL, tbound = tbnd)
}
## raw constructor of the empty object
.Pimage <- function(x, xbound, ybound, type = c("primary", "intermediate"),...) {
n <- length(x)
type <- match.arg(type)
## by default, we populate a primary bin with missing end times
endtime <- rep(as.numeric(NA), n)
x <- unclass(x)
if (type == "intermediate") {
## drop the last
n <- n - 1
endtime <- x[-1L]
x <- x[-length(x)]
}
## raw list of internal pimg objects
p0 <- vector("list", n)
for (i in seq_along(p0)) p0[[i]] <- .pimg(c(x[i], endtime[i]))
structure(list(p = p0,
xbound = xbound,
ybound = ybound,
tbound = range(x), ...),
class = "Pimage")
}
##' @rdname Pimage
##' @method Pimage default
##' @export
Pimage.default <- function(x, type = c("primary", "intermediate"),
pimg = NULL, grid = NULL, proj = NULL, ...) {
## this better be a SGAT fit object
tstT <- inherits(x$model$time, "POSIXct")
tstX <- (inherits(x$x, "list") &&
all(sapply(x$x,function(x) inherits(x, "array"))))
tstZ <- (inherits(x$z, "list") &&
all(sapply(x$z,function(z) inherits(z, "array"))))
if (!tstT | !tstX) stop("this does not appear to an SGAT fit object")
modeltimes <- x$model$time
type <- match.arg(type)
if(type == "intermediate") {
if (!tstZ) stop("type intermediate specified no z array list found")
chain <- chainCollapse(x$z)
} else {
if (!tstX) stop("type primary specified but no x array list found")
chain <- chainCollapse(x$x)
}
## set up tests for map projection, and rgdal availibility
## options
## 1) pimg is supplied (assume it's lonlat)
## 2) no pimg, grid is supplied
## 3) no pimg, no grid, proj dictates the Pimage space along with chain extents
## 4) no nothing, build default longlat grid from chain extents
projected <- FALSE
## 3 no pimg, no grid, proj used to build default
if (is.null(pimg) & is.null(grid) & !is.null(proj)) {
rextent <- .chaingrid(chain)
if (!isLonLat(proj)) {
.check_rgdal()
proj1 <- .process_proj(proj, rextent)
if (is.null(proj1)) stop(sprintf("cannot process PROJ.4 from \"%s\"", proj))
proj <- proj1
projected <- TRUE ## check for both longlat in rextent and proj, otherwise check rgdal avail and
grid <- projectExtent(rextent, proj)
} else {
grid <- rextent
}
}
## 4 no pimg, no grid, no proj
if (is.null(pimg) & is.null(grid) & is.null(proj)) {
grid <- .chaingrid(chain)
}
## 2 no pimg, grid is supplied, proj is ignored
if (is.null(pimg) & !is.null(grid)) {
pimg <- Pimage(modeltimes, grid = grid, type = type)
if (!isLonLat(raster(grid))) {
## I think isLonLat has changed its behaviour (2014-11-25) need to check the case
## where projection is NA, but couldBeLonLat should be TRUE
.check_rgdal()
proj <- projection(grid)
projected <- TRUE
}
}
##output
## projected(flag), proj(string), pimg(Pimage),
## 1 pimg is supplied, proj and grid ignored
## need checks to see that pimg matches chain ...
###############################################
## we have to transform
if (projected) {
.check_rgdal()
tocrs <- CRS(proj)
fromcrs <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
p <- coordinates(spTransform(SpatialPoints(cbind(as.vector(chain[,1,]),as.vector(chain[,2,])), fromcrs), tocrs))
chain[,1,] <- p[,1]
chain[,2,] <- p[,2]
}
check <- dim(chain)[1] == length(pimg)
if (!check) stop(sprintf("dimensions of chain, nrow = %i do not match the length of Pimage, length = %i", dim(chain)[1], length(pimg)))
if (type == "intermediate") {
weights <- .durations(pimg)
##units(weights) <- "hours" ## could allow user control here
##weights <- as.numeric(weights) ## R controls the units, and now we drop the class
} else {
weights <- NULL
}
for (k in seq_along(pimg)) {
binx <- t(chain[k, 1:2, ])
pimg[[k]] <- chain.bin(pimg[[k]], binx, weight = weights[k], previters = pimg$itersbin, ...)
##pimg[[k]]
}
pimg$itersbin <- pimg$itersbin + nrow(binx)
pimg
}
.chaingrid <- function(x) {
xrange <- range(x[,1,])
yrange <- range(x[,2,])
## buffer this out to be sure
xrange <- xrange + c(-1, 1) * (diff(xrange) / 10)
yrange <- yrange + c(-1, 1) * (diff(yrange) / 10)
## need to determine lon/lat aspect and modify default dims
raster(nrows = 300,
ncols = 300,
xmn = xrange[1L],
xmx = xrange[2L],
ymn = yrange[1L],
ymx = yrange[2L],
crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
)
}
## Try to load rgdal and fail otherwise.
.check_rgdal <- function() {
if (!("rgdal" %in% loadedNamespaces())) {
ns <- try(loadNamespace("rgdal"))
if (isNamespace(ns)) {
message("[loaded the rgdal namespace]")
} else {
msg <- paste("This method requires the rgdal package",
"but is unable to load rgdal namespace",
sep=",")
stop(msg)
}
}
invisible(NULL)
}
##' @importFrom MASS kde2d bandwidth.nrd
chain.bin <- function(pimg, xy, weight = NULL, method = c("bin", "kde"),
hscale = 0.7, previters = 0) {
method <- match.arg(method)
xbnd <- pimg$xbound
ybnd <- pimg$ybound
## Bin the locations into the global image coords
i <- ceiling(xbnd[3]*(xy[,1]-xbnd[1])/(xbnd[2]-xbnd[1]))
j <- ceiling(ybnd[3]*(xy[,2]-ybnd[1])/(ybnd[2]-ybnd[1]))
## Delete those points outside the global image
keep <- (i >= 1 & i <= xbnd[3] & j >= 1 & j <= ybnd[3])
if(any(keep)) {
i <- i[keep]
j <- j[keep]
## Expand image to new size
if(is.null(pimg$p[[1L]]$image)) {
irange <- range(i)
jrange <- range(j)
off <- c(irange[1],jrange[1])
img <- matrix(0,diff(irange)+1,diff(jrange)+1)
} else {
irange0 <- pimg$p[[1]]$offset[1] + c(0,nrow(pimg$p[[1L]]$image)-1)
jrange0 <- pimg$p[[1]]$offset[2] + c(0,ncol(pimg$p[[1L]]$image)-1)
irange <- range(i,irange0)
jrange <- range(j,jrange0)
off <- c(irange[1],jrange[1])
if(all(irange==irange0) && all(jrange==jrange0)) {
## Keep original image
img <- pimg$p[[1L]]$image
} else {
## Expand image
img <- matrix(0,diff(irange)+1,diff(jrange)+1)
img[(irange0[1]-off[1]+1):(irange0[2]-off[1]+1),
(jrange0[1]-off[2]+1):(jrange0[2]-off[2]+1)] <- pimg$p[[1L]]$image
}
if (!is.null(weight) & previters > 0) img <- (img / sum(img)) * previters
}
## Add binned points to new image
## img <- img + weight * tabulate(nrow(img) * (j - off[2]) + i + (1 - off[1]), nbins = prod(dim(img)))
## save the original scaling
img1 <- img + tabulate(nrow(img) * (j - off[2]) + i + (1 - off[1]), nbins = prod(dim(img)))
if (!is.null(weight)) img1 <- (img1 / sum(img1)) * weight
pimg$p[[1L]]$image <- img1
pimg$p[[1L]]$offset <- off
if (method == "kde") {
x <- as.local.pimg(pimg)
kde <- kde2d(xy[,1], xy[,2],
h = pmax(1, c(bandwidth.nrd(xy[,1]), bandwidth.nrd(xy[,2]))) * hscale,
n = c(length(x$x), length(x$y)),
lims = c(range(x$x), range(x$y)))
pimg$p[[1L]]$image <- kde$z
}
##class(pimg) <- c("pimg", "list")
}
pimg
}
##' Extract parts of Pimage.
##'
##' Extraction with \code{"["} can be done in the usual ways by
##' numeric or integer indexes, the result is returned as a
##' \code{\link[raster]{raster}}. There is no replacement method,
##' attemping to do so is an error.
##'
##' Use \code{[[.Pimage} to extract a single time step, and
##' \code{subset.Pimage} to extract multiple elements without coercion
##' to Raster.
##' @rdname extract-Pimage
##' @param x Pimage object
##' @param i numeric or logical vector
##' @param j ignored
##' @param drop ignored
##' @param ... ignored
##' @return RasterLayer, or Pimage, see Details
##' @method [ Pimage
##' @seealso \code{\link{cut.Pimage}} for creating temporal
##' partitions.
##' @export
"[.Pimage" <- function(x, i, j, drop = TRUE, ...) {
## When i is missing, return all
n <- length(x)
if (missing(i)) i <- seq_len(n)
## When i is character, return the elements of the top level list
proj <- x$projection
if (all(class(i) == "character")) {
class(x) <- NULL
return(x[i])
}
## When i is logical,
if (all(class(i) == "logical")) {
i <- which(i)
}
x$p <- x$p[i]
val <- as.image.Pimage(x)
val$z[!val$z > 0] <- NA
raster(val, crs = proj)
}
##' subset one or more elements from a Pimage, without coercion to
##' Raster
##'
##' See \code{[.Pimage}
##' @method subset Pimage
##' @rdname subset-Pimage
##' @param x Pimage
##' @param ... processed only for \code{i} for index
##' @examples
##' subset(Pimage(Sys.time() + 1:3), c(2, 3))
##' @export
subset.Pimage <- function(x, ...) {
args <- list(...)
if (length(args) == 0L) {
} else {
i <- args[[1L]]
x$p <- x$p[i]
}
x
}
projection.Pimage <- function(x, asText = TRUE) {
x$projection
}
##' Sub element replacement is disallowed
##'
##' @param x Pimage
##' @param ... ignored
##' @param value ignored
##' @usage \method{[}{Pimage}(x, \dots) <- value
##' @method [<- Pimage
##' @rdname replace-Pimage
##' @export
"[<-.Pimage" <- function(x, ..., value) {
## in the spirit of total clarity
stop("disallowed")
}
## internal
## @param exact ignored
## @param ... ignored
## export
##' Extract a part of a Pimage object.
##'
##'
##' @title Extract a single Pimage time step.
##' @param x Pimage
##' @param i integer index
##' @param ... ignored
##' @param exact ignored
##' @method [[ Pimage
##' @rdname EExtract-Pimage
##' @return Pimage
##' @export
"[[.Pimage" <- function(x, i, ..., exact = TRUE) {
cl <- oldClass(x)
class(x) <- NULL
## note this has to be the 1-element list, perhaps other i-s should be an error
if(!is.character(i)) {
x[["p"]] <- x[["p"]][i[1L]]
} else {
return(x[[i]])
}
class(x) <- cl
x
}
##' Convert Pimage to POSIXct
##'
##' Coercion method to extract times from Pimage object. The type is
##' inferred from the internal storage.
##' @title as.POSIXct for Pimage
##' @seealso \code{estelleMetropolis} and \code{Pimage}
##' @param x Pimage object
##' @param ... ignored
##' @export
as.POSIXct.Pimage <- function(x, ...) {
switch(.type(x),
primary = .Xtimes(x),
intermediate = .Ztimes(x))
}
##' cut.Pimage
##'
##' Cut a Pimage object based on a date-time input breaks character string, and return a multi-layer raster
##'
##' @title cut.POSIXt for Pimage
##' @rdname cut.Pimage
##' @param x Pimage
##' @param breaks an interval specification, see \code{\link{cut.POSIXt}}
##' @param ... pass arguments to \code{\link{cut.POSIXt}}
##' @method cut Pimage
##' @return RasterLayer or RasterBrick, for multiple or single level cut respectively
##' @importFrom raster brick xmin xmax ymin ymax projection getValues setZ
##' @export
cut.Pimage <- function(x, breaks, ...) {
r1 <- x[1]
datetimes <- .Xtimes(x)
ct <- cut.POSIXt(datetimes, breaks = breaks, ...)
## now rebuild the output
resarr <- array(0.0, c(dim(r1)[1:2], nlevels(ct)))
for (i in seq_len(nlevels(ct))) {
## this is the solution to the as.matrix namespace problem (don't use it) MDSumner 2013-07-18
resarr[,,i] <- getValues(x[ct == levels(ct)[i]], format='matrix')
}
dates <- as.POSIXct(levels(ct), tz = "GMT")
resr <- brick(resarr, xmn=xmin(r1), xmx=xmax(r1), ymn=ymin(r1), ymx=ymax(r1), crs=projection(r1))
names(resr) <- format(dates, "%b.%d_%Y")
setZ(resr, dates, name = "datetime")
}
##' General methods for Pimage
##'
##' These are the basics for print, plot, concatenation, coercion.
##' @param x Pimage
##' @param ... arguments passed along
##' @rdname Pimage-methods
##' @aliases print
##' @method print Pimage
##' @export
print.Pimage <- function(x, ...) {
## this needs to know the x/y/time range, and possibly the sizes of all images, whether any are NULL or funny
a <- x[1L]
ext <- bbox(a) # equivalent to extent(a)
trange <- format(range(.Xtimes(x)))
type <- .type(x)
cat("class :", class(x), type, "\nlength :", length(x), "\ntime :", trange, "\n")
##cat("Time Steps :")
##str(attr(x, "times"))
cat("extent : ", ext[1, 1], ", ", ext[1, 2], ", ", ext[2, 1], ", ", ext[2, 2], " (xmin, xmax, ymin, ymax)\n", sep = "")
cat("CRS :", projection(a))
cat("\n")
invisible(NULL)
}
##' @method str Pimage
##' @rdname Pimage-methods
##' @param object Pimage
##' @importFrom utils str
##' @export
str.Pimage <- function(object, ...) {
class(object) <- NULL
str(object)
}
##' @method length Pimage
##' @rdname Pimage-methods
##' @export
length.Pimage <- function(x, ...) {
length(x$p)
}
##' @rdname Pimage-methods
##' @param recursive ignored
##' @method c Pimage
##' @export
c.Pimage <- function(..., recursive = FALSE) {
obj <- list(...)
projections <- sapply(obj, attr, "projection")
if (!length(unique(projections)) == 1L) stop("inputs have non-matching projections")
x0 <- sapply(obj$p, function(x) x$xbound[1L])
x1 <- sapply(obj$p, function(x) x$xbound[2L])
nx <- sapply(obj$p, function(x) x$xbound[3L])
y0 <- sapply(obj$p, function(x) x$ybound[1L])
y1 <- sapply(obj$p, function(x) x$ybound[2L])
ny <- sapply(obj$p, function(x) x$ybound[3L])
unq <- function(x) length(x) == length(unique(x))
if (!(all(c(unq(x0), unq(x1), unq(nx), unq(y0), unq(y1), unq(ny))))) stop("inputs have different grid bounds")
x <- do.call("c", lapply(obj, as.list))
p <- obj[[1L]][[1L]]
p$p <- x
p
}
##' @rdname Pimage-methods
##' @method as.list Pimage
##' @export
as.list.Pimage <- function(x, ...) {
## drop the class and attributes
x$p
}
##' @rdname Pimage-methods
##' @method plot Pimage
##' @export
plot.Pimage <- function(x, ...) {
plot(x[], ...)
}
## internal
## method [[<- Pimage
## rdname [.Pimage
## export
"[[<-.Pimage" <- function(x, i, value) {
## need further checks that these are equivalent Pimage objects
if (inherits(value, "Pimage") & length(value) == 1L) {
x$p[[i]] <- value$p[[1L]]
} else {
stop(sprintf("no method to [[.<- for objects of class %s", class(x)))
}
x
}
.type <- function(x) {
c("intermediate", "primary")[all(is.na(sapply(x$p, function(x) x$tbound[2L]))) + 1]
}
.Xtimes <- function(x) {
## Return the smaller time bound
.POSIXct(sapply(x$p, function(x) x$tbound[1L]))
}
.Ztimes <- function(x) {
## Return the middle of the time bound
.POSIXct(sapply(x$p, function(x) {xt <- x$tbound; xt[1L] + diff(xt)/2}))
}
.durations <- function(x) {
sapply(x$p, function(x) diff(x$tbound))
}
## QUERY - should this return matrix of zeros?
## workers used by as.image.Pimage and each other
## these are old, but they work on new Pimage
as.matrix.pimg <- function(x) {
img <- matrix(0, x$xbound[3L], x$ybound[3L])
## this used to do the addition work, but not now
##if(!is.null(pimg$image)) {off <- pimg$offset;img[off[1]:(off[1]+nrow(pimg$image)-1), off[2]:(off[2]+ncol(pimg$image)-1)] <- pimg$image}
img
}
as.image.pimg <- function(pimg) {
img <- coords.pimg(pimg)
img$z <- as.matrix.pimg(pimg)
img
}
coords.pimg <- function(pimg) {
list(x = seq(pimg$xbound[1L], pimg$xbound[2L], length = pimg$xbound[3L]),
y = seq(pimg$ybound[1L], pimg$ybound[2L], length = pimg$ybound[3L]))
}
## this is used by the kde method
as.local.pimg <- function (pimg) {
img <- coords.pimg(pimg)
img$x <- img$x[pimg$p[[1L]]$offset[1L]:(pimg$p[[1L]]$offset[1L] + nrow(pimg$p[[1L]]$image) - 1)]
img$y <- img$y[pimg$p[[1L]]$offset[2L]:(pimg$p[[1L]]$offset[2L] + ncol(pimg$p[[1L]]$image) - 1)]
img$z <- pimg$p[[1L]]$image
img
}
as.image.Pimage <- function (pimgs) {
## should have checks elsewhere for these NULLs, do they persist when no mixing?
## bad <- unlist(lapply(pimgs, function(x) is.null(x$image)))
res <- as.image.pimg(pimgs[[1]])
if (all(sapply(pimgs$p, function(x) is.null(x$image)))) return(res)
for (i in seq_len(length(pimgs))) {
## bit weird, but we need p[[1L]]
img <- pimgs$p[[i]]
## here nothing to add (that is a problem needs checking apart
## from the empty case, when no binning has been done)
if (is.null(img$image)) next;
Xpos <- img$offset[1L]
Ypos <- img$offset[2L]
Xind <- Xpos:(Xpos + dim(img$image)[1] - 1)
Yind <- Ypos:(Ypos + dim(img$image)[2] - 1)
res$z[Xind, Yind] <- res$z[Xind, Yind] + img$image
}
res
}
.process_proj <- function(x, ext) {
check <- try(CRS(x), silent = TRUE)
iscrs <- !inherits(check, "try-error")
tokens <- strsplit(x, " ")[[1]]
tokens <- tokens[nchar(tokens) > 0]
## succeeds CRS, and is more than just a projname
if (iscrs & length(tokens) > 1) return(x)
## fails CRS, but isn't just a proj name
if (!iscrs & length(tokens) > 1) return(NULL)
## what if it's a proj, but there's no leading "+"
if (!iscrs) {
## this is otherwise not robust to "lonlat", "latlon" aliases of "longlat", sp says no
if (tokens[1L] %in% c("lonlat", "latlon")) tokens[1L] <- "longlat"
x <- sprintf("+proj=%s", tokens[1])
}
check <- try(CRS(x), silent = TRUE)
if (inherits(check, "try-error")) return(NULL)
## so we get this far, lazy wants central coordinates from ext
## just the centre in long/lat
sprintf("%s +lon_0=%s +lat_0=%s +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0", x, round(mean(c(xmin(ext), xmax(ext)))),
round(mean(c(ymin(ext), ymax(ext)))))
}
## not required
##".times<-" <- function(x, value) {
## attr(x, "times") <- value
## x
##}
## functions for merging global datasets with pimage
## .aligned(p[], global)
## plot(parent)
## points(xyFromCell(parent, cn.pimg(p[[i]])), pch = ".")
## plot(extent(trim(p[i], value = 0)), add = TRUE)
## check that extents do align, otherwise hint at rebuilding Pimage
## to match
#.aligned <- function(x, y, ...) {
# all.equal(extent(x), alignExtent(extent(x), y, snap = "out"))
#}
## function to convert pimg offset to cellnumbers of PARENT
## (then we can crop global to match parent and do direct transfers)
#cn.pimg <- function(x) {
# xbnd <- x$xbound
# ybnd <- x$ybound
# offs <- x$offset
# tl <- ((ybnd[3L] - (offs[2L] + ncol(x$image))) + 1L) * xbnd[3L] + offs[1L]
# sort(rep(seq(tl, by = xbnd[3], length = ncol(x$image)), each = nrow(x$image)) +
# rep(seq_len(nrow(x$image)) - 1, ncol(x$image)))
#}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.