ProTrackR: chiptunes in R (part two)

A short follow-up on my previous post. I have just released a new version of the ProTrackR package (available on the CRAN). It can now actually play ProTracker modules, yeah! The routine isn't perfect yet but it does a decent job at the moment. The example below shows how you can download a module from ModArchive and play it:

	# load the required packages
	require(ProTrackR)
	# set up a connection to 'cyberrid.mod':
	con <- url("http://api.modarchive.org/downloads.php?moduleid=40293", "rb")
	# download module:
	mod <- read.module(con)
	close(con)
	rm(con)
	# and play it:
	mod.wave <- playMod(mod,
	max.duration = 300,
	target.rate = 22050)
	# plot the waveform:
	plot(mod.wave)

view raw ProTrackR02.01.r hosted with ❤ by GitHub

ProTrackR: chiptunes in R (part one)

ProTracker is a popular music tracker to sequence music on a Commodore Amiga machine. As I still own a Commodore Amiga 500, I thought it would be a nice challenge to develop an R package that could handle ProTracker module files, and so I did. The first version of this package was just released on CRAN.

In this post I will briefly explain what ProTracker modules are and how they work. Then I will give a short demonstrations of the current sweetness the ProTrackR package has to offer. On this blog I will try to keep you updated on new versions of this package and its new features.

The module (MOD) is a computer file format used primarily to represent music. A MOD file contains a set of instruments in the form of samples, a number of patterns indicating how (which notes and effects) and when the samples are to be played, and a list of what patterns to play in what order.

The first release of the ProTrackR package contains the bare bones required to import and export such module files, and modify them any way you'd like (and the file format restrictions allow you to). This is a good time to start with an example. Let's start by loading the module 'Elekfunk' from The Mod Archive.

	# required packages:
	require(plot3D)
	require(ProTrackR)
	# Read the module 'elekfunk' from the Mod Archive:
	con <- url("http://api.modarchive.org/downloads.php?moduleid=41528", "rb")
	mod <- read.module(con)
	close(con)
	rm(con)
	# see how much memory the module object is eating:
	format(object.size(mod), "Kb")
	# calculate the ratio between the size of the object
	# in memory and the size of the mod.file:
	print(as.numeric(object.size(mod)/moduleSize(mod)))
	# note that the object is roughly four times the size as the original file.
	# This is probably because the Wave object (which the ProTrackerSample objec
	# enherits), stores 8 bit wave data as integer values. In R, integers are
	# actually 32-bit long, which explains the factor (32 / 8 =) 4 difference.

view raw ProTrackR01.01.r hosted with ❤ by GitHub

This module was used in a demo by Sanity to show off the Amiga's capabilities, which were impressive ate the time.

Playing routines in the ProTrackR package are currently limited to just playing the samples at specific notes (more fancy routines are in development). The examples below show how to play all samples in a module, a specific sample at different notes and applying different fine tunes. But at least you can export them as wave files.

	# play all samples in the module at note "C-3"
	playSample(mod)
	# play chromatic scales of octaves 2 and 3 using the 23rd sample
	data(period_table)
	for (i_note in as.vector(outer(names(period_table)[-1:-2],
	c(2:3), function(x, y) paste(x, y, sep = ""))))
	{
	print(i_note)
	playSample(PTSample(mod, 23), 0.1, T, i_note, 0.5)
	}
	# play all finetunes of C-3
	for (ft in -8:7)
	{
	print(ft)
	playSample(PTSample(mod, 23), 0.1, T, "C-3", 0.5, finetune = ft)
	}
	# export one of the samples as *.wav file
	writeWave(PTSample(mod, 25), "guitar slide.wav")

view raw ProTrackR01.02.r hosted with ❤ by GitHub

Combined with the power of R, all kinds of fancy analyses become possible. The module for instance contains a sample of a guitar slide. In R we can now generate a power spectrum (which would have been tricky to do real time on an original chipset Amiga).

	# select a sample:
	plot_samp <- PTSample(mod, 25)
	# generate power spectrum
	spec <- powspec(waveform(plot_samp),
	wintime = 0.1,
	steptime = 0.005)
	spec_cutoff <- quantile(log10(spec), 0.01)
	spec[log10(spec) < spec_cutoff] <- 10^spec_cutoff
	# fancy 3D plot of waveform and powerspectrum
	png("powerspec.png", 600, 300, type = "cairo",
	bg = "transparent", pointsize = 14)
	par(oma = c(0,0,0,0), mar = c(0.5,0.5,0.5,0.5))
	lines3D(x = 60*(0:(length(plot_samp) - 1))/(length(plot_samp) - 1),
	y = rep(15, length(plot_samp)),
	z = 3 + 11*waveform(plot_samp)/255,
	phi = 40, theta = 21, lwd = 0.15,
	xlim = c(0,60), ylim = c(0,15), zlim = c(-3, 14.5),
	xlab = "time", ylab = "Frequency", zlab = "Amplitude",
	col = "black", lwd = 0.5, scale = F, bty = "b2")
	persp3D(x = 60*((1:ncol(spec)) - 1)/(ncol(spec) - 1),
	y = 15*((1:(nrow(spec) - 3)) - 1)/(nrow(spec) - 4),
	z = t(log10(spec))[,-(1:3)],
	xlab = "time", ylab = "Frequency", zlab = "Amplitude",
	add = T, colkey = F, shade = 0.1)
	dev.off()

view raw ProTrackR01.03.r hosted with ❤ by GitHub

The power spectrum shows the shift in frequency which you would expect for a guitar slide. There are plenty other examples I could show you here but leave it at this for the moment. Hopefully you enjoyed these examples of the ProTrackR package, please stay tuned for updates on this package.

Geotags from JPEG images

Most modern camera's and telephones can geotag their pictures. Kind of neat that you can use your camera as GPS logger. Even better if we could do something useful with it in R. In this post I will show how GPS coordinates can be extracted from jpeg images and utilized in R. Although it sounds pretty simple, and it is, it is not very straight forward to achieve.

In this example I will use pictures I made on a day out. We made a small round trip by steamtram from Hoorn to Medemblik where we visited castle Radboud. I will show how to plot the locations where I took the pictures can be plotted onto a map. For privacy reasons I will not share the pictures. However, the example should also work with your own geotagged pictures.

	# required packages
	require(rgdal)
	require(OpenStreetMap)
	# select files for which geotags need to
	# be extracted and store filenames in data.frame
	file_info <- data.frame(file.name = choose.files(filters = Filters["jpeg",]))

view raw geotag01.r hosted with ❤ by GitHub

Geotags and other meta-information is stored in jpeg files in the so-called exchangeable image file format (or EXIF). Luckily it can be easily obtained using the function GDALinfo. This is where things get messy. The function returns the EXIF codes as character strings in the attributes of the object that is returned. So, we need to fetch the relevant attributes; find the relevant EXIF codes; and extract the information from the code and convert it into something workable.

The EXIF codes for GPS information will look something like this ‘EXIF_GPSLatitude=(52) (40) (57.45)’. I use regular expressions to extract the coordinate information from between the brackets. I'm no expert in the field of regular expressions, so there may be more efficient regular expressions to get the right information. If you have improvements, please feel free to leave them in the comments.

I've wrapped the routines to isolate the respective EXIF code and to extract the numerical coordinates from the literal EXIF character string in slick functions (‘get_EXIF_value’ for getting the literal string and ‘get_EXIF_coordinate’ to get the numerical coordinate from that string).

	# function to get a parameter value from EXIF code,
	# x = a vector of character strings representing the EXIF codes of a file
	# parameter = character string representing the parameter
	# for which you want to extract the value
	# the literal character string is returned for the requested parameter
	get_EXIF_value <- function(x, parameter)
	{
	line <- x[grepl(paste("EXIF_", parameter, "=", sep = ""), x)]
	if (!is.null(line) && length(line) > 0)
	{
	pos <- gregexpr("=", line, fixed = T)[[1]]
	pos <- pos[length(pos)]
	line <- substr(line, pos + 1, nchar(line))
	line <- gsub("[(]|[)]", "", unlist(strsplit(line, "([)]\\s+[(])|[)][()]")))
	}
	return(line)
	}
	# function to get the longitude and latitude from the
	# character string obtained with the function above
	# x = a character string with the coordinates
	get_EXIF_coordinate <- function(x)
	{
	lat <- get_EXIF_value(x, "GPSLatitude")
	lat <- as.numeric(lat)
	lat <- lat[1] + lat[2]/60 + lat[3]/(60*60)
	latref <- get_EXIF_value(x, "GPSLatitudeRef")
	lat <- ifelse(toupper(latref) == "S", -lat, lat)
	if (length(lat) == 0) lat <- NA
	lon <- get_EXIF_value(x, "GPSLongitude")
	lon <- as.numeric(lon)
	lon <- lon[1] + lon[2]/60 + lon[3]/(60*60)
	lonref <- get_EXIF_value(x, "GPSLongitudeRef")
	lon <- ifelse(toupper(lonref) == "W", -lon, lon)
	if (length(lon) == 0) lon <- NA
	return (c(lon, lat))
	}
	# for each file obtain the EXIF codes and extract the coordinates
	GPStags <- apply(array(file_info$file.name), 1, function(x) {
	EXIF_code <- NULL
	try (EXIF_code <- attr(GDALinfo(x), "mdata"), silent = T)
	if (is.null(EXIF_code)) return(c(NA, NA)) else return(get_EXIF_coordinate(EXIF_code))
	})

view raw geotag02.r hosted with ❤ by GitHub

After some rearranging of the data, we have a data.frame holding the filename and the longitude and latitude in decimal WGS84 degrees. I noticed that not all my pictures had a geotag (even though I had geotagging turned on, on my device). Apparently taking pictures is not a full proof method for logging locations (at least not on my device).

	# rearrang the data
	GPStags <- t(matrix(unlist(GPStags), 2))
	# store in the data.frame
	file_info$lon <- GPStags[,1]
	file_info$lat <- GPStags[,2]
	# check if any files lack GPS longitudes
	table(is.na(file_info$lon))
	# clean up after us
	rm(GPStags)

view raw geotag03.r hosted with ❤ by GitHub

Let's plot the coordinates on a map. I download some OpenStreetMap-based map tiles using the function ‘openmap’. This function conveniently allows you to specify the upperleft and lowerright coordinates of your map. That way, you can be certain that your locations will fit on the map.

As we cannot change the projection of the downloaded maptiles. We therefore have to transform the coordinates of our locations to the same projection as the maptiles. We do this by first creating a SpatialPoints object from our coordinates and then transform those coordinates to the proper projection by applying spTransform. Plotting is now simply a matter of calling some plot routines.

	# Download nice-looking map-tiles for the relevant location
	photo_map <- openmap(upperLeft = c(max(file_info$lat, na.rm = T) + 0.2, min(file_info$lon, na.rm = T) - 0.2),
	lowerRight = c(min(file_info$lat, na.rm = T) - 0.2, max(file_info$lon, na.rm = T) + 0.2),
	type = "stamen-watercolor")
	# we need to turn our coordinates into an sp object in order to
	# properly project the points onto the map.
	photo_points <- SpatialPoints(na.omit(file_info[,c("lon", "lat")]),
	CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
	# transform the lon lat coordinates to the same projection as the map tiles.
	photo_points <- spTransform(photo_points, photo_map$tiles[[1]]$projection)
	# plot the map and locations of the pictures
	png("jpg_geotags01.png", photo_map$tiles[[1]]$yres/3, photo_map$tiles[[1]]$xres/3, type = "cairo")
	plot(photo_map)
	plot(photo_points, add = T, pch = 19, cex = 0.8)
	plot(photo_points, add = T, pch = 19, cex = 0.5, col = rainbow(length(photo_points)))
	dev.off()

view raw geotag04.r hosted with ❤ by GitHub

To make things a bit more fancy, I used the rainbow palette to color the locations in the order in which the pictures were taken. Note that the use of geotagged images is limited. Proper GPS loggers for instance also log the precision of the position and your speed. But I will save that for a future post.

R Chemistry in 3D

The three dimensional model of the acetylcholine molecule shown below is interactive, when your browser supports WebGL. You can use your mouse to rotate and zoom the model. In this post I will show you the main steps required to create such an interactive 3D model in R: easy peasy (more or less).

Your browser does not support the HTML5 canvas element.

Drag mouse to rotate model. Use mouse wheel or middle button to zoom it.

There are several R packages available for cheminformatics. However, as far as I could tell, there aren't any packages that are directly able to create 3D visualisations of molecules. Therefore, we need the capabilities of the rcdk package to load molecule files and 3D rendering capabilities of the rgl package.

Although it is possible to derive a 3D structure from a SMILES code (see e.g., the SMILES translator from the Computer-Aided Drug Design Group), I couldn't find such functionality in R packages. Therefore, we use an existing 3D model of the acetylcholine molecule from PubChem in this example. The first part of the script downloads the molecule and then loads it into the working memory.

	# required packages
	require(rcdk)
	require(rgl)
	# download 3D structure of acetylcholine from Pubchem:
	download.file("https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/cid/187/record/SDF/?record_type=3d&response_type=save&response_basename=Structure3D_CID_187",
	"acetylcholine.sdf", mode = "wb")
	# load the downloaded file:
	molfile <- load.molecules("acetylcholine.sdf")[[1]]

view raw molecule3d01.r hosted with ❤ by GitHub

Subsequently, we define some atom properties, which we will need for plotting later on: the atom symbol, the colour to be used in the model, and the Van der Waals radii of the atoms as listed at Wikipedia. Note that we only define these properties for the four atoms present in acetylcholine. If you want to apply the code to other molecules, you have to make sure all relevant atoms are specified. In our example we continue by getting the 3D atom coordinates from the molecule file, sticking it in a data.frame and merge it with the previously defined atom properties.

	# create a data.frame with atom properties
	atom.props <- data.frame(
	symbol = c("O", "N", "C", "H"),
	vanderwaals = c(1.52, 1.55, 1.7, 1.2), #Van der Waals radii of atoms
	colour = c("red", "blue", grey(0.1), "white")
	)
	# get atoms from the molecule
	atoms <- get.atoms(molfile)
	# get 3D locations of atoms as read from the file
	atom_dat <- data.frame(do.call('rbind', lapply(atoms, get.point3d)))
	names(atom_dat) <- c("x", "y", "z")
	# merge with atom info
	atom_dat$symbol <- unlist(lapply(atoms, get.symbol))
	atom_dat <- merge(atom_dat, atom.props, all.x = T, all.y = F)

view raw molecule3d02.r hosted with ❤ by GitHub

Now, all we have to do is to plot spheres at the positions of the atoms, with the right sizes and colours. For this we use the spheres3d function. We could also define our own (smoother) mesh of a sphere, but that will result in huge html files when saving the model with writeWebGL. We therefore stick with the lighter spheres3d function. Unfortunately, in Internet Explorer strange shadows appear (a bug?), but in other browsers it looks fine. Just see for yourself at the top of this page.

	# simply plot the atoms as coloured spheres.
	# multiply the Van der Waals radii with 0.75
	# to reduce overlap of atoms.
	mapply(function(x, y, z, vdw, col){
	spheres3d(x, y, z, + vdw*0.75, col = col, add = T)
	}, atom_dat$x, atom_dat$y, atom_dat$z,
	atom_dat$vanderwaals, atom_dat$colour)
	# Save as an interactive model on web page
	writeWebGL()
	# Open in default browser:
	browseURL(paste("file://", file.path(getwd(), "webGL/index.html"), sep = ""))
	# Tested in Chrome on Windows 7 machine with success
	# Tested in IE on Windows 7 machine with limited success
	# Tested in Chrome on Android device without success

view raw molecule3d03.r hosted with ❤ by GitHub

Although writeWebGL produces an html file, which is ready to go and can be uploaded to your website directly, it may need some editing when embedding it in a custom template. This requires some knowledge of the html and JavaScript languages. I won't be explaining the modifications I made here.

As a bonus, if your more a fan of the ball and stick model, below the code to add sticks to your molecule model.

	#################################
	# ball and stick model
	#################################
	# get bond info
	bonds <- get.bonds(molfile)
	bonds.atoms <- data.frame(do.call('rbind', lapply(bonds, function(x){
	atoms <- get.atoms(x)
	start <- get.point3d(atoms[[1]])
	end <- get.point3d(atoms[[2]])
	return (list(start.x = start[1], start.y = start[2], start.z = start[3],
	end.x = end[1], end.y = end[2], end.z = end[3]))
	})))
	# open rgl device
	open3d()
	# plot sticks
	mapply(function(x1, y1, z1, x2, y2, z2){
	shade3d(cylinder3d(data.frame(x = c(x1, x2), y = c(y1, y2), z = c(z1, z2)),
	radius = 0.1, sides = 16,
	e2=rbind(c(1,0,0),c(1,0,0))), col="grey")
	}, bonds.atoms$start.x, bonds.atoms$start.y, bonds.atoms$start.z,
	bonds.atoms$end.x, bonds.atoms$end.y, bonds.atoms$end.z)
	# plot balls
	mapply(function(x, y, z, vdw, col){
	spheres3d(x, y, z, vdw*0.25, col = col, add = T)
	}, atom_dat$x, atom_dat$y, atom_dat$z,
	atom_dat$vanderwaals, atom_dat$colour)

view raw molecule3d04.r hosted with ❤ by GitHub