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

North Celestial Pole

In a LOFAR system, a "beamlet" is a pointing direction and a frequency channel (which is called a subband). In some ways, you can think of it as a  "radio pixel" on the sky. These can be static, pointing at a fixed direction, or tracking, where they follow a particular astronomical object as it rises, travels across the sky, and then sets.

This rising and setting, and the motion of the stars generally, is dominated by the rotation of the Earth. Although some objects move about in the sky, due to orbits or proper motion, it is the 24-hour cycle of the Earth that is the dominant effect.

However, there are two points on the sky which do not move very quickly nor very much (although they do still move!). These are the Celestial Poles. These are the points at infinity that are projections of the Earth's rotation axis. Check the Wikipedia article for a full description.

So, in principle, if you were to point a beamlet at the North Celestial Pole (which is the one we see from KAIRA), then you would expect the power level measured to be absolutely flat. There is radio power coming from that direction, but it should remain constant, as it is only rotating within that radio pixel.

It is sort of flat... but not quite!

Here's an example:

The variations come from several things. Tall vertical spikes are radio-frequency interference. And the fuzzy-noisy patches can be from scintillating sources going through the side lobes.

More generally, gradual variations can be from larger areas of Galactic emission and discrete sources passing through the sidelobes of the beam pattern. And then there are subtle variations caused by snow on the antennas, temperature variations, receiver stability, etc.. These are much smaller, but are still present.

As a result, it can be quite a challenge to calibrate data from the instrument.

Such is the fun of radio astronomy!  :-)

Receiver/USRP installation on the Kiruna Demonstrator Array

This week a team from KAIRA/SGO and LTU, together with our colleagues from EISCAT, are working to install a digital beamforming system on the Kiruna Demonstrator Array. During the next week, we'll be featuring a few photographs from this busy time.

The Kiruna Demonstrator Array (Photo: D. McKay-Bukowski)

Current configuration, 6x4 dual-polarisation directional antennas. (Photo: D. McKay-Bukowski)

Lars-Göran digging out the duct entry points. Our previous experience from KAIRA
was absolutely invaluable for this part of the work. (Photo: D. McKay-Bukowski)

Laying the new (bright yellow!) ducting for the network cables. (Photo: D. McKay-Bukowski)

How do KAIRA and LOFAR work? — Part 10 : Multi-beam operations

As we saw in the last web log entry in this series, digitising the signals gives some distinct benefits when it comes to adding a delay. Apart from being able to do it very precisely, you can change it quickly. However, digital data has another massive advantage.

It can be easily copied.

So, if you make a copy of your digital data, you can process each copy however you want. If you want to filter it differently (for example, to ‘listen’ at a different frequency) you can. However, you can also make a copy and apply different sets of delays to it. This means that if one copy of the data has one set of delays and another copy has a different set for the different antennas involved, then these two copies will be sensitive to different directions.

In other words the telescope will be looking in two different directions at once!

This technique of introducing delays and combining the signals to change the direction in which the telescope is sensitive is known as ‘beamforming’. And making digital copies and delaying them differently to sensitive in several directions is called ‘multiple beamforming’ (or sometimes ‘forming multiple beams’).


In the above diagram (which you can click to see a larger version), each received signal is digitised and split into two copies. The copies are fed into the different digital delay units (D1, D2, D3 & D4). Let’s say that D1 and D3 have no delay at all. The added signal from these will be sensitive to something from directly overhead. In this case, the original radio signal was coming in obliquely; it adds out-of-phase and no output is seen. On the other hand, by adding a delay to D2 (but not D4), the added signal from these will be sensitive to the oblique radio wave. It adds in phase and a strong output is seen.

Because the data is digital, you can copy and delay it as much as you want and there is no degradation of the signal. If you want to be sensitive in more directions, simply add more copies, delay units and adders. The limits are not related to the signal, but rather to issues of equipment, memory, power, space (... and ultimately the project budget).

In practice, systems like KAIRA and LOFAR are also limited by the amount of digital signal processing equipment. At the time of writing this, LOFAR is routinely operating with up to 30 or so different beams. However, some experimental observations have been carried out with over 250 beams.