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.

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).

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).

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.

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.

As the world turns: projecting a 2D image onto a 3D sphere

How difficult could it be two project a two dimensional image onto a three dimensional sphere in R? Well, it's not really that difficult, but it took me quite a while to figure out how to achieve this. Now why would we want to do this in the first place? Well, you could use such a technique to create a three dimensional model of Earth on which you can plot all sorts of interesting global phenomena. Well, I couldn't think of many other useful applications besides modelling fancy Christmas ornaments. For now we'll stick to the model of Earth.

This blog was rewritten several times, before publication, as I tried several strategies with varying degrees of success. I started with the persp function from the graphics package. Which just didn't work. This function can plot a single hemisphere, with limited control over the colours used for its surface.

Then I ended up using the spheresurf3D function from the plot3D package. Although I finally got it to work, it required a horrific workaround. The plot3D package uses a matrix of colour values that will be used in combination with a colour ramp to visualise on the surface. My workaround was to get all unique values from the target image and create a matrix with indices for those colours. I will not present this strategy in detail here due to this monstrosity and the other graphical limitations of this package. Don't get me wrong, I think plot3D is a great package, with many useful applications, it's just not the ultimate tool for what I want here.

The easiest way to create a three dimensional model of the globe was with the rgl package. But first, some code to prepare ourselves. The code below downloads a topographic map from NASA's website, which we will use to project onto a sphere. We convert the jpeg image into a png file, as the rgl package will only accept that file type. Note that free third party software (ImageMagick) is used for this step. We also specify a function to convert polar coordinates (longitude, latitude and altitude) into Cartesian coordinates. We also specify a matrix with Cartesian coordinates (‘world_coords’) which we will use to georeference the topographic image.

Next, we set up an rgl device and our model. We create an environment with a black background (space), and set up a light to represent the sun. With the persp3d function we create a model of Earth. We use the ‘texture’ argument to project the topographic image onto the sphere. We add a slightly larger transparent light blue sphere to represent the Earth's atmosphere. Finally, we add a small red sphere hovering above my home country. That's it! We have a simple model of our planet. The nice thing of rgl devices is that they are interactive. You can rotate your model and zoom in and out. With the writeWebGL you can even save your interactive model as a web page. For me, this web page worked offline, but due to restrictions of blogspot, I could get it to work on my blog.

In order to produce something that I can show here, I created an animation. Therefore I first changed the orientation of the model, and then animated a spinning world with the movie3d function.

Once you know how to do it, it's not that difficult to create an animation of a rotating planet in R. Of course there is always stuff that I could improve upon. For example, the Earth isn't a perfect sphere, but I had assumed so as a simplification. If you have any refinements for the approach shown here, please let me know in the comments.

Filtering images with Fast Fourier Transform

Fourier transformation can be a useful instrument in the analyses, processing and compression of data such as images or sound. There is a tremendous amount of information on two-dimensional discrete fast Fourier transformation available online. However, clearcut examples in R are scarce. It took my quite a while to piece together all information and figure out how to get it done in R. With this post I neither intend to explain the principles of fast Fourier transformation nor the mathematics behind it. However, I couldn't avoid talking about some technical stuff, without explaining it too much. Therefore, it is advisable to look around for some background info on the topic, before delving into this post. For starters you could have a look at Steven Lehar's page. Hopefully you can use the example presented here to your advantage.

For those readers I didn't scare off yet: let's get started. With a Fourier transform, information from discrete space (e.g. the brightnes of your image) is transformed to a frequency domain. This sounds very technical, which it is. If you took the time Googling a bit on the topic, you will understand that sinusoidal patterns in an image will show up as peaks in the frequency domain spectrum generated with a Fourier transform. One application of the Fourier transform is to detect and remove regular distortion patterns in images. In the example provided here, I will guide you through the following steps:

• Read an image and store as a matrix in the working memory
• Transform the image to the frequency domain
• Do some filtering magic
• Back transform the image to normal space and save the filtered image

Let's start with the image below. I've taken this picture years ago in the wonderful Dutch amusement park De Efteling. I've added some distortion, the kind that you will typically get when scanning a picture from a newspaper.

First things first. We first need to load the jpeg in memory such that we can process it as a matrix. Let's do this with the jpeg package. We also need to rotate the matrix as the jpeg uses a different origin for plotting:

Transforming the image to the frequency domain in R is actually ridiculously simple. Just one call to the fft function is enough. The multiplication with the ‘checkerboard’ is to ensure that the frequency spectrum is centred in the middle. The call to the fft function results in a matrix with the same dimensions as the original image. However, the matrix contains complex numbers. A complex number is composed of a real part and an imaginary part, which can be represented as a vector on two axes. The length of this vector (the modulus, calculated with the Mod function) gives the magnitude of frequencies in the spectrum. The angle of the vector (the argument, calculated with the Arg function) gives the phase shift of the sinusoidal waves. The latter is not very interesting for filtering of the image, but crucial in the reconstruction of the image after filtering.

When the magnitude of frequencies is plotted, the spectrum as shown below is obtained. Peaks in the spectrum are shown as white areas. The white blur in the centre of the spectrum contains much of the information of the original image. Other peaks circled in red are the result of the regular distortions in the image. Using the locator function we can now manually select the circled peaks with mouse clicks, so we have their location in the spectrum.

Next step is to filter out the selected (red circled) peaks from the spectrum. Basically, what we need to do is set the magnitude to zero at the location of the peaks. I've written a fancy function that does that, only slightly more complicated. The function multiplies the selected location with a bell-shaped curve (with the bell upside down). If you don't like this fancy approach, you set the ‘gaussian’ argument in the function to FALSE. When you do that, all values with a radius of r around the selected location in the spectrum are set to zero. We're all done filtering the unwanted peaks from the spectrum.

Note that we modified the magnitudes in the frequency spectrum. We left the phase shifts unchanged. We need to reconstruct the complex values from the modified magnitudes and the phase. This is not complex at all using the complex function (pardon the awful wordplay). Next we can apply the inverse Fourier transform, by simply setting the parameter ‘inv’ to TRUE in the fft function.

Just like that we have our filtered image which we can save with the writeJPEG function:

Before filter

After filter

I know, the image is a bit blurry, noisy and a bit pale, but at least we got rid of the nasty distortion (which I put there myself in the first place). As you can imagine, there are plenty of other applications possible with Fourier transform. Maybe I'll save those for a later blog...