Plotting single-station LOFAR data with the R programming language

The R programming language is a free  (GNU Public Licence) language and environment use for scientific computing. It has developed a strong reputation for being used for statistics and data analysis, although its use goes far beyond that. It is easily extensible with function packages and has an active community. Staff at Sodankylä Geophysical Observatory use R, particularly for signal processing and inverse problem methods.

I myself have never used R before; my native language is actually C, with PostScript, Python, FORTRAN as some of the others that I know. I've been aware of R because my colleagues use it, but other than that I've not had a chance to use it. I keep meaning to write something, but those opportunities rarely come up and, when they do, it is usually faster to hack something out in a language I know.

But a couple of days ago an good opportunity presented itself, so here's my first real programme beyond the "Hello world!".

Lassi has asked for some KAIRA data, in particular some beamlet statistics data. These are the data generated by any LOFAR single station. Not only did he want the data, but also a way to plot it... using R (which is what he uses). Fair enough. So here's my attempt:


R-code to plot a single beamlet from a LOFAR beamlet statistics (BST) file.

   # Set the name of the file that we wish to read
   filename = "20131220_150320_bst_00X.dat"

   # Set the beamlet number we wish to plot
   beamlet = 243

   # Create a connection to this file
   file = file(filename,"rb")

   # Typical maximum file size (244 beamlets x 86400 seconds/day)
   maxsize = 244 * 86400

   # Read the entire file as a raw data stream (up to maxsize)
   raw = readBin(file,double(),size=8,n=maxsize,endian="little")

   # Re-shape the data into matrix w/ 244 columns (beamlets)
   data = matrix(raw,ncol=244,byrow=TRUE)

   # Convert the data into a dB scale
   data_db = 10.0*log10(data)

   # Create an array of sample numbers for the X-axis
   sample_num = array(1:dim(data_db)[1])

   # Plot the data (we need "beamlet+1"; R counts from 1, not 0)
   plot(sample_num,data_db[,beamlet+1],
        xlab="Samples since file start time",
        ylab=expression("ADUs 10log"[10]),
        main=paste("KAIRA   ",filename,"   Beamlet=",beamlet) )

The plot looks like this:

I also plotted the same thing with KPB (KAIRA Plot Beamlet-statistcs). KPB is the python software that we normally use for plotting things.

For those curious, the command line settings for this were:

   %  kpb --db -y 243 20131220_150320_bst_00X.dat -g trace \

          -t "KAIRA   %n   bl=243" -o kpb_example.png

Of course I will continue to use Python and C code, as that is what I do. However, hopefully the above code fragment will help anyone either learning R or getting started with LOFAR single-station data.

Bear in mind that I am a novice though! So, if you're an R-expert and have suggestions on how I could improve that code, or make it more R-esque, do let us know in the comments below!


Links

• http://www.r-project.org/ -- R programming home page
• http://en.wikipedia.org/wiki/R_%28programming_language%29 -- Wikipedia