Happy new year to all ! For this first post of 2017, let's dig in some techincal details on graphs in GCDkit, for those of you who would like to develop their own templates. Interested? Read on....
One of the beauties of GCDkit is the possibility to access many
predefined diagrams of common use. They reside in the
"plots/classification" or "plots/geotectonic" menus. The difference
between the two types of diagrams is subtle; diagrams from
"classification" can be used for to classify samples according to where
they fall in the diagram, whereas those from the "geotectonic" group
cannot be used in this way. In some cases it is rather obvious (some of
the geotectonic diagrams actually include multiple panels, so the
classification would be a bit more difficult). In other cases the
difference is just a matter of how the function is written internally,
and converting a geotectonic diagram to a classification diagram would
be easy.
Diagram templates
Let's start by having a look at how
diagrams are defined in GCDkit. You may have noticed that, when you load
data, GCDkit prints a list of diagrams that can be used – i.e. the
program reads diagrams from files on the disk, and checks for each of
them whether the required data exists for this dataset (technically, if
WR[] has the appropriate column headings). So adding one diagram
template is actually a simple matter of writing a template file.
Template
files reside (on a windows 7 install) in C:\Program
Files\R\R-3.2.1\library\GCDkit\Diagrams\ . If you navigate there (or in
the …\library\GCDkit directory in your R installation, if you did put R
somewhere else, for some reason), you will notice 4 "spider" something
files (which we are not going to consider for now); a classification
directory; a geotectonic directory. The "classification" directory has 3
sub-directories for the 3 languages now supported. Each of the 3
language directories contains the same files, with just the names being
translated. In fact the master directory is the English one, and the
other two are automatically generated from a dictionary file (.dict.txt)
during compilation – not something you'll have access to, so for the
time being better stick to English.
Let's look at a diagram
template in … \library\GCDkit\Diagrams\Classification\English . For
instance a simple one, a Shand diagram. The corresponding file, Shand.r,
can be opened in a text editor of your choice (notepad is by default on
your computer but if you are serious about programming, you should look
for one of the myriad of more clever alternatives – notepad++,
ultraedit, winedt, tinn-R…). Shand.r is a 28-lines file, that we will
now examine.
The first line reads
Shand<-function(){
So
we are defining a new function, named
Shand(). The function takes no
arguments. The curly brace { closes at the very end of the file, line
28, so all of the file is enclosed in the function.
Next are two lines:
x.data<<-WR[,"A/CNK"]
y.data<<-WR[,"A/NK"]
Here the code defines two variables
(vectors) called
x.data and
y.data; they are a subset of the matrix
WR,
specifically
x.data gets the content of the column called
"A/CNK" and
y.data gets
"A/NK". Both are calculated by GCDkit from your major
elements data when loading a file, so you can trust they'll always be
there – otherwise you would have to do the calculation yourself. Note
that we use here a "super assignment",
<<- , to make sure the
variables are accessible even outside of the function (ok, you do not
really need to understand that if you are just up to hacking some
diagram templates).
Then comes the core of the graph function, the
template. The template will be used by GCDkit's internal plotting
system called "Figaro", which is a rather complex affair that allows
interactions with the graphs. You do not really need to understand the
specifics of it; all you need to know is that you have to supply an
object called
temp (for template, not temporary!) that will be a list
containing all the stuff you'd like Figaro to plot on your diagram (in
addition to the datapoints): lines, labels, etc.
Here, the template is assembled lines 21 to 25, with the following code:
if(getOption("gcd.plot.text")){
temp<-c(temp1,temp2)}
else{
temp<-temp1
}
Which looks at the value of the (Boolean) option
gcd.plot.text, and accordingly will either merge
temp1 and
temp2, or use
only
temp1, to create the list
temp.
temp2, as you may have guessed,
contains the text labels as we shall see latter on.
First, we define the lines with
temp1<-list(
clssf=list("NULL",use=2:5,rcname=c("metaluminous","peraluminous","peralkaline","undefined")),
lines1=list("lines",x=c(1,1,0.5,0.5,1),y=c(1,7,7,1,1)),
lines2=list("lines",x=c(1,2,2,1,1),y=c(1,1,7,7,1)),
lines3=list("lines",x=c(1,1,0.5,0.5,1),y=c(1,0,0,1,1)),
lines4=list("lines",x=c(1,2,2,1,1),y=c(1,1,0,0,1)),
lines11=list("abline",h=1,col=plt.col[2]),
lines12=list("abline",v=1,col=plt.col[2]),
lines13=list("lines",x=c(0,7),y=c(0,7),col=plt.col[2],lty="dashed"),
GCDkit=list("NULL",plot.type="binary",plot.position=14,plot.name="A/CNK - A/NK (Shand 1943)")
)
temp1 is a list of "things" to plot on the graph; each item in the list is a list itself. Here we define the following:
- clssf,
which contains the names to use for classification purposes. Here we
indicate that the items to use for classification will be number 2 to 5
(number 1 is clssf itself), and the names to be used are, in the same
order, "metaluminous"… etc.
- lines1 to lines4 are, as the first
item of the list tell us, "lines". Their coordinates are given as a
series of x values, then y values. So the first item, lines1, goes from
(1,1) to (1,7) to (0.5,7) to (0.5,7) back to (1,1). Note that since it
will be used for classification purpose, we must close it and return to
the origin – hence the first and last points are the same. Also note
that these lines have neither color nor pattern, so they will be plotted
but invisible (you may, however, encounter them if you edit the
resulting graphical file in Illustrator or similar, as you will get
these surprising invisible polygons – yes, that's why!
- lines11
and lines12 are "abline", the first one is horizontal at y=1 (h=1),
the second one is vertical. Both are plotted in colour #2 (which has
been defined somewhere else).
- lines13 is a dashed line, defined as a "lines" Figaro object with x and y values.
- GCDkit
is a list of type "NULL" (it's not a type known to Figaro), but
contains the info GCDkit needs to add it to the menu: the name under
which it will appear (plot.name = ) and the position it will assume in
the list of graphs (plot.position=14).
Temp2 contains more template info:
temp2<-list(
text1=list("text",x=0.55,y=6,text="Metaluminous",col=plt.col[2],adj=0),
text2=list("text",x=1.05,y=6,text="Peraluminous",col=plt.col[2],adj=0),
text3=list("text",x=0.55,y=0.5,text="Peralkaline",col=plt.col[2],adj=0)
)
As you can see, these are the text labels to draw. All are of
type
"text", and come with a plotting position, a
text,
colour, and
adj
(how the text will be "aligned" around its nominal position).
Both are then merged into temp, as we discussed; and the last thing to do is to call Figaro with
sheet<<-list(demo=list(fun="plot", call=list(xlim=c(0.5,1.5),ylim=c(0,7),main=annotate("A/CNK-A/NK plot (Shand 1943)"),col="green",bg="transparent",fg="black",xlab="A/CNK",ylab="A/NK",main=""),template=temp))
(or
in fact, to put all the info in a figaro object, that is actually drawn
when calling the function
plotDiagram, which is what the GUI does).
The
Figaro object is one more list of lists, which you use to pass the
graph template (
template=temp, at the end) as well as parameters such as
the plot limits (
xlim, ylim), the title of the graph (
main) and of the
axes (
xlab, ylab), etc.
The function that has just been defined
(called
Shand(), as you may remember) does nothing visible on its own.
But it can be used as an argument to more interesting things such as
plotDiagram(diagram="Shand") # asks lots of questions, and eventually plots the diagram
or
classify(diagram="Shand") # writes the classification of each sample
and these are the functions used in the GUI.
A
non-classifying diagram is very similar – except that its template has
no
clssf object; and this allows to have open fields (remember how here
the classification fields all have to loop back to the (x,y) of the
first point…)
Our own new diagram
 |
| La/Yb vs. Yb diagram (here from Moyen and Martin 2012 |
So
– we can now use our newly acquired knowledge to define a new diagram
template. Let's go for my old friend – Martin's La/Yb vs. Yb diagram for
TTGs vs non-TTGs. As you know this diagram has been seriously misused
and I suggest you read, for instance, Moyen (2009) and Moyen &
Martin (2012) before using it, but let's consider for the time being
that you know what you are doing, right ?
As shown here, this
diagram has a closed blue "TTG" field to the left, but in other work it
is actually left open to the top. Closing it does allow to use this
diagram as a classification diagram, so let's keep it like that.
Let's create a new R file. It starts with:
LaYb<-function() # The name of the function. {
x.data<<-WR[,"Yb"]/0.22
y.data<<-(WR[,"La"]/0.33)/(WR[,"Yb"]/0.22)
We are using
normalized values here, and not ppm. Also, we need to calculate a ratio
(La/Yb)
N, that is not available by default in WR, so we do the
calculation ourselves.
Now for the template :
temp1<-list(
clssf=list("NULL",use=2:3,rcname=c("TTG/adakite","ordinary")),
only two classification items – if we are somewhere else, "classify" will return the value "unclassified".
lines1=list("lines",x=c(0.28,1.55,2.55,5.3,8.6,8.6,3,1,0.28,0.28),y=c(150,150,70,40,25,5,5,15,35,150),col=plt.col[2]),
lines2=list("lines",x=c(4.5,4.5,20,20,4.5),y=c(2,30,10,2,2),col=plt.col[2]),
The
two polygons (TTG for the first one, ordinary calc-alkaline for the
second one). The most time consuming part is actually extracting
coordinated from the diagram. Here I did it by reading the cursor
coordinates in Illustrator… Also, we need not forget to define the
colour (colour number 2 here, the common GCDkit brown by default)
otherwise the polygons will be invisible.
GCDkit=list("NULL",plot.type="binary",plot.position=100,plot.name="La/Yb vs. Yb (Martin, 1986)")
Position
in the list of available GCD graphs. I put it at 100 to be (reasonably)
sure it does not interfere with anything else. If I define more
diagrams, I'll make them 101, 102, etc.
By the way don't forget to close the parenthesis here (see the final file).
temp2<-list(
text1=list("text",x=2.5,y=100,text="TTG/adakite",col=plt.col[2],adj=0),
text2=list("text",x=14,y=22,text="ordinary",col=plt.col[2],adj=0)
)
The text annotations. A bit of trial and error is required to get good-looking (x,y) coordinates.
if(getOption("gcd.plot.text")){
temp<-c(temp1,temp2)}
else{
temp<-temp1
}
This stays the same as in Sand's example before, I like the idea of a global option
defining whether to annotate the diagram or not (plus, I think the
un-annotated versions are used in plates).
sheet<<-list(
demo=list(
fun="plot",
call=list(
xlim=c(0,24),
ylim=c(0,175),
main=annotate("La/Yb vs. Yb (Martin, 1986)"),
col="green",
bg="transparent",
fg="black",
xlab=annotate("Yb[N]"),
ylab=annotate(("(La/Yb)[N]")),
main=""
),
template=temp
)
)
And the actual Figaro call (here reformatted because I
was getting lost in all the parentheses, now you can see what's going
on !). Note the use of
annotate and square brackets to get subscripts.
Again, a closing curly brace is required here to close the function.
Using it !
Now,
you can test your diagram. The easiest is to "manually" load the
function, maybe by "sourcing" its code (from file/source R code in the
GUI, or by using
source("path\\to\\file\\LaYb.r") from the command line)
or even by copy/paste to the console. Of course your diagram will not
show in the menus (it has not been loaded yet).
But you can still get it
(and, presumably, debug it !) with
plotDiagram(diagram="LaYb")
and the classification
classify(diagram="LaYb",silent=TRUE)
(we
need
silent=TRUE otherwise GCDkit tries to look into the list of known
diagrams to retrieve its name – which, of course, fails miserably since
the diagram has not been properly loaded yet).
 |
| The "classification" menu after adding the diagram in the right directory |
When you're happy
with your diagram, it's time to drop it in the …
\Diagrams\Classification\English folder, and load some data to see what
happens (until you load data, nothing will happen, as GCDkit only
searches for diagrams when loading data !).
 |
| The diagram in GCDkit. Yes, some rocks are not adakites. Even if they have high La/Yb. | | |
Enjoy your new diagram !
The whole code in a single R file: click
here.
