R/Pimage.R

In SWotherspoon/SGAT: Solar/Satellite Geolocation for Animal Tracking

##' 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)))
#}