Nothing
#' Set or get the ordering of a grid
#'
#' Specify or retrieve the ordering in which a grid is stored in a vector (or matrix).
#'
#' @param x a data container for the elements of the grid; the grid order is stored as an attribute to it
#' @param refpoint a string specifying which point of the grid corresponds
#' to the first element of `x`; see below
#' @param cycle a permutation of the integers `1:G` (see below)
#'
#' @return `setGridOrder(x,...)` returns the object `x` with the grid order description attached
#' as an attribute "gridOrder"; `getGridOrder(x)` retrieves this attribute and returns it.
#' @details A "gridOrder" attribute is a list consisting of two named elements:
#' \describe{
#' \item{refpoint}{one of "topleft", "bottomleft", "topright" or "bottomright" in 2D, or
#' also of "topleftsurf", "bottomleftsurf", "toprightsurf", "bottomrightsurf", "topleftdeep","bottomleftdeep",
#' "toprightdeep" or "bottomrightdeep" in 3D ("deep" is accessory, i.e. "topleft"== "topleftdeep"),
#' indicating the location on the grid of the first point of the object `x`}
#' \item{cycle}{a permutation of `1:G` indicating in which order run the dymensions, from faster to slower}
#' }
#'
#' Thus, a conventional ordering of a (nX*nY)-element vector into a matrix to plot with [graphics::image()]
#' corresponds to an `refpoint="bottomleft"` and `cycle=1:2`, i.e. start with the lower left corner
#' and run first by rows (eastwards), then by columns (northwards). This is constructed by
#' `gridOrder_array(G)`, and can be directly set to an object `x` by `setGridOrder_array(x,G)`;
#' `gridOrder_GSLib` is an alias for `gridOrder_array`.
#'
#' The grids from package "sp" (and many other in R), on the contrary follow the convention
#' `refpoint="topleft"` and `cycle=1:2`, i.e. start with the upper left corner
#' and run first by rows (eastwards), then by columns (**southwards**). This is constructed by
#' `gridOrder_sp(G)`, and can be directly set to an object `x` by `setGridOrder_sp(x,G)`; `gridOrder_gstat`
#' is an alias for `gridOrder_sp`.
#'
#' @export
#' @aliases getGridOrder
#' @seealso [sortDataInGrid()] for ways of reordering a grid
#' @examples
#' gt = sp::GridTopology(cellcentre.offset=c(1,11), cellsize=c(1,1), cells.dim=c(5,3))
#' sp::coordinates(sp::SpatialGrid(grid=gt))
#' gridOrder_sp(2)
setGridOrder = function(x, refpoint, cycle){
refpoints = gmApply(expand.grid(c("top","bottom"), c("left","right"), c("","surf", "deep")),1, paste, collapse="")
if(!(refpoint %in% refpoints)) stop(paste("setGridOrder: refpoint not understood, must be one of ", paste(refpoints, collapse=", ")))
if(any(sort(cycle)!=(1:length(cycle)))) stop("setGridOder: cycle must be a permutation of c(1,2) in 2D, or c(1,2,3) in 3D")
attr(x,"gridOrder") = list(refpoint=refpoint, cycle=cycle)
return(x)
}
#' @export
getGridOrder = function(x){
res = attr(x, "gridOrder")
if(is.null(res)){
warning("getGridOrder: no grid ordering info found, taking the default, from options()$gmGeostats$gridOrder")
res = options()$gmGeostats$gridOrder
}
return(res)
}
#' @describeIn setGridOrder Set or get the ordering of a grid
#' @param G number of geographic dimensions of the setting, typically `G=2` or `G=3`
#' @export
setGridOrder_sp = function(x, G=2){
attr(x,"gridOrder") = gridOrder_gstat(G=G)
return(x)
}
#' @describeIn setGridOrder Set or get the ordering of a grid
#' @export
setGridOrder_array = function(x, G=2){
attr(x,"gridOrder") = gridOrder_array(G=G)
return(x)
}
#' @describeIn setGridOrder Set or get the ordering of a grid
#' @export
gridOrder_sp <- function(G=2){
refpoint = ifelse(G==2, "topleft", "topleftdeep")
cycle = 1:G
return(list(refpoint=refpoint, cycle=cycle))
}
#' @describeIn setGridOrder Set or get the ordering of a grid
#' @export
gridOrder_gstat <- gridOrder_sp
#' @describeIn setGridOrder Set or get the ordering of a grid
#' @export
gridOrder_array <- function(G=2){
refpoint = ifelse(G==2, "bottomleft", "bottomleftdeep")
cycle = 1:G
return(list(refpoint=refpoint, cycle=cycle))
}
#' @describeIn setGridOrder Set or get the ordering of a grid
#' @export
gridOrder_GSLib <-gridOrder_array
#### internal gstatCokriging2something function ------
#' Reorganisation of cokriged compositions
#'
#' Produce compositional predictions out of a [gstat::gstat()] prediction
#'
#' @param COKresult output of a [gstat::predict.gstat()] cokriging,
#' typically of class "data.frame", [sp::SpatialPointsDataFrame()],
#' [sp::SpatialGridDataFrame()] or [sp::SpatialPixelsDataFrame()]
#' @param ... further arguments needed for nscore (**deprecated**)
#'
#' @return an (N,D)-object of class `c("spatialGridAcomp","acomp")`
#' with the predictions, together with an extra attribute "krigVar"
#' containing the cokriging covariance matrices in an (N, D, D)-array; here N=number of
#' interpolated locations, D=number of original components of the composition
#' @aliases gsi.gstatCokriging2rmult
#' @export
#' @seealso [image_cokriged.spatialGridRmult()] for an example
gsi.gstatCokriging2compo <- function(COKresult, ...) UseMethod("gsi.gstatCokriging2compo", COKresult)
#' @describeIn gsi.gstatCokriging2compo Reorganisation of cokriged compositions
#' @method gsi.gstatCokriging2compo default
#' @export
gsi.gstatCokriging2compo.default <- function(COKresult, ...){
stopifnot(is(COKresult, "Spatial"), "data" %in% slotNames(COKresult))
coord = sp::coordinates(COKresult)
res = gsi.gstatCokriging2compo(cbind(as.data.frame(coord), COKresult@data), ...)
if("grid" %in% slotNames(COKresult)){
G = ncol(coord)
dr = 1+(getGridTopology(COKresult)@cellsize<0)
rp = paste(c("top","bottom")[dr[1]], c("left","right")[dr[2]], sep="")
if(G==3) rp = paste(rp, c("deep", "surf")[dr[3]], sep="")
res = setGridOrder(res, refpoint = rp, cycle = 1:G)
}
return(res)
}
#' @describeIn gsi.gstatCokriging2compo Reorganisation of cokriged compositions
#' @method gsi.gstatCokriging2compo data.frame
#' @param V string or matrix describing which logratio was applied ("ilr", "alr",
#' or a matrix computing the ilr corrdinates; clr is not allowed!)
#' @param orignames names of the original components (optional, but recommended)
#' @param tol for generalized inversion of the matrix (**rarely touched!**)
#' @param nscore boolean, were the data normal score-transformed? (**deprecated**)
#' @param gg in the case that normal score transformation was applied, provide the gstat object! (**deprecated**)
#' @export
gsi.gstatCokriging2compo.data.frame = function(COKresult, # output of predict.gstat
V=NULL, # string or matrix describing which logratio was applied
orignames=NULL, # names of the original components (OPTIONAL)
tol=1e-12, # for generalized inversion of the matrix (RARELY USED)
nscore=FALSE, # were the data NS-transformed?
gg=NULL, # in the case that NS was applied, provide the gstat object!
... # further arguments needed for nscore
){
if(is.null(V)) stop("error! V must be either one of the strings 'alr', 'ilr' or 'clr', or else the matrix Psi of coordinate definition!")
cn = colnames(COKresult)
prednames = colnames(COKresult)[grep("pred", cn)]
varnames = colnames(COKresult)[grep("var", cn)]
covnames = colnames(COKresult)[grep("cov", cn)]
coordnames = colnames(COKresult)[-c(grep("cov", cn), grep("var", cn), grep("pred", cn))]
D = length(prednames)+1
if(is.null(orignames) & is.matrix(V)) orignames = rownames(V)
if(is.null(orignames)) orignames = paste("v", 1:D, sep="")
if(length(orignames)!=D) stop("names provided not consistent with number of logratio variables. Did you forget the rest?")
prefix = "ilr"
prediccions = COKresult[,prednames, drop=FALSE]
if(is.character(V)){
if(V=="ilr"){
V = ilrBase(prediccions)
}else if(V=="alr"){
V = rbind(diag(D-1), -1)
prefix = "alr"
}else if(V=="clr"){
V = (diag(ncol(prediccions))-matrix(1/D, ncol=D, nrow=D))[, -D]
prefix = "clr"
}
}
Vsvd = svd(V)
W = with(Vsvd, v %*% diag(ifelse(d>tol, 1/d, 0)) %*% t(u) )
colnames(W) = orignames
if(nscore){
## space to back-transform the predictions
#if(is.null(gg))stop("To apply a nscore backtransformation, the gstat object must be provided!")
#for(i in 1:(D-1)){
# nsc = list()
# nsc$trn.table = attr(gg$data[[i]]$data@data[,i],"trn.table")
# prediccions[,i] = backtr(scores=prediccions[,i], nscore=nsc,...)
#}
stop("gsi.gstatCokriging2compo: use of 'nscore' is deprecated")
}
rg = clrInv(as.matrix(prediccions) %*% W)
# add geographic coordinates as an attribute
attr(rg,"coords") = COKresult[,coordnames]
if(!nscore){
# add cokriging variance matrices as an attribute as well
noms = sub(".pred", "", prednames)
cvmat = array(0, dim=c(nrow(rg), D-1, D-1), dimnames=list(NULL, noms, noms))
vrs = COKresult[,varnames, drop=FALSE]
colnames(vrs) = sub(".var", "", varnames)
for(ivr in noms){
cvmat[ ,ivr, ivr] = vrs[,ivr]
}
cvs = COKresult[,covnames, drop=FALSE]
colnames(cvs) = sub("cov.", "", covnames)
for(ivr in noms){
for(jvr in noms){
if(ivr!=jvr){
dosnoms = c(paste(ivr, jvr, sep="."), paste(jvr, ivr, sep="."))
quin = dosnoms[dosnoms %in% colnames(cvs)]
cvmat[ ,ivr, jvr] = cvs[,quin]
}
}
}
attr(rg,"krigVar") = cvmat
}
class(rg) = c("spatialGridAcomp","acomp")
return(rg)
}
#' @export
gsi.gstatCokriging2rmult <- function(COKresult, ...) UseMethod("gsi.gstatCokriging2rmult", COKresult)
#' @describeIn gsi.gstatCokriging2compo Reorganisation of cokriged multivariate data
#' @method gsi.gstatCokriging2rmult default
#' @export
gsi.gstatCokriging2rmult.default <- function(COKresult, ...){
stopifnot(is(COKresult, "Spatial"), "data" %in% slotNames(COKresult))
coord = sp::coordinates(COKresult)
res = gsi.gstatCokriging2rmult(cbind(coord, COKresult@data), ...)
if("grid" %in% slotNames(COKresult)){
G = ncol(coord)
dr = 1+(getGridTopology(COKresult)@cellsize<0)
rp = paste(c("top","bottom")[dr[1]], c("left","right")[dr[2]], sep="")
if(G==3) rp = paste(rp, c("deep", "surf")[dr[3]], sep="")
res = setGridOrder(res, refpoint = rp, cycle = 1:G)
}
return(res)
}
#' @describeIn gsi.gstatCokriging2compo Reorganisation of cokriged multivariate data
#' @method gsi.gstatCokriging2rmult data.frame
#' @export
gsi.gstatCokriging2rmult.data.frame = function(COKresult, # output of predict.gstat
nscore=FALSE, # were the data NS-transformed?
gg=NULL, # in the case that NS was applied, provide the gstat object!
... # further arguments needed for nscore
){
cn = colnames(COKresult)
prednames = colnames(COKresult)[grep("pred", cn)]
varnames = colnames(COKresult)[grep("var", cn)]
covnames = colnames(COKresult)[grep("cov", cn)]
coordnames = colnames(COKresult)[-c(grep("cov", cn), grep("var", cn), grep("pred", cn))]
D = length(prednames)
noms = sub(".pred", "", prednames)
prediccions = rmult(COKresult[,prednames, drop=FALSE])
if(nscore){
## space to back-transform the predictions
#if(is.null(gg))stop("To apply a nscore backtransformation, the gstat object must be provided!")
#for(i in 1:(D-1)){
# nsc = list()
# nsc$trn.table = attr(gg$data[[i]]$data@data[,i],"trn.table")
# prediccions[,i] = backtr(scores=prediccions[,i], nscore=nsc,...)
#}
stop("gsi.gstatCokriging2compo: use of 'nscore' is deprecated")
}
rg = prediccions
colnames(rg) = noms
# add geographic coordinates as an attribute
attr(rg,"coords") = COKresult[,coordnames]
if(!nscore){
# add cokriging variance matrices as an attribute as well
cvmat = array(0, dim=c(nrow(rg), D, D), dimnames=list(NULL, noms, noms))
vrs = COKresult[,varnames, drop=FALSE]
colnames(vrs) = sub(".var", "", varnames)
for(ivr in noms){
cvmat[ ,ivr, ivr] = vrs[,ivr]
}
cvs = COKresult[,covnames, drop=FALSE]
colnames(cvs) = sub("cov.", "", covnames)
for(ivr in noms){
for(jvr in noms){
if(ivr!=jvr){
dosnoms = c(paste(ivr, jvr, sep="."), paste(jvr, ivr, sep="."))
quin = dosnoms[dosnoms %in% colnames(cvs)]
cvmat[ ,ivr, jvr] = cvs[,quin]
}
}
}
attr(rg,"krigVar") = cvmat
}
class(rg) = c("spatialGridRmult","rmult")
return(rg)
}
#### spatial grids for package "compositions" ------------
# (first attempt: ideally they should be extensions of Spatial)
#' Construct a regionalized composition / reorder compositional simulations
#'
#' Connect some coordinates to a composition (of hard data, of predictions
#' or of simulations); currently, the coordinates
#' are stored in an attribute and the dataset is given a complex S3 class.
#' This functionality **will** change in the future, to make use of package
#' "sp" classes.
#'
#' @param coords coordinates of the locations
#' @param compo (observed or predicted) compositional data set; or else array of
#' simulated compositions
#' @param dimcomp which of the dimensions of `compo` does correspond to the
#' parts of the compositon?
#' @param dimsim if `compo` contains simulations, which of its dimensions does
#' run across the realisations? leave it as NA if `compo` has observations or predictions.
#' @return A (potentially transposed/aperm-ed) matrix of class c("spatialGridAcomp","acomp")
#' with the coordinates in an extra attribute "coords".
#' @seealso [image_cokriged.spatialGridAcomp()] for an example; [gsi.gstatCokriging2compo()] to
#' restructure the output from [gstat::predict.gstat()] confortably
#' @export
spatialGridAcomp = function(coords, compo, dimcomp=2, dimsim=NA){
res = compo
if(is.na(dimsim)){
if(dimcomp==1){
res = t(res)
}
}else{
res = aperm(res, order(c(1,dimsim, dimcomp)))
}
attr(res,"coords") = coords
class(res) = c("spatialGridAcomp","acomp")
return(res)
}
#' Construct a regionalized multivariate data
#'
#' Connect some coordinates to a multivariate data set (of hard data, of predictions
#' or of simulations); currently, the coordinates
#' are stored in an attribute and the dataset is given a complex S3 class.
#' This functionality **will** change in the future, to make use of package
#' "sp" classes.
#'
#' @param coords coordinates of the locations
#' @param data (observed or predicted) rmult or matrix data set; or else array of
#' simulated rmult /real-valued multivariate data
#' @param dimcomp which of the dimensions of `data` does correspond to the
#' variables?
#' @param dimsim if `data` contains simulations, which of its dimensions does
#' run across the realisations? leave it as NA if `data` has observations or predictions.
#' @return A (potentially transposed/aperm-ed) matrix of class c("spatialGridAcomp","acomp")
#' with the coordinates in an extra attribute "coords".
#' @seealso [image_cokriged.spatialGridRmult()] for an example; [gsi.gstatCokriging2rmult()] to
#' restructure the output from [gstat::predict.gstat()] confortably
#' @export
spatialGridRmult = function(coords, data, dimcomp=2, dimsim=NA){
res = data
if(is.na(dimsim)){
if(dimcomp==1){
res = t(res)
}
}else{
res = aperm(res, order(c(1,dimsim, dimcomp)))
}
attr(res,"coords") = coords
class(res) = c("spatialGridRmult","rmult")
return(res)
}
#### image functions for cokriging results ---------
#' Plot an image of gridded data
#'
#' Plot an image of one variable (possibly, one realisation) of output
#' of cokriging or cosimulation functions.
#'
#' @param x object with the interpolated / simulated data; currently there are methods
#' for "spatialGridAcomp" and "spatialGridRmult", but the default method is able to
#' deal with "SpatialPointsDataFrame", "SpatialPixelsDataFrame" and "SpatialGridDataFrame"
#' objects, and with the "data.frame" output of [gstat::predict.gstat()] and
#' [predict_gmSpatialModel]
#' @param ... generic functionality, currently ignored
#'
#' @return Invisibly, a list with elements `breaks` and `col` containing the breaks
#' and hexadecimal colors finally used for the z-values of the image. Particularly
#' useful for plotting other plotting elements on the same color scale.
#' @export
#' @importFrom stats quantile
#' @importFrom graphics layout image title axis abline matplot matlines matpoints
#'
#' @examples
#' \dontrun{
#' getTellus(cleanup=TRUE, TI=TRUE)
#' load("Tellus_TI.RData")
#' head(Tellus_TI)
#' coords = as.matrix(Tellus_TI[,1:2])
#' compo = compositions::acomp(Tellus_TI[,3:7])
#' dt = spatialGridAcomp(coords=coords, compo=compo)
#' image_cokriged(dt, ivar="MgO") # equi-spaced
#' image_cokriged(dt, ivar="MgO", breaks = NULL) # equi-probable
#' }
image_cokriged <- function(x, ...) UseMethod("image_cokriged", x)
#' @describeIn image_cokriged Plot an image of gridded data
#' @param ivar which variable do you want to plot?
#' @param breaks either the approximate number of breaks, or the vector of
#' exact breaks to use for the plotting regions of the chosen variable
#' @param col vector of colors to use for the image
#' @param legendPropSpace which proportion of surface of the device should be used
#' for the legend? trial and error might be necessary to adjust this to your needs
#' @param legendPos where do you want your legend? one of c("top","left","right","bottom")
#' @param main main title for the plot
#' @method image_cokriged default
#' @export
image_cokriged.default <- function(x, ivar=3,
breaks=quantile(as.data.frame(x)[,ivar], probs=c(0:10)/10, na.rm=TRUE),
col = spectralcolors(length(breaks)-1),
legendPropSpace = 0.2, legendPos="top",
main=ifelse(is.character(ivar),ivar,colnames(x)[ivar]),
...){
if(is(x,"SpatialGridDataFrame")){
coords = sp::coordinates(x)
# choose variables
z = x@data[,ivar]
}else if(is(x,"SpatialPixelsDataFrame")){
coords0 = as(x, "SpatialPoints")
gt = getGridTopology(x)
idx = getGridIndex(cc=sp::coordinates(coords0), grid=gt)
# choose variables
z = rep(NA, prod(gt@cells.dim))
z[idx] = x@data[,ivar]
coords = sp::coordinates(gt)
}else if(is(x,"SpatialPointsDataFrame")){
coords0 = as(x, "SpatialPoints")
gt = points2grid(coords0)
idx = getGridIndex(cc=sp::coordinates(coords0), grid=gt)
# choose variables
z = rep(NA, prod(gt@cells.dim))
z[idx] = x@data[,ivar]
coords = sp::coordinates(gt)
}else{
coords = x[,1:2]
# choose variables
z = x[,ivar]
}
if(length(breaks)==1) breaks = pretty(z, n = breaks)
if(length(breaks)-length(col)!=1) col = colorRampPalette(col)(length(breaks)-1)
if(is.logical(legendPos)){
if(legendPos) legendPos = "top"
}
if(!is.logical(legendPos)){
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
# make sure where the color legend goes
par(oma=c(1,1,1,1))
if(legendPos %in% c("top","left")){
ord = 1:2
space = c(legendPropSpace,1-legendPropSpace)
}else if(legendPos %in% c("right","bottom")){
ord = 2:1
space = c(1-legendPropSpace,legendPropSpace)
}else stop("legend must be one of 'top', 'left', 'bottom' or 'right'")
if(legendPos %in% c("top","bottom")){
ord = matrix(ord,ncol=1)
layout(ord,widths=1, heights=space)
a = ifelse(legendPos=="top", 4,1)
par(mar=c(4.5-a,5,a,0.1))
image(breaks, c(0,1), cbind(breaks, breaks), breaks=breaks, col=col, xlab="", ylab="",
yaxt="n", xaxs="i", xlim=range(breaks), oldstyle=TRUE)
par(mar=c(4.5,5,1+a,0.1))
}else{
ord = matrix(ord,nrow=1)
layout(ord,widths=space, heights=1)
par(mar=c(4.5,5,2,0.1))
image(c(0,1), breaks, rbind(breaks, breaks), breaks=breaks, col=col, xlab="", ylab="", xaxt="n",
yaxs="i", ylim=range(breaks), oldstyle=TRUE)
par(mar=c(4.5,5,2,0.1))
}
}
# sort the data so that it can be plotted with image
oo = order(coords[,2], coords[,1])
coords = coords[oo,]
z = z[oo]
# dimensionalizse
xx = sort(unique(coords[,1]))
yy = sort(unique(coords[,2]))
dim(z) = c(length(xx),length(yy))
# plot!
image(xx,yy,z, breaks=breaks, col=col,
xlab=colnames(coords)[1],ylab=colnames(coords)[2], asp=1)
title(main=main, outer=TRUE, line=-1)
invisible(list(breaks=breaks, col=col))
}
#' @describeIn image_cokriged method for spatialGridRmult objects
#' @param isim in case of simulated output, which simulation?
#' @param mask optional mask object if `x` is of class "spatialGridAcomp" or
#' "spatialGridRmult", and they have been masked (see [setMask()])
#' @method image_cokriged spatialGridRmult
#' @export
image_cokriged.spatialGridRmult<- function(x, ivar=1, isim=NULL, breaks=10, mask=attr(x, "mask"),
col = spectralcolors(length(breaks)-1),
legendPropSpace = 0.2, legendPos="top",
main=ifelse(is.character(ivar),ivar,dimnames(x)[[length(dimnames(x))]][ivar]),
...){
# define data, and if necessary, breaks:
if(is.null(isim)){
X = x[,ivar]
}else{
X = x[,isim, ivar]
}
if(!is.null(mask) & inherits(mask,"mask")){
if(!is.null(attr(mask, "fullgrid"))){
X = unmask(data.frame(X), mask=mask)[,1]
coords = sp::coordinates(attr(mask, "fullgrid"))
}else stop("mask should be a 'mask' object with stored fullgrid!")
}else{
coords = attr(x, "coords")
}
oo = order(coords[,2], coords[,1])
X = X[oo]
xx = sort(unique(coords[,1]))
yy = sort(unique(coords[,2]))
dim(X) = c(length(xx),length(yy))
if(is.null(breaks))
breaks = quantile(X, probs=c(0:10)/10, na.rm=TRUE)
if(length(breaks)==1) breaks = pretty(X, n = breaks)
if(length(breaks)-length(col)!=1) col = colorRampPalette(col)(length(breaks)-1)
# make sure where the color legend goes
par(oma=c(1,1,1,1))
if(legendPos %in% c("top","left")){
ord = 1:2
space = c(legendPropSpace,1-legendPropSpace)
}else if(legendPos %in% c("right","bottom")){
ord = 2:1
space = c(1-legendPropSpace,legendPropSpace)
}else stop("legend must be one of 'top', 'left', 'bottom' or 'right'")
if(legendPos %in% c("top","bottom")){
ord = matrix(ord,ncol=1)
layout(ord,widths=1, heights=space)
a = ifelse(legendPos=="top", 4,1)
par(mar=c(4.5-a,5,a,0.1))
image(breaks, c(0,1), cbind(breaks, breaks), breaks=breaks, col=col, xlab="", ylab="",
yaxt="n", xaxs="i", xlim=range(breaks))
par(mar=c(4.5,5,1+a,0.1))
}else{
ord = matrix(ord,nrow=1)
layout(ord,widths=space, heights=1)
par(mar=c(4.5,5,2,0.1))
image(c(0,1), breaks, rbind(breaks, breaks), breaks=breaks, col=col, xlab="", ylab="", xaxt="n",
yaxs="i", ylim=range(breaks))
par(mar=c(4.5,5,2,0.1))
}
image(xx, yy, X, breaks=breaks, col=col,
xlab=colnames(coords)[1],
ylab=colnames(coords)[2], asp=1)
title(main=main, outer=TRUE, line=-1)
invisible(list(breaks=breaks, col=col))
}
#' @describeIn image_cokriged method for spatialGridAcomp objects
#' @param isim in case of simulated output, which simulation?
#' @param mask optional mask object if `x` is of class "spatialGridAcomp" or
#' "spatialGridRmult", and they have been masked (see [setMask()])
#' @method image_cokriged spatialGridAcomp
#' @export
image_cokriged.spatialGridAcomp <- image_cokriged.spatialGridRmult
### grid sorting --------
#' Reorder data in a grid
#'
#' Reorder the data in a compact grid, changing between ordering specifications
#'
#' @param x gridded data
#' @param grid grid topology underlying
#' @param orderIn current ordering description (see [setGridOrder()])
#' @param orderOut desired output ordering description (see [setGridOrder()])
#'
#' @return the data from `x` (typically a matrix), but reordered as `orderOut`
#' @export
#' @seealso [setGridOrder()] for ways of specifying the grid ordering
#' @importClassesFrom sp Spatial
#' @importFrom sp getGridTopology proj4string point.in.polygon
#' @importFrom sp getGridIndex points2grid
#' @examples
#' \dontrun{
#' getTellus(cleanup=TRUE, TI=TRUE)
#' load("Tellus_TI.RData")
#' coords = as.matrix(Tellus_TI[,1:2])
#' compo = compositions::acomp(Tellus_TI[,3:7])
#' dt = spatialGridAcomp(coords=coords, compo=compo)
#' image_cokriged(dt, ivar="MgO", breaks = NULL)
#' x = sort(unique(coords[,1]))
#' y = sort(unique(coords[,2]))
#' x0 = c(min(x), min(y))
#' Ax = c(mean(diff(x)), mean(diff(y)))
#' n = c(length(x), length(y))
#' gr = sp::GridTopology(cellcentre.offset=x0, cellsize=Ax, cells.dim=n)
#' dt0 = sortDataInGrid(Tellus_TI, grid=gr, orderIn=gridOrder_array(2),
#' orderOut=list(refpoint="bottomright", cycle=2:1))
#' coords = as.matrix(dt0[,1:2])
#' compo = compositions::acomp(dt0[,3:7])
#' spatialGridAcomp(coords=coords, compo=compo)
#' }
sortDataInGrid = function(x, grid=attr(x, "grid"),
orderIn=attr(x, "gridOrder"),
orderOut=list(refpoint="bottomleft", cycle=1:2)){
if(is.null(grid)) stop("sortDataInGrid: grid argument cannot be null")
if(is.null(orderIn)) orderIn = list(refpoint = "topleft", cycle = 1:min(dim(as.data.frame(grid))))
if(is(grid,"SpatialGrid")) return(sortDataInGrid(x=x, grid=getGridTopology(grid), orderIn = orderIn, orderOut = orderOut))
if(is(x, "Spatial")) return(sortDataInGrid(x=x@data, grid=grid, orderIn = orderIn, orderOut = orderOut))
if(length(dim(x))==0) x = matrix(x, nrow=length(x), ncol=1)
# indexes array
G = length(orderIn$cycle)
ocycle = (1:G)[orderIn$cycle]
i = aperm( array(1:nrow(x), dim = grid@cells.dim[orderIn$cycle] ), perm = ocycle)
# flip
dd = dim(i)
irow = 1:(dd[1])
icol = 1:(dd[2])
if(length(grep("right",orderIn$refpoint))>0) irow = rev(irow)
if(length(grep("top",orderIn$refpoint))>0) icol = rev(icol)
if(grid@cellsize[1]<0) irow = rev(irow)
if(grid@cellsize[2]<0) icol = rev(icol)
if(G==3){
ideep = 1:(dd[3])
if(length(grep("surf",orderIn$refpoint))>0) ideep = rev(ideep)
if(grid@cellsize[3]<0) ideep = rev(ideep)
i = i[irow, icol, ideep]
}else{
i = i[irow, icol]
}
y = x[i,, drop=T]
if(all(orderOut$cycle==1:G) & orderOut$refpoint=="bottomleft") return(y)
ocycle = (1:G)[orderOut$cycle]
j = aperm( array((1:length(i))[i], dim = grid@cells.dim[orderOut$cycle] ), perm = ocycle)
dd = dim(j)
jrow = 1:(dd[1])
jcol = 1:(dd[2])
if(length(grep("right",orderOut$refpoint))>0) jrow = rev(jrow)
if(length(grep("top",orderOut$refpoint))>0) jcol = rev(jcol)
if(G==3){
jdeep = 1:(dd[3])
if(length(grep("surf",orderOut$refpoint))>0) jdeep = rev(jdeep)
j = j[jrow, jcol, jdeep]
}else{
j = j[jrow, jcol]
}
y = x[j,, drop=T]
return(y)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.