From 425a170e0029cfe16bf5b79c1ee991450e881153 Mon Sep 17 00:00:00 2001
From: Raimon Tolosana <r.tolosana@hzdr.de>
Date: Thu, 13 Apr 2023 15:49:02 +0200
Subject: [PATCH] corrected gmGeostats 0.11.3 version
---
NAMESPACE | 17 ++
NEWS.md | 1 +
R/compositionsCompatibility.R | 30 +-
R/exploratools.R | 1 +
R/gmAnisotropy.R | 2 +-
R/gmSpatialModel.R | 2 +
R/gstatCompatibility.R | 18 +-
R/mask.R | 1 +
R/spSupport.R | 550 +++++++++++++++++-----------------
R/variograms.R | 16 +-
10 files changed, 345 insertions(+), 293 deletions(-)
diff --git a/NAMESPACE b/NAMESPACE
index 84575aa..7dc1d88 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -36,14 +36,29 @@ S3method(as.DataFrameStack,array)
S3method(as.DataFrameStack,data.frame)
S3method(as.DataFrameStack,list)
S3method(as.LMCAnisCompo,CompLinModCoReg)
+S3method(as.LMCAnisCompo,LMCAnisCompo)
+S3method(as.LMCAnisCompo,gmCgram)
S3method(as.LMCAnisCompo,gstat)
S3method(as.LMCAnisCompo,variogramModelList)
S3method(as.array,DataFrameStack)
S3method(as.function,gmCgram)
+S3method(as.gmCgram,LMCAnisCompo)
S3method(as.gmCgram,default)
+S3method(as.gmCgram,variogramModel)
+S3method(as.gmCgram,variogramModelList)
+S3method(as.gmEVario,gstatVariogram)
+S3method(as.gmEVario,logratioVariogram)
+S3method(as.gmEVario,logratioVariogramAnisotropy)
+S3method(as.gmSpatialModel,default)
+S3method(as.gmSpatialModel,gstat)
+S3method(as.gstatVariogram,default)
+S3method(as.gstatVariogram,gmEVario)
S3method(as.gstatVariogram,logratioVariogram)
S3method(as.gstatVariogram,logratioVariogramAnisotropy)
S3method(as.list,DataFrameStack)
+S3method(as.logratioVariogram,gmEVario)
+S3method(as.logratioVariogram,gstatVariogram)
+S3method(as.logratioVariogram,logratioVariogram)
S3method(as.logratioVariogramAnisotropy,default)
S3method(as.logratioVariogramAnisotropy,logratioVariogram)
S3method(as.logratioVariogramAnisotropy,logratioVariogramAnisotropy)
@@ -96,6 +111,7 @@ S3method(noSpatCorr.test,data.frame)
S3method(noSpatCorr.test,default)
S3method(noSpatCorr.test,matrix)
S3method(noStackDim,default)
+S3method(pairsmap,SpatialPointsDataFrame)
S3method(pairsmap,default)
S3method(plot,accuracy)
S3method(plot,gmCgram)
@@ -130,6 +146,7 @@ S3method(swath,default)
S3method(swath,rcomp)
S3method(unmask,DataFrameStack)
S3method(unmask,SpatialPixels)
+S3method(unmask,SpatialPoints)
S3method(unmask,data.frame)
S3method(validate,LeaveOneOut)
S3method(validate,NfoldCrossValidation)
diff --git a/NEWS.md b/NEWS.md
index 3360145..cf2cff7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,6 @@
# gmGeostats 0.11.3
+* (2023-04-13) batch of hidden methods exported to the namespace (mostly for `as.*()` functions converting between different representations of empirical variograms and of geostatistical models)
* (2023-04-04) C routines called as symbols and not name strings; non-ASCII characters removed
# gmGeostats 0.11.2
diff --git a/R/compositionsCompatibility.R b/R/compositionsCompatibility.R
index 551bcc0..7af1848 100644
--- a/R/compositionsCompatibility.R
+++ b/R/compositionsCompatibility.R
@@ -356,14 +356,14 @@ image.logratioVariogramAnisotropy = function(x, jointColor=FALSE, breaks=NULL,
# plot it
if(jointColor){
if(is.function(col)) col=col(length(breaks)-1)
- image.polargrid(r,phi2,z,add=TRUE, breaks=breaks, col=col)
+ image_polargrid(r,phi2,z,add=TRUE, breaks=breaks, col=col)
}else{
if(is.null(breaks)){
if(is.function(col)) col=col(length(breaks)-1)
- image.polargrid(r,phi2,z,add=TRUE, col=col)
+ image_polargrid(r,phi2,z,add=TRUE, col=col)
}else{
if(is.function(col)) col=col(length(probs)-1)
- image.polargrid(r,phi2,z,add=TRUE, col=col, probs=probs)
+ image_polargrid(r,phi2,z,add=TRUE, col=col, probs=probs)
}
}
}
@@ -1102,10 +1102,14 @@ variogramModelPlot.logratioVariogram <- function(vg, model = NULL, # gstat or
## as.gmEVario (empirical) -----
-# @describeIn as.gmEVario
+#' @describeIn as.gmEVario logratioVariogram method not yet available
+#' @method as.gmEVario logratioVariogram
+#' @export
as.gmEVario.logratioVariogram = function(vgemp, ...) stop("not yet available")
-# @describeIn as.gmEVario
+#' @describeIn as.gmEVario logratioVariogramAnisotropy method not yet available
+#' @method as.gmEVario logratioVariogramAnisotropy
+#' @export
as.gmEVario.logratioVariogramAnisotropy = function(vgemp, ...) stop("not yet available")
@@ -1123,11 +1127,17 @@ as.gmEVario.logratioVariogramAnisotropy = function(vgemp, ...) stop("not yet ava
#'
#' @return the same model in the new format.
#' @export
-#' @aliases as.logratioVariogram.logratioVariogram as.logratioVariogram.gmEVario
+# @aliases as.logratioVariogram.logratioVariogram as.logratioVariogram.gmEVario
as.logratioVariogram <- function(vgemp,...) UseMethod("as.logratioVariogram",vgemp)
+#' @describeIn as.logratioVariogram identity method
+#' @method as.logratioVariogram logratioVariogram
+#' @export
as.logratioVariogram.logratioVariogram <- function(vgemp, ...) vgemp
+#' @describeIn as.logratioVariogram method for gmEVario objects (not yet available)
+#' @method as.logratioVariogram gmEVario
+#' @export
as.logratioVariogram.gmEVario = function(vgemp, ...){
stop("not yet available")
}
@@ -1196,6 +1206,7 @@ as.LMCAnisCompo <- function(m, ...) UseMethod("as.LMCAnisCompo",m)
#' @describeIn as.LMCAnisCompo Recast compositional variogram model to format LMCAnisCompo
#' @method as.LMCAnisCompo LMCAnisCompo
+#' @export
as.LMCAnisCompo.LMCAnisCompo <- function(m, ...) m
@@ -1203,6 +1214,7 @@ as.LMCAnisCompo.LMCAnisCompo <- function(m, ...) m
#' @method as.LMCAnisCompo gmCgram
#' @param V eventually, a specification of the way `m` is presently represented
#' @param orignames eventually, vector of names of the components, if `V` is provided and it does not have rownnames
+#' @export
as.LMCAnisCompo.gmCgram = function(m, V=NULL, orignames=rownames(V), ...){
stop("not yet available")
}
@@ -1267,7 +1279,11 @@ gsi.extractCompLMCstructures = function(m){
## as.gmCgram (LMC) -------
-# @describeIn as.gmCgram
+#' @describeIn as.gmCgram method for "LMCAnisCompo" variogram model objects
+#' @param V original logratio matrix used in the definition of "LMCAnisCompo"
+#' @param orignames original variable names of the composition
+#' @export
+#' @method as.gmCgram LMCAnisCompo
as.gmCgram.LMCAnisCompo = function(m, V=attr(m,"contrasts"),
orignames=rownames(m["sill",1]$sill), ...){
# produce constants
diff --git a/R/exploratools.R b/R/exploratools.R
index 468f7fe..fb33fac 100644
--- a/R/exploratools.R
+++ b/R/exploratools.R
@@ -195,6 +195,7 @@ pairsmap <- function(data, ...) UseMethod("pairsmap", data)
#' @describeIn pairsmap Multiple maps
#' @method pairsmap SpatialPointsDataFrame
+#' @export
pairsmap.SpatialPointsDataFrame <- function(data, ...){
pairsmap.default(data@data, loc=sp::coordinates(data), ...)
}
diff --git a/R/gmAnisotropy.R b/R/gmAnisotropy.R
index 8aaede8..bc20e80 100644
--- a/R/gmAnisotropy.R
+++ b/R/gmAnisotropy.R
@@ -66,7 +66,7 @@ as.AnisotropyScaling.AnisotropyScaling = function(x){
#' @describeIn as.AnisotropyScaling from a vector of numbers
#' @method as.AnisotropyScaling numeric
#' @export
-as.AnisotropyScaling.double <-as.AnisotropyScaling.numeric <- function(x){
+as.AnisotropyScaling.numeric <- function(x){
if(length(dim(x))==0){
if(length(x)==2) return(anis2D_par2A(ratio=x[2], angle=x[1], inv=TRUE))
if(length(x)==5) return(anis3D_par2A(ratios=x[4:5], angles=x[1:3], inv=TRUE))
diff --git a/R/gmSpatialModel.R b/R/gmSpatialModel.R
index b11abdc..c5a1d38 100644
--- a/R/gmSpatialModel.R
+++ b/R/gmSpatialModel.R
@@ -357,6 +357,7 @@ as.gmSpatialModel <- function(object, ...) UseMethod("as.gmSpatialModel", object
#' @describeIn as.gmSpatialModel Recast spatial object to gmSpatialModel format
#' @method as.gmSpatialModel default
+#' @export
as.gmSpatialModel.default = function(object, ...) object
@@ -365,6 +366,7 @@ as.gmSpatialModel.default = function(object, ...) object
#' were used to express it? Thsi can be either one string of "alr", "ilr" or "clr", or else a
#' (Dx(D-1))-element matrix of logcontrasts to pass from compositions to logratios
#' @method as.gmSpatialModel gstat
+#' @export
as.gmSpatialModel.gstat = function(object, V=NULL, ...){
stop("as.gmSpatialModel.gstat: not yet implemented")
}
diff --git a/R/gstatCompatibility.R b/R/gstatCompatibility.R
index 90d8ea9..7d0e6eb 100644
--- a/R/gstatCompatibility.R
+++ b/R/gstatCompatibility.R
@@ -344,13 +344,21 @@ variogramModelPlot.gstatVariogram =
## as gmEVario
-# @describeIn as.gmEVario
-# @export
+#' @describeIn as.gmEVario gstatVariogram method not yet available
+#' @method as.gmEVario gstatVariogram
+#' @export
as.gmEVario.gstatVariogram = function(vgemp, ...) stop("not yet available")
## as.logratioVariogram (empirical) -------
# transforms a gstat empirical variogram into a logratioVariogram object (evtl. with anisotropy)
# @describeIn as.logratioVariogram
+#' @describeIn as.logratioVariogram method for gstatVariogram objects
+#' @method as.logratioVariogram gstatVariogram
+#' @export
+#' @param V matrix or name of the logratio transformation used
+#' @param tol tolerance for generalized inverse (eventually for clr case; defaults to 1e-12)
+#' @param orignames names of the original component (default NULL)
+#' @param symmetrize do you want a whole circle of directions? (default FALSE)
as.logratioVariogram.gstatVariogram = function(vgemp, # gstatVariogram object, emprical logratio variogram
V=NULL, # matrix or name of the logratio transformation used
tol=1e-12, # tolerance for generalized inverse (eventually for clr case)
@@ -443,11 +451,13 @@ as.gstatVariogram <- function(vgemp, ...) UseMethod("as.gstatVariogram", vgemp)
#' @describeIn as.gstatVariogram Represent an empirical variogram in "gstatVariogram" format
#' @method as.gstatVariogram default
+#' @export
as.gstatVariogram.default <- function(vgemp, ...) vgemp
-#' @describeIn as.gstatVariogram Represent an empirical variogram in "gstatVariogram" format
+#' @describeIn as.gstatVariogram Represent an empirical variogram in "gstatVariogram" format (not yet available)
#' @method as.gstatVariogram gmEVario
+#' @export
as.gstatVariogram.gmEVario <- function(vgemp,...) stop("not yet available")
#' @describeIn as.gstatVariogram Represent an empirical variogram in "gstatVariogram" format
@@ -707,6 +717,7 @@ as.LMCAnisCompo.variogramModelList <-
#' @describeIn as.gmCgram Convert theoretical structural functions to gmCgram format
#' @method as.gmCgram variogramModelList
+#' @export
as.gmCgram.variogramModelList = function(m, ...){
as.gmCgram(as.LMCAnisCompo(m, ...), ...)
}
@@ -715,6 +726,7 @@ as.gmCgram.variogramModelList = function(m, ...){
#' @describeIn as.gmCgram Convert theoretical structural functions to gmCgram format
#' @method as.gmCgram variogramModel
+#' @export
as.gmCgram.variogramModel = function(m, ...){
# extract nugget
isNugget = m$model=="Nug"
diff --git a/R/mask.R b/R/mask.R
index 53a1b38..c0744b0 100644
--- a/R/mask.R
+++ b/R/mask.R
@@ -458,6 +458,7 @@ unmask.SpatialPixels <- function(x, mask=NULL, fullgrid =attr(mask, "fullgrid"),
#' @describeIn unmask.data.frame Unmask a masked object
#' @method unmask SpatialPoints
#' @importClassesFrom sp SpatialPoints
+#' @export
unmask.SpatialPoints <- function(x, mask=attr(x@data,"mask"),
fullgrid = attr(mask, "fullgrid"),
forceCheck=FALSE, ...){
diff --git a/R/spSupport.R b/R/spSupport.R
index 0f83bf2..7987dde 100644
--- a/R/spSupport.R
+++ b/R/spSupport.R
@@ -1,275 +1,275 @@
-par_remove_readonly = function(x, args=c("cin", "cra", "csi", "cxy", "din", "fig", "page")){
- x[args]<- NULL
- return(x)
-}
-
-
-oneimage = function(x,y,Z,isim,ivar,asp=1,...){
- aux = Z[,ivar,isim]
- dim(aux) = c(length(x),length(y))
- image(x,y,aux,asp=asp,...)
-}
-
-
-#' Write a regionalized data set in GSLIB format
-#'
-#' Write a regionalized data set in plain text GSLIB format
-
-#' @param x regionalized data set
-#' @param file filename
-#' @param header the first line of text for the file, defaults to filename
-#'
-#' @return The status of closing the file, see \code{\link{close}}
-#' for details, although this is seldom problematic. This function is basically called
-#' for its side-effect of writing a data set in the simplified Geo-EAS format that is
-#' used in GSLIB.
-#'
-#' @seealso \url{http://www.gslib.com/gslib_help/format.html}
-#'
-#' @export
-#' @importFrom utils write.table
-#'
-#' @examples
-#' data("jura", package="gstat")
-#' \dontrun{write.GSLib(jura.pred, file="jurapred.txt")}
-write.GSLib <- function(x,file, header=basename(file)) {
- if(is(x, "Spatial")){
- if("data" %in% slotNames(x)) y = x@data else y=NULL
- x = cbind(data.frame(sp::coordinates(x)),y)
- }
- x = as.data.frame(lapply(x,as.numeric))
- x[is.na(x)] = -1e99
- f<-file(file,"w")
- cat(file=f,header,"\n")
- cat(file=f,ncol(x),"\n")
- write(file=f,colnames(x),sep="\n",append=TRUE)
- write.table(file=f,x,append=TRUE,col.names=FALSE,row.names=FALSE)
- close(f)
-}
-
-
-ISATISrotation = function(Az, Ay, Ax){
- myfun = function(a, i, j){
- m = diag(3)
- a = pi*a/180
- m[c(i,j),c(i,j)] = matrix(c(cos(a), sin(a), -sin(a), cos(a)), ncol=2, nrow=2)
- return(m)
- }
- M = myfun(Ax, 2,3) %*% myfun(Ay, 3,1) %*% myfun(Az, 1,2)
- #M = myfun(Az, 1,2) %*% myfun(Ay, 3,1) %*% myfun(Ax, 2,3)
- return(M)
- }
-
-
-
-
-
-
-#################################################
-## function to generate a polarplot representing
-## the number of pairs of data in each lag class
-pointpairs2polargrid = function(loc, #
- maxdist=max(dist(loc))/2, nbins=10, # h maximal value and nr of classes
- dists = seq(0,maxdist,length.out=nbins+1), # h classes to consider
- azimuths=(0:11)*30, # azimuths to consider
- plotit=TRUE){
- # calculate the distances between all points, in radius and in angle
- dr = as.matrix(dist(loc))
- da = outer(1:nrow(loc), 1:nrow(loc), function(i,j) atan2(x=loc[j,2]-loc[i,2], y=loc[j,1]-loc[i,1]) )
- # count the number of radial distances in each class
- distmids = dists[-1] - diff(dists)/2
- comparisonR = sapply(c(dr), function(dd){
- aux = abs(dd-distmids)
- which(aux==min(aux))[1]
- })
- # count the number of angular distances in each class
- angdist = function(a,b) gmApply(cbind(abs(a-b), abs(a-b+2*pi), abs(a-b-2*pi)), 1, min) # construct an angular distance function
- azimuths = azimuths * pi/180
- comparisonA = outer(c(da), azimuths, "-" ) %% (2*pi) # angular distances are modulo 2pi
- comparisonA = gmApply(comparisonA, 1, function(x) which(x==min(x)) )
- # result
- tt = table(comparisonR, comparisonA)
- colnames(tt) = azimuths
- rownames(tt) = distmids
- if(plotit){
- dfaz = diff(azimuths)[1]
- aux = c(azimuths[1]-0.5*dfaz, azimuths+0.5*dfaz)
- aux = pi/2-aux
- image.polargrid(r = dists, phi = aux, tt, breaks = unique(quantile(tt[tt!=0], probs=seq(0,1,0.1))))
- }
- return(tt)
-}
-#################################################
-
-#################################################
-## ensure that azimuths are numeric
-gsi.azimuth2angle = function(x) {
- if(is.numeric(x)) return(x)
- # if "N"
- x = sub("N","",x)
- return(as.numeric(x))
-}
-#################################################
-
-
-
-#################################################
-## internal alternative to paint one single pair
-
-imageOneAnisVario = function(avg, v1=NULL, v2=NULL){
- # construct the polar grid
- r = attr(avg, "dists")
- rlim = range(r)
- phi = as.double(colnames(avg))
- dphi = mean(diff(phi))
- phi = phi - dphi/2
- phi = c(phi,phi[length(phi)]+dphi)*pi/180
- phi2 = pi/2-phi
- # extract the dimension of the composition
- D = 1
- # extract common z levels
- aux = c(unlist(sapply(1:ncol(avg), function(i) avg["vg",i][[1]][,v1,v2])))
- bks = quantile(aux[aux!=0], probs=seq(0,1,0.1), na.rm=TRUE)
- # set the matrix of figures
- plot(c(-1,1)*rlim[2], c(-1,1)*rlim[2],type="n", asp=1, xlab="", ylab="", bty="n", xaxt="n", yaxt="n", ann=FALSE, xaxs="i",yaxs="i")
- # extract the directional variogram for the logratio k vs j
- z = sapply(1:ncol(avg), function(i){
- avg["vg",i][[1]][,v1,v2]
- })
- z[is.nan(z)]=NA
- dim(z) = c(length(r)-1, length(phi)-1)
- # plot it
- image.polargrid(r,phi2,z,add=TRUE)
-}
-#################################################
-
-
-
-#################################################
-## Utility function to
-#### generate a colored polar sector diagram ---------
-image.polargrid = function(
- r = seq(0, 1, length.out = nrow(z)), # radii
- phi = seq(0, 2*pi, length.out = ncol(z)), # angles
- z, # elevation
- zlim = range(z[is.finite(z)], na.rm=TRUE),
- rlim = range(r, na.rm=TRUE), philim = range(phi, na.rm=TRUE),
- col = spectralcolors(length(breaks)-1), # colors to use; by default blue to red
- add = FALSE, xaxs = "i", yaxs = "i",
- probs = seq(0,1,0.1),
- breaks=quantile(z[is.finite(z)], probs=probs), ... # elevation cuts to use; by default, deciles
-){
- #requireNamespace("DescTools", quietly = TRUE)
- # create a plotting region
- if(!add) plot(c(-1,1)*rlim[2], c(-1,1)*rlim[2],type="n", asp=1, xlab="", ylab="")
- # define the grid of polar coordinate indices
- polargrid = expand.grid(1:(length(r)-1), 1:(length(phi)-1))
- # choose the color of each sector
- mycol <- col[cut(c(z),breaks=breaks)]
- dim(mycol) = dim(z)
- # check which is the direction of the series of angles
- ss = sign(diff(phi))
- if(!all(ss>0) & !all(ss<0)) stop("neither clockwise nor counterclockwise angular series")
- # draw counterclockwise
- if(all(ss>0)){
- res <- gmApply( polargrid,1,
- function(i) gsi.DrawAnnulusSector(x = 0, y = 0, radius.in = r[i[1]],
- radius.out = r[i[1]+1],
- angle.beg = phi[i[2]],
- angle.end = phi[i[2]+1],
- col = mycol[i[1],i[2]],
- border=NA)
- )}
- # draw clockwise
- if(all(ss<0)){
- res <- gmApply( polargrid,1,
- function(i) gsi.DrawAnnulusSector(x = 0, y = 0, radius.in = r[i[1]], radius.out = r[i[1]+1],
- angle.end = phi[i[2]], angle.beg = phi[i[2]+1], col = mycol[i[1],i[2]],
- border=NA)
- )}
- invisible(res)
-}
-#################################################
-
-gsi.DrawAnnulusSector = function(x, y, radius.in,
- radius.out,
- angle.beg,
- angle.end,
- ...){
- if(requireNamespace("DescTools", quietly=TRUE)){
- DescTools::DrawCircle(x=x, y=y, r.out = radius.out, r.in=radius.in, theta.1 = angle.beg,
- theta.2 = angle.end, ...)
- # }else if(requireNamespace("circlize", quietly=TRUE)){
- # circlize::draw.sector(
- # start.degree = 0,
- # end.degree = 360,
- # rou1 = 1,
- # rou2 = NULL,
- # center = c(0, 0),
- # clock.wise = TRUE,
- # col = NA,
- # border = "black",
- # lwd = par("lwd"),
- # lty = par("lty"))
- }else stop("DrawAnnulusSector: one of the packages 'Desctools' or 'circlize' is needed. Install one!")
-
-}
-
-gsi.DrawCircle = function(x, y, r.in,
- r.out,
- theta.1,
- theta.2,
- ...){
- if(requireNamespace("DescTools", quietly=TRUE)){
- DescTools::DrawCircle(x=x, y=y, r.out = r.out, r.in=r.in, theta.1 = theta.1,
- theta.2 = theta.2, ...)
- # }else if(requireNamespace("circlize", quietly=TRUE)){
- # circlize::draw.sector(
- # start.degree = theta.2,
- # end.degree = theta.1,
- # rou1 = r.out,
- # rou2 = r.in,
- # center = c(x, y),
- # clock.wise = FALSE,
- # ...)
- }else stop("DrawCircle: one of the packages 'Desctools' or 'circlize' is needed. Install one!")
-}
-
-
-#### general internal functions -------------------------------
-
-gsi.powM <- function(A,alpha=1/2) {
- with(svd(A),u%*%diag(d^alpha)%*%t(v))
-}
-
-
-
-gsiDiag <- function(d) {
- if( length(d)>1)
- return(diag(d))
- return( structure(d,dim=c(1,1)))
-}
-
-gsiInv <- function(A,tol=1E-15) {
- with(svd(A),v%*% gsiDiag(ifelse(abs(d)/max(d)>tol,1/d,0))%*%t(u))
-}
-
-#### commons for C-R mixed code ---------------------------
-
-checkInt <- function(x,n) {
- x <- as.integer(x)
- if( !missing(n) ) stopifnot(length(x)==n)
- x
-}
-
-checkDouble <- function(x,n) {
- x <- as.numeric(x)
- if( !missing(n) ) stopifnot(length(x)==n)
- x
-}
-
-checkLogical <- function(x,n) {
- x <- as.logical(x)
- if( !missing(n) ) stopifnot(length(x)==n)
- x
-}
+par_remove_readonly = function(x, args=c("cin", "cra", "csi", "cxy", "din", "fig", "page")){
+ x[args]<- NULL
+ return(x)
+}
+
+
+oneimage = function(x,y,Z,isim,ivar,asp=1,...){
+ aux = Z[,ivar,isim]
+ dim(aux) = c(length(x),length(y))
+ image(x,y,aux,asp=asp,...)
+}
+
+
+#' Write a regionalized data set in GSLIB format
+#'
+#' Write a regionalized data set in plain text GSLIB format
+
+#' @param x regionalized data set
+#' @param file filename
+#' @param header the first line of text for the file, defaults to filename
+#'
+#' @return The status of closing the file, see \code{\link{close}}
+#' for details, although this is seldom problematic. This function is basically called
+#' for its side-effect of writing a data set in the simplified Geo-EAS format that is
+#' used in GSLIB.
+#'
+#' @seealso \url{http://www.gslib.com/gslib_help/format.html}
+#'
+#' @export
+#' @importFrom utils write.table
+#'
+#' @examples
+#' data("jura", package="gstat")
+#' \dontrun{write.GSLib(jura.pred, file="jurapred.txt")}
+write.GSLib <- function(x,file, header=basename(file)) {
+ if(is(x, "Spatial")){
+ if("data" %in% slotNames(x)) y = x@data else y=NULL
+ x = cbind(data.frame(sp::coordinates(x)),y)
+ }
+ x = as.data.frame(lapply(x,as.numeric))
+ x[is.na(x)] = -1e99
+ f<-file(file,"w")
+ cat(file=f,header,"\n")
+ cat(file=f,ncol(x),"\n")
+ write(file=f,colnames(x),sep="\n",append=TRUE)
+ write.table(file=f,x,append=TRUE,col.names=FALSE,row.names=FALSE)
+ close(f)
+}
+
+
+ISATISrotation = function(Az, Ay, Ax){
+ myfun = function(a, i, j){
+ m = diag(3)
+ a = pi*a/180
+ m[c(i,j),c(i,j)] = matrix(c(cos(a), sin(a), -sin(a), cos(a)), ncol=2, nrow=2)
+ return(m)
+ }
+ M = myfun(Ax, 2,3) %*% myfun(Ay, 3,1) %*% myfun(Az, 1,2)
+ #M = myfun(Az, 1,2) %*% myfun(Ay, 3,1) %*% myfun(Ax, 2,3)
+ return(M)
+ }
+
+
+
+
+
+
+#################################################
+## function to generate a polarplot representing
+## the number of pairs of data in each lag class
+pointpairs2polargrid = function(loc, #
+ maxdist=max(dist(loc))/2, nbins=10, # h maximal value and nr of classes
+ dists = seq(0,maxdist,length.out=nbins+1), # h classes to consider
+ azimuths=(0:11)*30, # azimuths to consider
+ plotit=TRUE){
+ # calculate the distances between all points, in radius and in angle
+ dr = as.matrix(dist(loc))
+ da = outer(1:nrow(loc), 1:nrow(loc), function(i,j) atan2(x=loc[j,2]-loc[i,2], y=loc[j,1]-loc[i,1]) )
+ # count the number of radial distances in each class
+ distmids = dists[-1] - diff(dists)/2
+ comparisonR = sapply(c(dr), function(dd){
+ aux = abs(dd-distmids)
+ which(aux==min(aux))[1]
+ })
+ # count the number of angular distances in each class
+ angdist = function(a,b) gmApply(cbind(abs(a-b), abs(a-b+2*pi), abs(a-b-2*pi)), 1, min) # construct an angular distance function
+ azimuths = azimuths * pi/180
+ comparisonA = outer(c(da), azimuths, "-" ) %% (2*pi) # angular distances are modulo 2pi
+ comparisonA = gmApply(comparisonA, 1, function(x) which(x==min(x)) )
+ # result
+ tt = table(comparisonR, comparisonA)
+ colnames(tt) = azimuths
+ rownames(tt) = distmids
+ if(plotit){
+ dfaz = diff(azimuths)[1]
+ aux = c(azimuths[1]-0.5*dfaz, azimuths+0.5*dfaz)
+ aux = pi/2-aux
+ image_polargrid(r = dists, phi = aux, tt, breaks = unique(quantile(tt[tt!=0], probs=seq(0,1,0.1))))
+ }
+ return(tt)
+}
+#################################################
+
+#################################################
+## ensure that azimuths are numeric
+gsi.azimuth2angle = function(x) {
+ if(is.numeric(x)) return(x)
+ # if "N"
+ x = sub("N","",x)
+ return(as.numeric(x))
+}
+#################################################
+
+
+
+#################################################
+## internal alternative to paint one single pair
+
+imageOneAnisVario = function(avg, v1=NULL, v2=NULL){
+ # construct the polar grid
+ r = attr(avg, "dists")
+ rlim = range(r)
+ phi = as.double(colnames(avg))
+ dphi = mean(diff(phi))
+ phi = phi - dphi/2
+ phi = c(phi,phi[length(phi)]+dphi)*pi/180
+ phi2 = pi/2-phi
+ # extract the dimension of the composition
+ D = 1
+ # extract common z levels
+ aux = c(unlist(sapply(1:ncol(avg), function(i) avg["vg",i][[1]][,v1,v2])))
+ bks = quantile(aux[aux!=0], probs=seq(0,1,0.1), na.rm=TRUE)
+ # set the matrix of figures
+ plot(c(-1,1)*rlim[2], c(-1,1)*rlim[2],type="n", asp=1, xlab="", ylab="", bty="n", xaxt="n", yaxt="n", ann=FALSE, xaxs="i",yaxs="i")
+ # extract the directional variogram for the logratio k vs j
+ z = sapply(1:ncol(avg), function(i){
+ avg["vg",i][[1]][,v1,v2]
+ })
+ z[is.nan(z)]=NA
+ dim(z) = c(length(r)-1, length(phi)-1)
+ # plot it
+ image_polargrid(r,phi2,z,add=TRUE)
+}
+#################################################
+
+
+
+#################################################
+## Utility function to
+#### generate a colored polar sector diagram ---------
+image_polargrid = function(
+ r = seq(0, 1, length.out = nrow(z)), # radii
+ phi = seq(0, 2*pi, length.out = ncol(z)), # angles
+ z, # elevation
+ zlim = range(z[is.finite(z)], na.rm=TRUE),
+ rlim = range(r, na.rm=TRUE), philim = range(phi, na.rm=TRUE),
+ col = spectralcolors(length(breaks)-1), # colors to use; by default blue to red
+ add = FALSE, xaxs = "i", yaxs = "i",
+ probs = seq(0,1,0.1),
+ breaks=quantile(z[is.finite(z)], probs=probs), ... # elevation cuts to use; by default, deciles
+){
+ #requireNamespace("DescTools", quietly = TRUE)
+ # create a plotting region
+ if(!add) plot(c(-1,1)*rlim[2], c(-1,1)*rlim[2],type="n", asp=1, xlab="", ylab="")
+ # define the grid of polar coordinate indices
+ polargrid = expand.grid(1:(length(r)-1), 1:(length(phi)-1))
+ # choose the color of each sector
+ mycol <- col[cut(c(z),breaks=breaks)]
+ dim(mycol) = dim(z)
+ # check which is the direction of the series of angles
+ ss = sign(diff(phi))
+ if(!all(ss>0) & !all(ss<0)) stop("neither clockwise nor counterclockwise angular series")
+ # draw counterclockwise
+ if(all(ss>0)){
+ res <- gmApply( polargrid,1,
+ function(i) gsi.DrawAnnulusSector(x = 0, y = 0, radius.in = r[i[1]],
+ radius.out = r[i[1]+1],
+ angle.beg = phi[i[2]],
+ angle.end = phi[i[2]+1],
+ col = mycol[i[1],i[2]],
+ border=NA)
+ )}
+ # draw clockwise
+ if(all(ss<0)){
+ res <- gmApply( polargrid,1,
+ function(i) gsi.DrawAnnulusSector(x = 0, y = 0, radius.in = r[i[1]], radius.out = r[i[1]+1],
+ angle.end = phi[i[2]], angle.beg = phi[i[2]+1], col = mycol[i[1],i[2]],
+ border=NA)
+ )}
+ invisible(res)
+}
+#################################################
+
+gsi.DrawAnnulusSector = function(x, y, radius.in,
+ radius.out,
+ angle.beg,
+ angle.end,
+ ...){
+ if(requireNamespace("DescTools", quietly=TRUE)){
+ DescTools::DrawCircle(x=x, y=y, r.out = radius.out, r.in=radius.in, theta.1 = angle.beg,
+ theta.2 = angle.end, ...)
+ # }else if(requireNamespace("circlize", quietly=TRUE)){
+ # circlize::draw.sector(
+ # start.degree = 0,
+ # end.degree = 360,
+ # rou1 = 1,
+ # rou2 = NULL,
+ # center = c(0, 0),
+ # clock.wise = TRUE,
+ # col = NA,
+ # border = "black",
+ # lwd = par("lwd"),
+ # lty = par("lty"))
+ }else stop("DrawAnnulusSector: one of the packages 'Desctools' or 'circlize' is needed. Install one!")
+
+}
+
+gsi.DrawCircle = function(x, y, r.in,
+ r.out,
+ theta.1,
+ theta.2,
+ ...){
+ if(requireNamespace("DescTools", quietly=TRUE)){
+ DescTools::DrawCircle(x=x, y=y, r.out = r.out, r.in=r.in, theta.1 = theta.1,
+ theta.2 = theta.2, ...)
+ # }else if(requireNamespace("circlize", quietly=TRUE)){
+ # circlize::draw.sector(
+ # start.degree = theta.2,
+ # end.degree = theta.1,
+ # rou1 = r.out,
+ # rou2 = r.in,
+ # center = c(x, y),
+ # clock.wise = FALSE,
+ # ...)
+ }else stop("DrawCircle: one of the packages 'Desctools' or 'circlize' is needed. Install one!")
+}
+
+
+#### general internal functions -------------------------------
+
+gsi.powM <- function(A,alpha=1/2) {
+ with(svd(A),u%*%diag(d^alpha)%*%t(v))
+}
+
+
+
+gsiDiag <- function(d) {
+ if( length(d)>1)
+ return(diag(d))
+ return( structure(d,dim=c(1,1)))
+}
+
+gsiInv <- function(A,tol=1E-15) {
+ with(svd(A),v%*% gsiDiag(ifelse(abs(d)/max(d)>tol,1/d,0))%*%t(u))
+}
+
+#### commons for C-R mixed code ---------------------------
+
+checkInt <- function(x,n) {
+ x <- as.integer(x)
+ if( !missing(n) ) stopifnot(length(x)==n)
+ x
+}
+
+checkDouble <- function(x,n) {
+ x <- as.numeric(x)
+ if( !missing(n) ) stopifnot(length(x)==n)
+ x
+}
+
+checkLogical <- function(x,n) {
+ x <- as.logical(x)
+ if( !missing(n) ) stopifnot(length(x)==n)
+ x
+}
diff --git a/R/variograms.R b/R/variograms.R
index 21b6b54..dbe3e5f 100644
--- a/R/variograms.R
+++ b/R/variograms.R
@@ -1002,11 +1002,13 @@ plot.gmEVario = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo=
#' (for compositional data) and \code{gstatVariogram}
#' (from package \code{gstat}).
#' @export
-#' @aliases as.gmEVario.gstatVariogram as.gmEVario.logratioVariogram
-#' as.gmEVario.logratioVariogramAnisotropy
-#'
#' @family gmEVario functions
+# @aliases as.gmEVario.gstatVariogram as.gmEVario.logratioVariogram
+# as.gmEVario.logratioVariogramAnisotropy
as.gmEVario <- function(vgemp,...){ UseMethod("as.gmEVario",vgemp)}
+
+#' @describeIn as.gmEVario default method
+#' @method as.gmEVario default
as.gmEVario.default <- function(vgemp,...) vgemp
@@ -1182,11 +1184,11 @@ print.directionClass = function(x, complete=TRUE, ...){
-
-as.directorVector <- function(x){ UseMethod("as.directorVector",x) }
+# @export
+as.directorVector <- function(x, ...){ UseMethod("as.directorVector",x) }
#' @method as.directorVector default
-as.directorVector.default = function(x,...) x
+as.directorVector.default = function(x) x
#' @method as.directorVector azimuth
as.directorVector.azimuth = function(x, D=2){
@@ -1199,7 +1201,7 @@ as.directorVector.azimuth = function(x, D=2){
}
#' @method as.directorVector azimuthInterval
-as.directorVector.azimuthInterval = function(x, D=2){
+as.directorVector.azimuthInterval = function(x, D=2, ...){
res = (x[[1]]+x[[2]])/2
return(as.directorVector.azimuth(res))
}
--
GitLab