I'd initially put together a plot to visualize surface mining targets based on native drill hole information. I thought at the time it would be fairly easy to extend this to a block model with cut-offs, and so here's an early iteration of it:
There's a zip file of a spinnable version here. You'll need an OpenGL compatible browser like Chrome, and remember to unzip it (htm file and *_files folder to the same directory). Also remember to download from Google Drive you need to click File>>Download, or Ctrl+S
I extended the original simple plot in the following ways:
1) Tweaked the data generator to make the simulated data look "right"
2) Built a basic bucketing function which splits the area into blocks, percolates in some surrounding data, and assigns a projected grade using approx() - clearly this is simplistic, but it served to demonstrate the point while I decide on the best interp method (currently looking at gstat)
3) Modified the alpha function so that I could add a "cut-off" grade, i.e. hide any block with a projected grade below a certain amount.
4) Added a level function in the color call, so that I could use 2 separate gradients across the range.
5) Used shade3d() to draw individual cubes representing each block
The commented code is below:
######################################################### ## PRELIMINARIES ## ######################################################### ## load required packages require(rgl) # 3d package require(RColorBrewer) # color package for ramp creation ######################################################### ## MAKE TEST DATA ## ######################################################### ## make some data that looks like mining data for drilling of a surface ore body ## this represents 144 holes drilled at 25m intervals on a 300 m square grid ## each hole is 50 m long (sampled at 2m intervals) ## each hols in drilled on a bearing of 30%, and at 55% from surface ## obviously in the real world you'd have real data which is why these repeat ## the ore grade is set to 0 as a placeholder dg<-data.frame(expand.grid(seq(0,300,25),seq(0,300,25),seq(0,-100,-3)),0,30,45) colnames(dg)<-c("X","Y","Z","grade","bearing","downholedeg") ######################################################### ## FUNCTIONS TO PREPARE IN ADVANCE ## ######################################################### genGrade<-function(vec,g){ ## This function generates grade data, clustering higher values ## in the middle of the series (what you'd expect in a mine), ## and adding some random noise ## Also allows shift on x, y & z with input g ## basically it's just making dirty x^y curves and centering them vec<-vec-((max(vec)-min(vec))/2)-min(vec) vec<-abs(vec)+g vec<-1/vec^g vec<-vec*runif(length(vec),1,2) vec<-vec/max(vec) vec<-(vec^1.5)+runif(length(vec),0,5) return(vec) } scaleColors <- function(colors, vec, floor=0, floorCol=NULL) { ## This function generates block colors, based on a scale ## By default it uses one gradient ## But there's an option to a subordinate "floor" series scaletoUnit <- (vec - min(vec))/diff(range(vec)) floor<-max((floor-min(vec))/diff(range(vec)),0) floor<-min(floor,max(scaletoUnit)) colFrame<-data.frame("val"=scaletoUnit,"1"=0,"2"=0,"3"=0) colFrame[which(colFrame$val>=floor),2:4] <- colorRamp(colors)(colFrame[which(colFrame$val>=floor),1]) ifelse(is.null(floorCol), print("NULL"), colFrame[which(colFrame$val<=floor),2:4] <- colorRamp(colors)(colFrame[which(colFrame$val<=floor),1])) rgb(colFrame[,2], colFrame[,3], colFrame[,4], maxColorValue = 256) } scaleAlpha <- function(vec,maxa=1,mina=0,pwr=1,floor=NULL) { ## This function creates a discrete scale which maps to alpha values ## typically 0 to 1, although it allows for a superset or subset ## and has a raise to power factor to increase or decrease the hi/lo differential ## it also has a "floor", block below this are assigned an alpha of 0 to hide them ## which shows a visual "cut-off" ifelse(is.null(floor),,vec[which(vec<=floor)]<-min(vec)) scaletoUnit <- (vec - min(vec))/diff(range(vec)) x<-min(mina,maxa)+(scaletoUnit*abs(maxa-mina)) x<-x^pwr return(x) } scaleSize <- function(vec,base=0,factor=1,pwr=1) { ## This function creates a discrete scale for size ## Again the power factor allows you to de/emphasize the tail ## Not used in the block model scaletoUnit <- (vec - min(vec))/diff(range(vec)) x<-base+(scaletoUnit*factor)^pwr return(x) } absCoords<-function(xrel,yrel,zdepth,hangle,zangle,xabs=0,yabs=0){ ## This converts the x(rel),y(rel),z(depth) numbers to useful "real" co-ordinates realz<-zdepth*sin(zangle) realx<-(-realz*sin(zangle)*cos(hangle))+xrel+xabs realy<-(-realz*sin(zangle)*sin(hangle))+yrel+yabs return(data.frame(X=realx,Y=realy,Z=realz)) } blockGrade<-function(i, percFactor=1){ ## This function finds all of the grade points in each block ## or which sit within a percolation envelope outside ## and generates a grade based on the approx() function ## which is essentially an average of 50 fitted points ## to the implied grade distribution ## bit of a hack, I know pf<-percFactor # percolation factor must be >=1 derivedGrade<-0 lx<-min(dg$X)+((blockGrid$X[i]-1)*blockSize[1])-(pf-1)*blockSize[1] ly<-min(dg$Y)+((blockGrid$Y[i]-1)*blockSize[2])-(pf-1)*blockSize[2] lz<-min(dg$Z)+((blockGrid$Z[i]-1)*blockSize[3])-(pf-1)*blockSize[3] ux<-lx+blockSize[1]*pf uy<-ly+blockSize[2]*pf uz<-lz+blockSize[3]*pf gradeBucket<-dg[ which((dg$X<ux) & (dg$X>lx) & (dg$Y<uy) & (dg$Y>ly) & (dg$Z<uz) & (dg$Z>lz)),4] if (length(gradeBucket)>=2) {derivedGrade<-mean(approx(gradeBucket)$y)} return(derivedGrade) } ######################################################### ## MAKE THE PLOT DATA FOR THE GRAPH ## ######################################################### ## Call the function to generate grade series dg$grade<-genGrade(dg$X,2)*genGrade(dg$Y,3)*genGrade(dg$Z,4) ## Then scale it to 0-20 maximum (this is arbitrary, still creating data) dg$grade<-abs(dg$grade)*(20/max(dg$grade)) ## Then call the function to make the co-ordinates abolute dg[,1:3]<-absCoords(xrel=dg$X,yrel=dg$Y,zdepth=dg$Z, hangle=(2*pi*dg$bearing/360),zangle=(2*pi*dg$downholedeg/360), xabs=208,yabs=1102) ## Now create the blocks blockSize<-c(30,30,10) # <<<<<< define block size ####################################### xBlocks<-diff(range(dg$X))/blockSize[1] yBlocks<-diff(range(dg$Y))/blockSize[2] zBlocks<-diff(range(dg$Z))/blockSize[3] ## Make a data frame to hold the block data blockGrid<-data.frame( expand.grid(seq(1,xBlocks+1,1), seq(1,yBlocks+1,1), seq(1,zBlocks+1,1)) ,0,0,"") colnames(blockGrid)<-c("X","Y","Z","Grade","Alpha","Color") ## Then fill the blocks with the interpreted grades, colors and alpha ## GRADES for(i in 1:nrow(blockGrid)){blockGrid$Grade[i]<-blockGrade(i,percFactor=1.2)} ## COLOR blockGrid$Color<-scaleColors(c("yellow","red"), blockGrid$Grade, floor=4,floorCol=c("white","yellow")) ## ALPHA blockGrid$Alpha<-scaleAlpha(blockGrid$Grade,pwr=2.2,floor=2.5) blockGrid$Z<-(-blockGrid$Z) ######################################################### ## NOW CALL THE PLOTS ## ######################################################### ## THE BASIC PLOT plot3d(x=blockGrid$X,y=blockGrid$Y,z=blockGrid$Z, type="p", size=1, xlab="NS Axis - block co-ordinate", ylab="EW Axis - block co-ordinate", zlab="Vertical depth (blocks)") ## DRAW THE CUBES for(i in 1:nrow(blockGrid)) { shade3d(translate3d(cube3d(), blockGrid$X[i], blockGrid$Y[i], blockGrid$Z[i]), col=blockGrid$Color[i], alpha = blockGrid$Alpha[i]) } ## CHANGE INITIAL ORIENTATION rgl.viewpoint( theta = -15, phi = -85, scale=c(blockSize[1],blockSize[2],blockSize[3]), zoom=0.5)
No comments:
Post a Comment