Dynamically creating and manipulating images in ggplot2

This post was inspired by a question on stackoverflow, which wondered whether it was possible to plot a stacked bar chart using rendered images/labels which would overlay their corresponding bar.

The example given was a subset of 'movies' data from IMDB - looking at movies with over 100,000 rating votes and grouping them by year and rating bucket. A basic bar chart representation was shown below:

When I looked at the question (and whether the data could be reshaped into a 'motif' chart), it seems like an interesting exercise to explore a few things:

• Importing and using custom fonts in R
• Dynamically generating images from a list of variables (in this case, years)
• Creating custom grobs which could then be rendered in ggplot

The chart I managed to produce is here, and the commented code below:

require(ggplot2)
require(png)
require(plyr)
require(grid)
require(extrafont)
 
#font_import(pattern="Show") RUN THIS ONCE ONLY
#load the fonts
loadfonts(device="win")
 
#create a subset of data with big votes
big_votes_movies = movies[movies$votes > 100000,]
 
#create a custom palette and append to a table of the unique years (labels) 
years<-data.frame(year=unique(big_votes_movies$year))
palette(rainbow(nrow(years)))
years$col<-palette()
 
#function to create the labels as png files
writeYear<-function(year,col){
 
  png(filename=paste(year,".png",sep=""),width=440,height=190,bg="transparent")
  im<-qplot(1,1,xlab=NULL,ylab=NULL,geom="blank") + 
    geom_text(label=year,size=70, family="Showcard Gothic", color=col,alpha=0.8) +
    theme(axis.text.x = element_blank(),axis.text.y = element_blank()) +
    theme(panel.background = element_rect(fill = "transparent",colour = NA), 
          plot.background = element_rect(fill = "transparent",colour = NA), 
          panel.grid.minor = element_line(colour = "transparent"), 
          panel.grid.major = element_line(colour = "transparent"),
          axis.ticks=element_blank())
  print(im)
  dev.off()
}
 
#call the function to create the placeholder images
apply(years,1,FUN=function(x)writeYear(x["year"],x["col"]))
 
#summarize the data, and create bins manually
summarydata<-big_votes_movies[,c("year","rating","votes")]
summarydata$rating<-cut(summarydata$rating,breaks=c(0,8,8.5,9,Inf),labels=c(0,8,8.5,9))
 
aggdata <- ddply(summarydata, c("year", "rating"), summarise, votes  = sum(votes) )
aggdata<-aggdata[order(aggdata$rating),]
aggdata<-ddply(aggdata,.(rating),transform,ymax=cumsum(votes),ymin=c(0,cumsum(votes))[1:length(votes)])
#identify the image placeholders
aggdata$imgname<-apply(aggdata,1,FUN=function(x)paste(x["year"],".png",sep=""))
ymax<-max(aggdata$ymax)
 
#do the basic plot
z<-qplot(x=10,y=10,geom="blank",xlab="Rating",ylab="Votes \n",main="Big Movie Votes \n") + 
  theme_bw() +
  theme(panel.grid.major = element_line(colour = "transparent"),
        text = element_text(family="Kalinga", size=20,face="bold")        
        ) +
  scale_x_continuous(limits=c(8,9.5)) + 
  scale_y_continuous(limits=c(0,ymax))  
 
#creat a function to create the grobs and return annotation_custom() calls
callgraph<-function(df){
  tiles<-apply(df,1,FUN=function(x)return(annotation_custom(rasterGrob(image=readPNG(x["imgname"]),
                                                      x=0,y=0,height=1,width=1,just=c("left","bottom")),
                                               xmin=as.numeric(x["rating"]),xmax=as.numeric(x["rating"])+0.5,ymin=as.numeric(x["ymin"]),ymax=as.numeric(x["ymax"]))))
  return(tiles)
}
#add the tiles to the base chart
z+callgraph(aggdata)

Created by Pretty R at inside-R.org

Using plot3d (rgl) to render a simple mining block model

This is a follow up to my previous post.

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)

Created by Pretty R at inside-R.org

Quick interpretation of drilling data with R & rgl 3D sphere plots

I have some interest in small scale exploratory mining projects. When you're doing these things (on a private budget anyway) what you're trying to do is build interpreted information as quickly and cheaply as possible. If the project's a dud (or great), you'd like to know as soon as possible.

One of the most useful tools in decision making, especially for feedback to the ongoing project approach, is being able to see a lot of different visualizations of the most current information set you have. I've been looking at rgl for producing 3D plots, and was wondering how it would fit with typical 3d drill data. Up until now this has been possible with specialist proprietary software (or using expensive consultants), but this approach is not optimal in terms of cost, responsiveness or turnaround time.

So I very quickly generated some synthetic drill data, and used plot3d from the rgl library to try and produce a grade map of the sort that would be useful for target indication. It took almost no time at all:

Something I found handy was the ability to tinker with the color, size scaling, and alpha by adding basic parameters to curve or shift the scaling factors.

Playing around with these massively helped visualization (for me anyway) by allowing observation of subtly different summaries of the data.

If you have an OpenGL supporting browser, you can download a spinnable version here. Remember to unzip it.

So as a rudimentary tool for visualization, it's been great, and took no time. As a  follow up to this, I think that I could pretty quickly do the following (and will publish here when I do):

• Add a Shiny Wrapper to make it interactive
• Add some interpolative and sampling functions to mimic the continuous block model generation of commercial packages
• Add other information to map in faults, economic cut-offs etc.

So  in summary I think that in the same way R has been a bootstrap for quants, scientists and modellers, it can really perform the same function for entrepreneurial projects; having only started looking at it a short time ago, I'm amazed at the improved abilities it gives you to handle information.

Here's the commented code for the map (and data generation) for those who are interested:

#########################################################
##  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,-50,-2)),0,30,55)
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+runif(length(vec),-0.1,0.5)
    return(vec)
}
 
scaleColors <- function(colors, vec) {
    ## This function scales a set of colors to any numeric vector
    ## basically creating a discrete scale
    ## and mapping to rgb colors for plot3d to understand
    scaletoUnit <- (vec - min(vec))/diff(range(vec))
    x <- colorRamp(colors)(scaletoUnit)
    rgb(x[,1], x[,2], x[,3], maxColorValue = 255)
}
 
scaleAlpha <- function(vec,maxa=1,mina=0,pwr=1) {
    ## 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
    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
    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))
}
 
#########################################################
##  FUNCTION TO MAKE THE PLOT DATA FOR THE GRAPH       ##
#########################################################
## First in the right shape
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 CALL THE PLOTS                                 ##
#########################################################
plot(dg$grade) # This is a simple 2d plots so I can eyball the distributions
 
myDrillBlock<-plot3d  (x = dg$X, y = dg$Y, z = dg$Z, 
                      col = scaleColors(c("yellow","red"),dg$grade), 
                      size=scaleSize(dg$grade,base=2.4,factor=2.5,pwr=1.5), 
                      alpha=scaleAlpha(dg$grade,pwr=0.8), 
                      xlab="Distance N from reference point (m)",
                      ylab="Distance E from reference point (m)",
                      zlab="Vertical depth from surface (m)",
                      type="s")
 
#########################################################
##  AND HERE YOU CAN CALL OUT TO A LOCAL URL           ##
#########################################################
## This is a complete html OpenGL file format
## If you want to persist or share the graph
## you can save as "Web Page Complete" and share it
## with anyone who has a Chrome browser
 
browseURL(
  paste(
    "file://", 
    writeWebGL(dir=file.path(tempdir(), 
    "myDrillBlock"), width=700),
  sep="")
)

Created by Pretty R at inside-R.org