Merge branch 'master' of gitlab.hzdr.de:geomet/gmGeostats

R/gmAnisotropy.R

@@ -4,22 +4,22 @@
 
 
 #' Convert to anisotropy scaling matrix
-#' 
+#'
 #' Convert an anisotropy specification to a scaling matrix
 #'
 #' @param x an matrix to be tagged as anisotropy scaling matrix
 #'
-#' @return An anisotropy scaling matrix \eqn{A} is such that for any 
-#' lag vector \eqn{h}, the variogram model turns isotropic in terms 
+#' @return An anisotropy scaling matrix \eqn{A} is such that for any
+#' lag vector \eqn{h}, the variogram model turns isotropic in terms
 #' of \eqn{u'=h'\cdot A}. This function does not check any special
-#' property for this matrix! You should probably be using `anis_GSLIBpar2A()` 
+#' property for this matrix! You should probably be using `anis_GSLIBpar2A()`
 #' isntead, and leave `AnisotropyScaling()` for internal uses.
 #'
-#'    
+#'
 #' @export
 #' @family anisotropy
 #' @examples
-#' ( l = anis_GSLIBpar2A(angles=30, ratios=0.5)) 
+#' ( l = anis_GSLIBpar2A(angles=30, ratios=0.5))
 #' ( ll = unclass(l) )
 #' AnisotropyScaling(l)
 AnisotropyScaling = function(x){
@@ -28,28 +28,28 @@ AnisotropyScaling = function(x){
 }
 
 #' Convert to anisotropy scaling matrix
-#' 
+#'
 #' Convert an anisotropy specification to a scaling matrix
 #'
 #' @param x an object convertible to an anisotropy scaling matrix; see details
 #'
-#' @return A matrix \eqn{A} such that for any lag vector \eqn{h}, the variogram model turns 
-#'  isotropic in terms of \eqn{u'=h'\cdot A}. 
+#' @return A matrix \eqn{A} such that for any lag vector \eqn{h}, the variogram model turns
+#'  isotropic in terms of \eqn{u'=h'\cdot A}.
 #'
-#' @details Method `as.AnisotropyScaling.numeric()` expects a vector of two numbers in 2D, 
-#' or a vector of 5 numbers in 3D. These are in 2D, the azimuth of maximum continuity (in 
-#' degrees, clockwise from North) and the anisotropy ratio of short/long range. In 3D 
-#' these are: 1,2) the azimuth and the dip of the direction of maximal continuity; 3) the 
-#' angle of rotation around the axis of the first direction; 4,5) the anisotropy ratios of 
+#' @details Method `as.AnisotropyScaling.numeric()` expects a vector of two numbers in 2D,
+#' or a vector of 5 numbers in 3D. These are in 2D, the azimuth of maximum continuity (in
+#' degrees, clockwise from North) and the anisotropy ratio of short/long range. In 3D
+#' these are: 1,2) the azimuth and the dip of the direction of maximal continuity; 3) the
+#' angle of rotation around the axis of the first direction; 4,5) the anisotropy ratios of
 #' the ranges of the second/first and third/first directions of maximal continuity. All angles
 #' are given in degrees, all ratios must be smaller or equal to 1. This follows gstat (and hence
 #' GSlib) conventions; see gstat::vgm() for details.
-#'    
+#'
 #' @family anisotropy
 #' @export
 #' @aliases as.AnisotropyScaling.numeric
 #' @examples
-#' ( l = anis_GSLIBpar2A(angles=30, ratios=0.5)) 
+#' ( l = anis_GSLIBpar2A(angles=30, ratios=0.5))
 #' ( ll = unclass(l) )
 #' as.AnisotropyScaling(ll)
 #' @export
@@ -60,7 +60,7 @@ as.AnisotropyScaling <- function(x){ UseMethod("as.AnisotropyScaling", x) }
 #' @export
 as.AnisotropyScaling.AnisotropyScaling = function(x){
   x
-} 
+}
 
 
 #' @describeIn as.AnisotropyScaling  from a vector of numbers
@@ -91,9 +91,9 @@ as.AnisotropyScaling.AnisotropyRangeMatrix = function(x){
 
 
 #' Force a matrix to be anisotropy range matrix,
-#' 
-#' Force a matrix M to be considered an anisotropy range matrix, i.e 
-#' with ranges and orientations, 
+#'
+#' Force a matrix M to be considered an anisotropy range matrix, i.e
+#' with ranges and orientations,
 #' such that \eqn{u = sqrt(h' * M^{-1} * h)} allows to use an isotropic
 #' variogram.
 #'
@@ -120,9 +120,9 @@ AnisotropyRangeMatrix = function(x, checkValidity=TRUE){
 }
 
 #' Force a matrix to be anisotropy range matrix,
-#' 
-#' Force a matrix M to be considered an anisotropy range matrix, i.e 
-#' with ranges and orientations, 
+#'
+#' Force a matrix M to be considered an anisotropy range matrix, i.e
+#' with ranges and orientations,
 #' such that \eqn{u = sqrt(h' * M^{-1} * h)} allows to use an isotropic
 #' variogram.
 #'
@@ -156,24 +156,23 @@ as.AnisotropyRangeMatrix.AnisotropyScaling <- function(x) AnisotropyRangeMatrix(
 
 
 #' Produce anisotropy scaling matrix from angle and anisotropy ratios
-#' 
-#' Produce anisotropy matrix (as the transposed of the Cholesky 
+#'
+#' Produce anisotropy matrix (as the transposed of the Cholesky
 #' decomposition) from angle and anisotropy ratios
-#' 
-#' @param ratios vector of two values between 0 and 1 giving the anisotropy ratios of 
+#'
+#' @param ratios vector of two values between 0 and 1 giving the anisotropy ratios of
 #' medium/largest smallest/largest ranges
-#' @param angles, as defined in gstat::vgm (and indeed GSLIB)
-
+#' @param angles as defined in gstat::vgm (and indeed GSLIB). For `anis2D_par2A` 'angle' is the direction of maximum range, i.e. largest spatial continuity, measured clockwise from North
 #' @param inv boolean or integer, see `return` for details
-#' 
-#' @return a 3x3 matrix of anisotropy. 
-#' 
-#' If `inv=TRUE` (or 1) the output is a matrix `A` such that `norm(h %*% A)` 
-#' allows to use isotropic variograms, being `h = c(hx, hy, hz)` the lag vectprs. 
-#' 
-#' If `inv=FALSE` (or 0) the output is a matrix `A` such that `norm(h %*% solve(A))` 
-#' allows to use isotropic variograms.  
-#' 
+#'
+#' @return a 3x3 matrix of anisotropy.
+#'
+#' If `inv=TRUE` (or 1) the output is a matrix `A` such that `norm(h %*% A)`
+#' allows to use isotropic variograms, being `h = c(hx, hy, hz)` the lag vectors.
+#'
+#' If `inv=FALSE` (or 0) the output is a matrix `A` such that `norm(h %*% solve(A))`
+#' allows to use isotropic variograms.
+#'
 #' Other values are meaningless.
 #'
 #' @family anisotropy
@@ -194,11 +193,11 @@ as.AnisotropyRangeMatrix.AnisotropyScaling <- function(x) AnisotropyRangeMatrix(
 
 #' @describeIn anis_GSLIBpar2A 2D case
 #' @param ratio an anisotropy ratio (min/max range)
-#' @param angle direction of maximum range, i.e. largest spatial continuity, measured 
+#' @param angle direction of maximum range, i.e. largest spatial continuity, measured
 #' clockwise from North
 #' @export
 #' @examples
-#' ## ratio=0.5, azimuth 30?? (i.e. direction 60??) 
+#' ## ratio=0.5, azimuth 30?? (i.e. direction 60??)
 #' A = anis2D_par2A(1, 30)
 #' A
 #' AAt = A %*% t(A)
@@ -217,7 +216,7 @@ anis2D_par2A = function(ratio, angle, inv=FALSE){
   ## transposed Cholesky decomp of the rotation matrix from angle of maximum
   #     spatial continuity + ratio min/max ranges
   a = angle * pi/180
-  
+
   # matrix of coordinates of the new basis (in columns)=
   #    =matrix of transformation from new to old coordinate system
   R = matrix(c(sin(a), cos(a), cos(a), -sin(a)), ncol=2)
@@ -235,23 +234,23 @@ anis2D_par2A = function(ratio, angle, inv=FALSE){
 
 #' @describeIn anis_GSLIBpar2A 3D case
 #' @export
-#' @examples 
+#' @examples
 #' c60 = cos(60*pi/180)
 #' s60 = sin(60*pi/180)
 #' c30 = cos(30*pi/180)
 #' s30 = sin(30*pi/180)
 #' #  in the new coordinates, 60cwN is (0,1,0)
 #' R60p = anis3D_par2A(ratios=c(1,1), angles=c(60,0,0))
-#' c(s60, c60, 0) %*% R60p 
+#' c(s60, c60, 0) %*% R60p
 #' R6030 = anis3D_par2A(ratios=c(1,1), angles=c(60,30,0))
 #' # the original X axis is positive on newX and newY, but negative on newZ
-#' c(1,0,0) %*% R6030 
+#' c(1,0,0) %*% R6030
 #' # rotate first direction 60 degrees azimuth, then dip 30degrees upwards
-#' c( c60*c30, -s60*c30, s30) %*% R6030 
+#' c( c60*c30, -s60*c30, s30) %*% R6030
 #' (Ranis = anis3D_par2A(ratios=c(0.5,0.25), angles=c(60,30,0)) )
 anis3D_par2A = function(ratios, angles, inv=FALSE){
   # A should be a matrix such that for h=[hx hy, hz] lag vectpr
-  #   hn = norm( h %*% A[1:ncol(h),]) allows to use the 
+  #   hn = norm( h %*% A[1:ncol(h),]) allows to use the
   #   isotropic correlogram
   ## transposed Cholesky decomp of the rotation matrix from angle of maximum
   #     spatial continuity + ratio min/max ranges

R/variograms.R

@@ -505,12 +505,12 @@ is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){
 
 
 #' Check for any anisotropy class
-#' 
-#' Check that an object contains a valid specification of anisotropy, in any form 
+#'
+#' Check that an object contains a valid specification of anisotropy, in any form
 #'
 #' @param x object to check
 #'
-#' @return a logical, TRUE if the object is an anisotropy specification; FALSE otherwise 
+#' @return a logical, TRUE if the object is an anisotropy specification; FALSE otherwise
 #' @export
 #' @family anisotropy
 #'
@@ -520,7 +520,7 @@ is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){
 #' is.anisotropySpecification(a)
 is.anisotropySpecification = function(x){
   length(grep("Anisotropy", class(x))) |  inherits(x, "Anisotropy")
-} 
+}
 
 
 
@@ -558,7 +558,7 @@ variogram_gmSpatialModel <-  function(object, methodPars=NULL, ...){
 
 # Variogram calculations
 #
-# Compute empricial variograms out of a spatial data object
+# Compute empirical variograms out of a spatial data object
 #
 # @param object spatial data container
 # @param ... further parameters for variogram calculation

vignettes/SimulatingMicrostructures.Rmd

@@ -29,7 +29,7 @@ library(data.table) # only needed, because Solveig is faster with melting in dat
 
 Granites or Gneiss consist mainly of Quartz, Feldspars and micas, e.g. Biotite and/or Muscovite.
 Due to the brittle characteristics these types of rocks have often fractures or joint systems.
-With `gmGeostats` mineral composition, anisotrpy of certain mineral phases and features like joints and fractures can be simulated, e.g. if simulations of crystalline rocks are required for testing code or work flows.
+With `gmGeostats` mineral composition, anisotropy of certain mineral phases and features like joints and fractures can be simulated, e.g. if simulations of crystalline rocks are required for testing code or work flows.
 
 ## General idea - Turning bands algorithm
 
@@ -151,13 +151,13 @@ The function `setCgram` generates these variograms.
 It has four parameters: 
 
 - The `type` provides the variogram model, e.g. Gaussian or Exponential.
-Gaussian models provide smooth transitions to other minerals, and may be used if bands of minerals in equilibrium should be modelled. 
+Gaussian models provide smooth transitions to other minerals, and may be used if bands of minerals in equilibrium should be modeled. 
 Exponential models provide rather abrupt changes between mineral phases and should be preferred model for fractal contacts, such as non-equilibrium contacts, dissolution/porosity, alteration, etc.
-In this example, Quarts and Feldspar are modelled by Gaussion types, and Biotite and the fractures by the Exponential type.
+In this example, Quarts and Feldspar are modeled by Gaussian types, and Biotite and the fractures by the Exponential type.
 - The `nugget` gives the nugget effect, but in the case of microstructures it should be omitted.
 - For the `sill`, the partial sill of the correlation function, the microstructure require simply the matrices with one `1` in the diagonal.
 We produced these matrices in the first part, where the number of phases was defined.
-- The `anisRages` use the information from the anistropy matrices generated in the previous part.
+- The `anisRages` use the information from the anisotropy matrices generated in the previous part.
 This anisotropy range has to be downscaled, therefore a factor is applied prior to the the anisotropy matrices:
 
 ```{r}
@@ -257,7 +257,7 @@ legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)
 
 
 
-If one of the phases should be less pronounced, e.g. here in the example the joints are very frequent, then the potentisal can be regulated by weights.
+If one of the phases should be less pronounced, e.g. here in the example the joints are very frequent, then the potential can be regulated by weights.
 The weights simply increase or diminish the potential of each phase.
 For example:
 ```{r}
@@ -268,7 +268,7 @@ dim(rocksims) = c(resol, resol)
 
 Then the the line is pushed down:
 ```{r, echo=F}
-rock_melt = melt(data.table(rock)[, joints := joints -1], id.vars = c("x", "y", "z"), 
+rock_melt = melt(data.table(rock)[, joints := joints + ph_weights[5]][, qu := qz + ph_weights[1]][, fsp := fsp + ph_weights[2]], id.vars = c("x", "y", "z"), 
                             variable.name = "phase", 
                             value.name = "relation")
 rock_melt[, phase := factor(phase, levels = phases)]
@@ -332,7 +332,7 @@ legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)
 In this plot, the fracture seem to be too big.
 Again, this can be easily resolved by downweighing the fracturesa (and joints) and increase a bit the Quartz:
 ```{r}
-ph_weights = c(.1,0,0,-1.2,-1) # diminish joints and fracture, increase quartz
+ph_weights = c(.5,0,0,-1.2,-1) # diminish joints and fracture, increase quartz
 rocksims = apply(rock2[, -(1:3)],1, function(x) which.max(x + ph_weights))
 dim(rocksims) = c(resol, resol)
 ```