#21 – Find significant relationships in data with a CoCo Matrix

Screen Shot 2013-08-19 at 08.47.17

The CoCo Matrix (correlation coefficient matrix) is a script for R that takes a table headed with multiple variables and calculates the correlation coefficients between each of the variables, determines which are statistically significant, and represents them visually in a grid-plot. I created the CoCo Matrix to cross correlate a table with a large number of variables to quickly assess where important correlations could be found.

Screen Shot 2013-08-19 at 08.47.27

Using the CoCo Matrix

The R file can be downloaded here or copied from the textbox at the end of this post.

  1. If you know the number of samples in your dataset (n) then degrees of freedom (df) = n-2. Use this table to find the R value above which significant values lie. In the code, at the top you should change the value of “p” as per the value you just looked up. If you don’t know the value for n then run the code once and type “n” into the console.
  2. If you want, customise the colours in the customisation area of the code
  3. Run the code. A dialogue box will request a file. Alternatively replace the code to direct to the file you want to use.
  4. Voila!

This is a very rough script I wrote, and I intend to make it a lot better at some point when I have the time. If you have any suggestions  for improvements then please comment below or get in touch with me.

# CoCo Matrix version 1.0
# Written by Darren J. Wilkinson
# wilkinsondarren.wordpress.com
# d.j.wilkinson@ed.ac.uk
# The "CoCo Matrix" visualises the correlation coefficients for a given set of data.
# Like-Like correlations are given NA values (e.g. Height vs Height = NA). For the moment
# duplicates such as Height vs. Weight and Weight vs. Height remain. At some point I'll 
# provide an update that removes duplicates like that.
# Please feel free to edit the code, and if you make any improvements please let me know
# either on wilkinsondarren.wordpress.com or send me an email at d.j.wilkinson@ed.ac.uk

# Packages -------
library (cwhmisc)
library (ggplot2)
library (grid)
library (scales)
# ----------------

# Plot Customisation ----------------------------------------------------------
# (for good colour suggestions visit colourlovers.com)
col.significant = "#556270"			# Colour used for significant correlations
col.notsignificant = "lightgrey"		# Colour used for non-significant correlations
col.na = "white"						# Colour used for NA values
e1 = c("nb", "ta", "ba", "rb", "hf", "zr", "yb", "y", "th", "u")   #  p) {s = "Significant"}
		if (temp < p) {s = "Not Significant"}
		if (temp == 1) {s = NA}
		if (temp == 1) {temp = NA}
		results[h,i] = temp
		plot.data[r,4] = s
		plot.data[r,3] = temp
		plot.data[r,2] = h
		plot.data[r,1] = i


# Open new quartz window
dev.new (
	width = 12, 
	height = 9

# Plot the matrix
ggplot (data = plot.data, aes (x = x, y = y)) + 

geom_point (aes (colour = sig), size = 20) + 

scale_x_continuous (labels = e1, name = "", breaks = c(1:n.e1)) +

scale_y_continuous (labels = e1, name = "", breaks = c(1:n.e1)) +

scale_colour_manual (values = c(col.notsignificant, col.significant, col.na)) +

labs (title = "CoCo Matrix v1.0")+

theme (
	plot.title = element_text (vjust = 3, size = 20, colour = "black"), #plot title
	plot.margin = unit (c(3, 3, 3, 3), "lines"), #adjust the margins of the entire plot
	plot.background = element_rect (fill = "white", colour = "black"),
	panel.border = element_rect (colour = "black", fill = F, size = 1), #change the colour of the axes to black
	panel.grid.major = element_blank (), # remove major grid
	panel.grid.minor = element_blank (),  # remove minor grid
	panel.background = element_rect (fill = "white"), #makes the background transparent (white) NEEDED FOR INSIDE TICKS
	legend.background = element_rect (colour = "black", size = 0.5, fill = "white"),
	legend.justification = c(0, 0),
	#legend.position = c(0, 0), # put the legend INSIDE the plot area
	legend.key = element_blank (), # switch off the rectangle around symbols in the legend
	legend.box.just = "bottom",
	legend.box = "horizontal",
	legend.title = element_blank (), # switch off the legend title
	legend.text = element_text (size = 15, colour = "black"), #sets the attributes of the legend text#
	axis.title.x = element_text (vjust = -2, size = 20, colour = "black"), #change the axis title
	axis.title.y = element_text (vjust = -0.1, angle = 90, size = 20, colour = "black"), #change the axis title
	axis.text.x = element_text (size = 17, vjust = -0.25, colour = "black"), #change the axis label font attributes
	axis.text.y = element_text (size = 17, hjust = 1, colour = "black"), #change the axis label font attributes#
	axis.ticks = element_line (colour = "black", size = 0.5), #sets the thickness and colour of axis ticks
	axis.ticks.length = unit(-0.25 , "cm"), #setting a negative length plots inside, but background must be FALSE colour
	axis.ticks.margin = unit(0.5, "cm") # the margin between the ticks and the text

# Print data tables in the console

#20 – In Progress: Project Lithograph

unsorted 26I’m building. Secretly building. It’s going on behind the scenes, and until now, it’s been nowhere near reality. But now I’m getting a tad excited that one day soon I’ll be ready! I’m talking about a project that has tickled my ticker for some time…

Back when I started my degree in geology, I started to collect so-called classic texts. Looking around the charity shops of affluent parts of Edinburgh, a surprisingly populous treasure trove of old geology books awaits the saturday bargain hunter! And now, I have shelves and shelves of these books looking all… well, old! But what’s the point? Sure they’re fun to read, but some of them are so old and so few that not many people have the chance to enjoy them. And how up to date is their content? Is it going to teach you something? Sure it will!

unsorted 28

But that’s not the point! Pick up anything written by Geikie, Miller or even Holmes, and you’ll find they’re written with utter elegance. Written like the musings of Mr Darcy in a world now so alien to most of us. Modern science writing is the aspergic younger brother of classic science writing, you know… back in the days when your papers started with “On the….” and you conducted your fieldwork in a suit followed by a trail of lackeys.

But it’s not the writing style that I enjoy most. It’s the illustrations. Maybe I’ve been looking at undergraduate field notebooks for too long, but the standard of drawing in these old books is no question, art. Also, since Edinburgh and Scotland are the home of Geology, there are so many great lithographs from around this great nation that it’s all the more fantastic to see something or somewhere that you are familiar with.

unsorted 33

So, I’m planning a resource. An online gallery of these images so that all may marvel at their beauty. And who knows, you may even find some illustrations that could be made use out of. I’m slowly scanning in image after image, hand-drawn maps, double spread lithographs, meticulous creations of geologists that had the time to sit down and spend half the day sketching a geologically significant scene.

Got some old texts yourself? Would you like to get involved? Get in touch and let’s do it!

#19 Update: Should you trust your digital compass?


In the wake of my last post on digital compass clinometer functions, I was left feeling rather dissatisfied that I had found out the whole truth about the accuracy of digital compasses found in most modern smartphones and tablets. I set aside half an hour over the weekend to gather some data to see for myself what was the case.

Here I present a very rough experiment using:
1 x iPhone (built-in compass app)
2 x Analogue compass clinometers (1 x Suunto MC-2, 1 x Silva type-15)

Data Gathering
Here we are testing the ability to measure magnetic north, so any readings must reflect that. The analogue compasses used are “working” compasses so have had the baseplate correction set to zero. The iPhone compass was also set to measure magnetic north.

To ensure meaningful data, both compasses must always be measuring the same direction of travel, and must not be too close as to magnetically interfere with each other. To achieve this I used a 130 mm wide piece of plastic with parallel sides (an empty DVD dox) to separate the two devices, whilst ensuring they are pointed in the same direction. I also used the right and side of the iPhone as this side has no buttons, and thus the straight side of the device should ensure the “top” of the phone is in the correct line.

I also needed to make sure that there were no significant external magnetic interferences. Since it was pissing it down with rain outside, I opted to do it in the centre of a large room. There were no metal/electronic objects within 2m in a horizontal radius. Assuming there is no large neodymium magnet under my floorboards (?) then I also assumed there would be no interference from above or below. I was not wearing a watch, necklace, bangle/bracelet, earrings or any other piercings, or chastity belt on my person that could also interfere.


On a side note, I am very confident that compass 1 (an by association, compass 2) is indeed accurate. This is because only a matter of weeks ago I used the compass on fieldwork and on a following hiking holiday where it was used countless times for triangulation. Most of the time I was showing students how to triangulate, locating ourselves on a known location on a map. These positions were also verified by a high-accuracy mobile GPS. In short, if my compass was inaccurate or damaged, I would be well aware of the fact. Interestingly, whilst on that fieldwork a colleague of mine discovered she had fallen victim to the phenomenon whereby the magnetic polarity of your compass becomes reversed due to prolonged proximity to a smartphone.

Data was gathered in two job-lots. Using firstly compass 1 (C1, Sunnto) and then compass 2 (C2, Sylva), I sequentially adjusted the baseplate in increments of 10 degrees from 0 to 180. I had the digital and analogue compass on my lap, separated by the plastic rectangle. Sat on my swivel chair, I then shwiffled (shuffled and swivelled) round until the magnetic needle of my compass overlay the north arrow on the baseplate as well as I could get it. My lap provided much needed stability. I then recorded the reading on the iPhone as a counterpart measurement.

The Results
I must say that I was alarmed at the apparent error in the iPhone. The mean difference in reading was 4.8 degrees (s.d. 9.2) and 3.2 degree (s.d. 6.8) for compass 1 and compass 2 respectively. Interestingly, however, there seems to be a pattern in the data. Both compasses seem to show a shift in positive to negative/less positive error right around 90-110 degrees. While this pattern is more symmetrical for compass 1, both seem to have a similar negative excursion over approximately 70-80 degrees.


Theoretically, any static disturbance in the room should cause a consistent error in the measurement as neither the needle, point in the room of measurement, nor the site of any potential disturbance actually move. Even if there was something in the room causing the interference then one would expect that the analogue compass would have been affected too, and thus there would be a homogenous error in the two measurements.

So what do these errors indicate? Is this error acceptable? Of course the answer is no. Sure, the internal magnetometer in the iPhone or any other device which uses a similar piece of hardware is probably good enough for giving you a bearing on a city street, but any remotely consequential use of that information would surely lead to undesirable results. Yet, although the vast majority of people out on mountainsides around the world are unlikely to be sorts of individuals who would use an iPhone for navigation, my concern surrounds the use of “smart” devices as geological tools. Indeed in my last post I was very positive about the app GeoId. That’s not to say that the app itself is not useful. With the prospect that the app’s data-gathering capabilities is hampered by hardware accuracy, it seems like the app GeoId is probably more useful and reliable if you were to collect data with an analogue compass clinometer and then input it into the program. What a laborious and time-consuming task that would be. A task that negates the entire convenience and quirkiness of the app. If you were going to sit there for a couple of hours and input that data into something, why wouldn’t you rather input it into a piece of software such as those discussed in this blog post on structuralgeology.org.

I encourage you to gather some data yourself and let me know what you find. I’m a big advocate of using technology wherever possible, especially in geology, so I’m really interested in knowing more about the limitations of such technology.

# 18 – Quick Review of GeoID and Digital Compass Clinometers


GeoID is a smart device app (available in iTunes, Google play, and CNET) which is a powerful tool for the field geologist, from student to experienced professional. It combines a variety of useful functions from a compass clinometer to plotting and analysis of 3D planes and poles. In this post I’m going to review some key features of the app, after having field tested it over two field trips in the North West Highlands of Scotland.

Background on the App


GeoID was developed by Jin Son, an independent developer and Energy Systems Engineering student based at Seoul National University. I’m not sure how long ago he developed the app, but it has to have been around more more than 12 months. It has not been rated at all in iTunes, and has had only a few ratings in Google’s App Store (4/5 stars).

Compass Clinometer Feature
There are already a range of apps out there that can effectively replace your traditional compass clinometer (e.g. GeoCompass, Strike and Dip), but none I have seen have the utility and design to match GeoID.

But the question to ask here is not which app looks better, but rather is a compass clinometer on your phone/tablet more accurate than a traditional compass-clino? Well the accuracy, and indeed precision, of the digital compass and digital clinometer is the product of two things: the hardware, and the coding. The compass part of the software uses both the in-built magnetometer AND the accelerometer. The magnetometer measures a heading from the Earth’s magnetic field, and combines this with data from the accelerometer on the orientation of the phone, to give you a direction of travel on your screen. Most compass applications on smart devices provide up-to-date corrections for magnetic vs. true north, so this shouldn’t be an issue with well coded apps. The general, better-informed, consensus in many online forums is that the accuracy of newer digital compasses is negligible compared with traditional needle compasses. This is based on direct comparisons between digital and analogue compasses reported online, and conducted by myself.

The clinometer seems to be more precise than analogue clinometers in most compass-clino’s. In a blog post I found, the iPhone’s clinometer was compared, by a carpenter, with a tilt box (a common, high precision clinometer used in many professional workshops). He wrote:

“Tried Clinometer last night and compared it to my tilt box, its most likely as accurate as the tiltbox – I used it as a reference to find out if two parts of a railing I am building were co-planar, the difference between the two planes was within 0.1 degrees between the tiltbox and the iphone app. On an absolute scale (i.e. I didn’t zero the reading on the first plane before moving it to the second plane), the two readings were within 0.1 degrees of each other, that’s probably as accurate as it gets in my shop.”

So it seems that the compass clinometer on a smartphone may well be a healthy rival to the conventional analogue device. Its power multiplied by the ease and speed with which it can make readings (which is not always a benefit, as I’ll discuss later).

GeoID further increases the power through clever coding. Dip/strike or Dip/Direction readings can be averaged over a given time interval to make readings more stable and precise. So when recording a surface there is small delay before a “stable” reading can be taken. If you record a reading before it is stable you are warned with a sound.


Data recording & real time plotting
So, the compass clinometer is a quick, accurate and powerful tool. On top of that, in GeoID the compass clinometer button allows you to quickly gather and record lots of structural data on poles and planes, all whilst seeing it plotted in real time onto a stereonet of your choice. Plus, if to have GPS in your device, each datapoint is geotagged!


I tested the app out on a thrust-related fold in Durness limestone on the shore of Loch Assynt, North West Highlands of Scotland. In under an hour I was able to collect 300 dip and dip-direction measurements. Although a really nice set of data, I noticed that the ability to gather such a large volume of data must be used with care. The limestone here is heavily carstified, and so measuring a large number of surfaces resulted in quite a large scatter on each of the fold limbs. Arguably, had GeoID an in-built averaging function this would have been less of a problem, but alas it was.

Once the data is collected, you can see the stereonet full-screen and change from equal-area to equal-angle projections. One thing that would be really useful here would be the ability to have layers of data so that planes may be fitted to each of the limbs to come up with a mean plane per limb, such that you could read off things such as the mean plunge and plunge direction.

I had also given my iPhone to another geologist to gather some extra data. Although smaller, it could be placed on a clipboard to average the bedding plane surfaces. The app allows you to then share the data via email, bluetooth, or export to file, however there is no ability to collaborate on the same project, so I was thus unable to amalgamate the two datasets.

I really like the GeoID app in its simplest form, as a tool for gathering dip/direction data quickly and accurately. I think that it has a lot of potential on the analytical side of things, although I appreciate the complexity in implementing such features into the program. What I think this app needs is more users, more rafters, and for us to give some constructive feedback to Jin Son on how he can make this app even more awesome!

– Compass clinometer is incredibly useful and accurate
– Ability to quickly gather and store pole/plane structural data
– Real-time plotting of data allows on-the-fly interpretation
– Very user friendly
– Easily send/export data

– Limited to one layer of data
– No structural analysis relevant to folding etc
– Cannot collaborate on, or amalgamate projects

#17 Homage to the Sausage – Be the Boudin


Geologist Tom Challands and his unique interpretation of “be the boudin”.

In the North West Highlands of Scotland, there’s a stretch of coastline graced with the, often understated, natural beauty of creamy soft sands and tempting turquoise water, and bedrock that tells a story far more captivating than man has ever written for himself.

Much of the lowlands here are made of grey metamorphic “Lewisian” gneisses that record some of the earliest parts of Earth’s history available in the rock record. A world so old and foreign that it’s amazing we can tell anything about it at all.


Suilven mountain, from the road to Clachtoll. Sulven is made of near horizontally bedded arkosic sandstones, deposited onto the hummocky lowlands of Lewisian gneiss.

A common feature of the gneisses are textures that tell a story of stretching in conditions deep and hot enough to reach the point at which rock deforms like soft toffee. Harder, darker layers made of minerals rich in iron and magnesium minerals (a.k.a. mafic) respond differently to deformation compared to the lighter layers made of minerals richer in silica (a.k.a. felsic). The contrast in response to extensional deformation result in the breaking up of the darker layers into chains of sausages (aka. boudinage). Each individual block, often lensoidal in shape, is known as a boudin (i.e. sausage).

These boudinaged dark layers in the gneisses are common not just here, but in geologically similar areas around the world, and can be seen on scales of a few millimetres to many tens of meters.


Boudinaged layers in gneiss. Note the darker layers have begun to break up into blocks in a disarticulated chain.

Boudinaged mafic layers in gneiss, Kangerlussuaq, Greenland.

Boudinaged mafic layers in gneiss, Kangerlussuaq, Greenland.

At the beach of Achmelvich (58.172915, -5.302604, [map]), arguably one of Britain’s most stunning beaches, you find gneisses which show clear evidence of the type of deformation mentioned above. There, the gneisses are folded, re-folded, and folded again into a complex, contorted basement which has been the subject of intense unravelling by geologists for over a hundred years.


Achmelvich Beach (58.172915 N, -5.302604 W), arguably one of Britain’s most beautiful beaches.

As you turn you back to the  shining blue sea and the clean cream sands, you find a crag with a whole set of these boudinaged dark layers. The brunt of the sea has, over thousands of years, gouged and gouged at this face. The pummelling and battering of the crag has plucked out one of these boudins, a pip of darker rock, leaving a pseudomorph of the former rocky sausage.

The half-egg shaped recession in the cliff is actually a surprisingly perfect perch for the derrière of most beach dwellers or geological guests, and has over the years been used for such a purpose.


Myself having a rather traumatic boudinage birthing experience.

Myself having a rather traumatic boudinage birthing experience.

Geologist Madeleine Berg during her first boudinisation.

Geologist Madeleine Berg during her first boudinisation.

But alas, the boudin is nowhere to be seen, and the perch is there like a memorial to the former knocker which ironically got knocked out. Almost in tribute to the boudin, geologists visiting the beach began paying their respects by assuming the form of the boudin for a brief, but immensely symbolic moment. For not only is this a sign of geological respect, homage to the “saus-age”, it has become an important ritual of initiation which marks a milestone in ones geological career. A milestone which will divide ones life into pre- “birth by boudin” (BBB) and post-BBB.

The first geologists to “be the boudin” were Nigel Kelly and Simon Harley. Simon is a regular visitor to the beach and overseas and encourages in the ritual with new geologists every year.

So, should you wonder out yonder, please have a ponder, on the act of birth by a boudin. And when you have, drop me a line and a photo so that we may appreciate your “be the boudin”!

#16 Colour: A quick guide to its use in informative graphics

1.0 Introduction
The most fundamental employment of colour in [qualitatively/quantitatively] informative graphics is to allow the observer to easily distinguish elements of the information displayed. The strict definition of the term “colour” must be cast aside here, as just as useful are black, greys and white. Much of what will be discussed in this post may seem intuitive, but it is the guidelines and their employment in practice which unfortunately escapes the many.

The most important point to take away from this post is: colour used well can enhance and clarify visual information; and colour used badly will likely obscure and confuse.

2.0 Principals of Colour in design

2.1 Hue, value and chroma

Colour designers use three variables to describe any particular colour:

  1. Hue – Fig. 2.1 – The name of the colour (e.g. Red, or Blue)
  2. Value – Fig. 2.2 – The lightness or darkness of that colour
  3. Chroma – Fig. 2.3 – Equivalent to saturation (reducing to zero gives equivalent grey value)

2.1 – Hue values

2.2 – Value, or brightness.


2.3 – Chroma, or saturation.

2.2 Legibility

If something is legibleit is clear to read. Although this extends to literal clarity, here we’re talking about optical clarity to the human observer. Hue and chroma do not contribute much to legibility, it is the luminance contrast (or contrast in value) of the background and the foreground. Fig. 2.4 shows varying degrees of legibility, due to differences in contrast between the red and black text, and the background.


2.4 – Legibility is a result of the difference in value (brightness).

Legibility, therefore can be achieved in a number of ways. Adequate contrast may be achieved in a either a monochromatic or multi-chromatic image. As far as monochromatic images are concerned (Fig. 2.5), the case is simple and value contrast is the only thing one needs to consider. For colour images (e.g. Fig. 2.7), in the strictest sense of the word, contrast surely gives us legibility, but unlike in monochromatic images, legibility is then not our only concern!


2.5 – Monochromatic image achieving legibility through contrast in greys.


2.6 – Image shows high analogy of colours and low contrast.


2.7 – Image shows high contrast and two analogous colour groups.

3.0 Selecting a colour palette

3.1 Review of concepts

So by now you should understand the important colour principals should be employed in the following manner:

  • Select COLOUR to suite FUNCTION
    e.g. I’ll select a red hue for high temperatures and a blue hue for low ones. Intuitive, eh?

3.1 – Colours here have been selected for intuitive function.

  • Use colour ANALOGY to delineate GROUPS
    e.g. Set A = Orange, Set B = Green. Each set may have different chroma values.

3.2 – Here, the red and blue palettes have analogous variations of chroma (i.e. saturation).

    i.e. New data vs. old data.

3.3 – Contrast here is used to highlight the “new” data over the “old” data.

    Often useful when labeling
Reducing the diaphragm on the microscope gives better colour reproduction for some camera phones. It also increases the apparent relief od the minerals.

3.4 – White text on a black background gives high contrast, which translates into high legibility.

In most situations where information is being displayed, it is best to limit your colour palette to two or three different hues. Using a limited hue palette, you are forced to increase variation by changing the value and chroma of your core hues. Examples showing theses limited palettes have already been discussed.

3.2 Colour Palettes vs. Symbols

Of course a lot of researchers need to display information in monochromatic images, due to constraints applied by the journal in which they wish to publish the data. The de facto choice for displaying multiple sets of information is to use different symbology (e.g. dashed/dotted lines, different point symbols). However, studies have shown that the human eye is only capable of understanding a maximum of 7 variables in any given set. This means that if you have more than seven symbols or line types, your plot is going to be hard to decipher. As discussed, the same applies to colour. So when it comes to symbols, it is also advisable to limit the number used to around three. This thinking is employed in the much used GGPLOT package for R. The package, written by Hadley Wickham limits the automatically assigned symbols to a plot (PCH values) to a maximum of seven.

3.3 Complex palettes

The idea of restricting your palette to only three hue values, or indeed the number of symbols used, is an ideal one. However, there are often many situations when you simply need to use more. More symbols or more colour? Well the advice from the experts is to use colour. It is much less work, and the mind is much better at distinguishing colours rather than a series of [more complex] shapes. So should you need to display 7 sets of data then colour is the way to go. Consensus on the ground is, that should you need to display more than 7 sets of data, then you need to rethink the way you are visualising your information.

3.4 Colour palette resources

Now you understand a lot more about using your colours, you could go out there and compose your own colour palettes… OR you could use any number of great websites to provide the palettes for you! AND if you use GGPLOT, the first of these (ColorBrewer) can be integrated with your plots! Remember, these palettes are not only great for plotting data, they are also great for providing colour schemes for posers and presentations!

4.0 A word on background colours

Some people prefer using background colours to their plots, commonly pale yellows or greys. The problem is that most colours (particularly preset colour palettes) are chosen to be printed on white paper, so it stands to reason that they should be displayed digitally on a white background too. Not only this, our brain is designed to constantly normalise our palette to a local reference for white. For example, if you’ve worn yellow ski goggles all day, when you take them off everything appears blue. This is also why when you take a photo indoors under a tungsten bulb, the image often comes out orange, more orange than you see it. This is because the camera’s white balance has not been set to match your eyes.

3.5 – Most camera users encounter white balance issues when shooting indoors. This is because the cameras white value is not set correctly. In the human brain, white balance is constantly being corrected so we would see the right hand image.

3.6 – Early computer screens didn’t have enough brightness to show black text on a white background. So the text had to be white and the background black to increase legibility.

Dark backgrounds are only ever recommended when the information is being displayed in a dark environment. The use of light text on a dark background reduces problems of legibility, allowing the viewer to see the information without being distracted by the background. The reason that early computer screens and slides from the earliest projection systems had dark backgrounds was because the displays weren’t that bright, and consequently needed to be used in rooms with reduced lighting. Now, however, projectors and screens are bright enough that this isn’t the case.

#15 Alkali Silica Template

Alkali Silica

Does what it says on the tin.


#-------- INFORMATION ---------
# Plotting points from Hugh
# Rallinson's "Using Geochemical
# Data" book. Code compiled by
# Darren J. Wilkinson,
# Grant Inst. Earth Science
# The University of Edinburgh
# d.j.wilkinson@ed.ac.uk

# -------- CONTROLS ----------
y.max = 16
x.min = 35
x.max = 80
lab.size = 5
save = "/Users/s0679701/Desktop/"
filename = "test.png"

# -------- LIBRARIES ----------
library (grid)
library (scales)
library (ggplot2)

# -------- DON'T EDIT ----------
a = c(1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 8, 8)
x = c(41.0,41.0,52.5,50.0,52.5,57.6,63.0,63.0,45.0,45.0,61.0,63.5,45.0,52.0,57.0,63.0,69.0,45.0,49.4,52.0,52.0, 48.4,53.0,57.0,57.0,69.0,69.0,76.6,41.0,45.0)
y = c(0.500,7.000,14.000,15.126,14.000,11.700,7.000,0.500,0.500,5.000,13.500,14.830,5.000,5.000,5.900, 7.000,8.000,9.400,7.300,5.000,0.500,11.500,9.300,5.900,0.500,12.500,8.000,0.500,3.000,3.000)
lines = data.frame (a, x, y)

b = c("Picro-Basalt", "Basalt", "Basaltic Andesite", "Andesite", "Dacite", "Basanite", "Trachy Basalt", "Basaltic Trachyandesite", "Trachyandesite", "Trachydacite", "Rhyolite", "Tephrite", "Phonotephrite", "Tephriphonolite", "Trachyte", "Foidite", "Phonolite")
c = c(1:17)
x = c(43, 48.5, 54.5, 60, 68, 43, 48.75, 53, 57.5, 65, 75, 46, 49, 53, 65, 45, 58)
y = c(2, 2.5, 3, 3.5, 4, 6, 5.5, 6.5, 8.5, 9, 8, 8, 9.5, 11.5, 12, 13, 14)
labels = data.frame (b, x, y)

x = c(39.8,65.5) #MacDonald (1968)
y = c(0.35, 9.7) #MacDonald (1968)
sub = data.frame (x, y)

# -------- BEGIN PLOT ----------

ggplot (lines, aes (x=x, y=y)) +

# Field Boundaries
geom_line (
aes(linewidth = factor (a))
) +

# Alkaline-Tholeiitic Line
geom_line (
data = sub,
aes (x=x, y=y),
linetype = "longdash"
) +

# Field Labels
geom_text (
data = labels,
aes(x = x, y = y, label = b),
size = lab.size
) +

scale_x_continuous (
name = (expression(paste("SiO"["2"], " (wt. %)"))),
limits = c(x.min, x.max)
) +

scale_y_continuous (
name = (expression(paste("Na"["2"], "O + K"["2"], "O", " (wt. %)"))),
breaks = c(seq(0, y.max, 2)),
limits = c(0, y.max)
) +

theme (
plot.title = element_text (vjust = 3, size = 20), #plot title
plot.margin = unit (c(3, 3, 3, 3), "lines"), #adjust the margins of the entire plot
panel.border = element_rect (colour = "black", fill = F, size = 2), #change the colour of the axes to black
panel.grid.major = element_blank (), # remove major grid
panel.grid.minor = element_blank (), # remove minor grid
panel.background = element_rect (fill = "white"), #makes the background transparent (white) NEEDED FOR INSIDE TICKS
legend.background = element_rect (fill = "white"),
legend.justification=c(1, 1),
legend.position = c(1, 1), # put the legend INSIDE the plot area
legend.key = element_blank (), # switch off the rectangle around symbols in the legend
legend.title = element_blank (), # switch off the legend title
legend.text = element_text (size = 15), #sets the attributes of the legend text
axis.title.x = element_text (vjust = -2, size = 20), #change the axis title
axis.title.y = element_text (vjust = -0.1, angle = 90, size = 20), #change the axis title
axis.text.x = element_text (size = 17, vjust = -0.25, colour = "black"), #change the axis label font attributes
axis.text.y = element_text (size = 17, hjust = 1, colour = "black"), #change the axis label font attributes
axis.ticks = element_line (colour = "black", size = 0.5), #sets the thickness and colour of axis ticks
axis.ticks.length = unit(-0.25 , "cm"), #setting a negative length plots inside, but background must be FALSE colour
axis.ticks.margin = unit(0.5, "cm") # the margin between the ticks and the text

ggsave (paste (save, filename), height = 12, width = 18, dpi = 75)