Filling SpatialPolygons with an image

Plotting country outlines is really easy in R, but making those plots a bit more fancy can be frustrating. I thought it would be nice to fill the country outlines with an image rather than with a solid colour. How hard could it be? After some googling it appears that there is very little documentation on this topic. The only document that I found was by Paul Murrel (2011). He plots a black filled country shape first and then captures the rasterized shape using the grid.cap function from the grid package. This captured raster is then used as a mask on the image that needs to be plotted in the shape of a country.

I'm having difficulties with this approach for two reasons: one, I can't seem to get grid.cap working as it should; two, I don't like having to plot the shape of a country first in order to create a mask. The first matter is probably my own wrongdoing, but all the grid.cap function returns me is a matrix of white pixels. I grew tired of figuring out what I was doing wrong. For the second aspect, I think a more direct approach should be possible by using the ‘over’ function from the sp package.

So these are the steps that I took to tackle the problem:

1. Download the country shape
2. Download a suitable png image
3. Georeference the image, such that it matches the location of the country
4. Use the ‘over’ function to determine which pixels are inside the country shape
5. Plot the pixels of the image that are inside the country and plot the outline

In this post I will show you how to create the image shown below, where each country is filled with an image of their respective flag. Below I will explain in a bit more detail how this image was created, following the steps listed above. I've also provided the script used to create the image below, where each of the steps are also clearly marked.

The country shapes for step 1 are downloaded from GADM using the getData function from the raster package. Of course you can use any SpatialPolygons object from any other source.

The country flag png images are downloaded from wikipedia (step 2). So this example shows how to download png image and turn it into a usable format. If you like, you can use any image format like jpg or tif, but that will require some modification of the code presented here.

To understand the next step (3), I'm using georeferencing and I may need to explain what that means. Basically, this is telling where and how the image should be positioned in the coordinate reference system (CRS) of the SpatialPolygons objects (the country outlines). In this step I'm instructing the system to stretch the image to the bounding box of the respective country. This also means that the aspect ratio of the original image is probably messed up. If you would like to keep the original aspect ratio, you need to modify this step. The image was downloaded as an array with the red, green and blue component in separate dimensions. With the brick function from the raster package this array is turned into a correctly georeferenced brick object.

Now that we have properly positioned the image, we need to determine which pixels are located inside the country outlines (step 4). For this purpose, we use the ‘over’ function from the sp package. This function will not accept a raster brick object as input. The raster brick object is therefore cast into a SpatialGrid format. The function will return a dataframe with the country shape element that matches with a specific grid cell. From this information it can be derived whether the pixel is situated inside or outside the country outline. The pixel values for those situated outside the country outlines are set to NA. Note that this can be a time-consuming step. Especially when the resolution of your image is high, and/or the country outline is complex (i.e., contains details, a lot of islands and/or holes), so either be patient. Alternatively, you can speed up the process by either lowering the resolution of your image (hint: focal) or simplifying the country outline (hint: gSimplify).

All there is left to do is to plot the image and the outline (step 5). I use the plotRGB function from the raster package to plot the image. Don't forget to set the ‘bgalpha’ argument to 0, to ensure that NA values are plotted as transparent pixels. Otherwise, it will plot white pixels over anything that has previously been plotted.

I've wrapped four of the five steps in a function, such that I could easily repeat these steps for different countries. Hopefully this post will help you to create your own cool graphics in R. Good luck, and please let me know in the comments if any of the steps are not clear...

2D particle tracking in R, part 2: kernel density

In my previous post, I showed how particle movement can be simulated in the marine environment. Now that we have tracked the particle movement, how can we use this information to translate this to something more meaningful? This can be done by determining the particle density on the map.

The most simple way would be to lay a grid over the particles and simply count the number of particles in each grid cells. However, I would like to go for something a little bit complicated that gives some smoother results. I will use kernel density estimates. This means that I will consider the particles as a density function, were the position of the particle gives the highest density for that particle. You can get a density map of all particles by simply adding up the probabilities of the density functions of each individual particle. Fortunately, the KernSmooth package will take care of all this.

Remember that in the original simulation, I released 5,000 particles within a specific time-frame. With the kernel density estimated, I can estimate the number of particles per square kilometer at any time and location. This will become even more interesting if these particles represent something in the real world, let's say oil. If we know how much cubic meter of oil is represented by each particle, the estimates can be transformed in to the thickness of the oil slick. Keep in mind that this is just an initial test simulation that does not include all processes in order to properly simulate oil spills (e.g., oil evaporation, degradation and dispersion is not simulated).

Here's a video created with ffmpeg of the particle density movement over time.

And here's the fully commented script used to create the frames for the video.

2D particle tracking in R

How hard is it nowadays to simulate particle movement in the sea and track them in R? Well I discovered that it is not all that difficult. Within a couple of hours I had a simple simulation model running in the North Sea. All you need are time dependent water currents and wind data, both of which can be obtained from marine.copernicus.eu.

The main steps that you will have to follow:

• load the water currents and wind data
• specify particle release details
• specify simulation details
• loop all simulation time steps and for each step do the following with the particles that are released:
• for each particle determine the local water current and wind vector at the current time index
• transform the particle coordinates to a projection in meters
• move the particles with the water current (m/s) multiplied by the time step duration (s)
• move the particles with the wind speed (m/s) multiplied by the time step duration (s) and a drift factor (between 0 and 1)
• move the particles in random directions (to simulate diffusion)
• transform the particle coordinates back to the original projection
• save particle information every n^th simulation step, where n is any positive integer value

Below you will find the script that does all this (each of the steps above are marked in the script). In this simulation I simulate the release of 5,000 particles with a constant release rate during 5 days and followed them for a couple of days. I found the trickiest bit not the simulation itself, but reading the the NetCDF files such that it could be used in a convenient manner. I create a video of this simulation with ffmpeg:

So what's the practical use of all this? Well obviously you can simulate the trajectory of any floating materials, such as oil. However, you have to keep in mind that the simulation is only as good as the underpinning data and parameterisation. For instance, the water current data used here is not that accurate near estuaries. Of course you can always do multiple runs of your simulation with different parameter settings (or make it even a Monte Carlo simulation). Also, here I only simulate in 2D. But there is nothing that would stop you from modifying it to include the third dimension. What also should be added is some detection that particles have beached. But that should not be very difficult when you have coastline data.

Why in R? Well, commercial particle tracking software is generally expensive. Whereas free software are often rigid and don't allow the user to include custom modifications. In R the simulation is free, and you can take the exact implementation into our own hands... Hopefully, this post will help you get started with your own simulations. Good luck!

Here's the script:

Edit: if your interested, you might want to check out the follow-up post on particle density estimates of this simulation.

ContourLines to SpatialPolygons

Having a very large gridded data set (e.g., raster or SpatialGrid) can be very useful for analyses. However, for visualisation it can be tedious to have a very large spatial grid, consuming a lot of memory and processing time. Converting such gridded data to contour lines could save a lot of time, especially when you want to reuse the data. Also, it is much easier to reproject polygons to a different coordinate system than gridded data.

Although there are several proper function to create contour lines from matrix data, creating outlines in the form of a SpatialPolygonsDataFrame for the purpose of reusing it is a pain. This post shows the problems I ran into and how I solved them. Improvements and suggestions are welcome in the comments below.

I've tried several strategies but I only present the strategy I got to work here. The example uses the volcano matrix. Although converting gridded data to SpatialPolygons becomes particularly useful for large grids, I use the small volcano data set as prove of principle.

My first attempt uses the functions rasterToContour, SpatialLines2PolySet, and PolySet2SpatialPolygons functions for conversion steps.

Of course, it couldn't have been that easy... The above script produces this:

Apparently, the PolySet2SpatialPolygons does not close the polygon rings correctly at the edge of the volcano matrix as I would have hoped. Now, I have to confess, I work a lot with ‘quick and dirty’ solutions, as is the case here. In order to fix the problem with the polygons that intersect with the edge of the volcano matrix, I simply modify the original matrix. I simply add an extra row and column to each side of the matrix, with a value smaller than the minimum value in the original matrix (in this case I use the value 60). This creates an artificial ridge in the dataset that forces the routines to correctly close the polygons:

The resulting plot shows that we are almost there:

There is still a problem that is hardly visible. Take a closer look at the crater of the volcano. The crater should be deeper, which isn't shown by the colour gradient. Some of the polygons should contain holes, which aren't automatically detected and converted by the PolySet2SpatialPolygons function. Maybe there are better solutions, but I chose to manually hack the SpatialPolygons object. I did this by looping through all polygons and when a polygon is inside another polygon of the same layer, it should be a hole. In that case, the ‘hole’-flag should be set to TRUE and the coordinate order should be reversed (hole ring direction should be opposite to that of the parent polygon). Not elegant, but hey, it works:

Almost there again:

The only thing left to do is to remove the artificial edge that we created earlier. This can easily be done with the function gIntersection:

There you are. We have successfully converted the grid into a much lighter (at least compared to large grids) SpatialPolygons object which we can save, load and reuse in plotting routines. I realise that the SpatialPolygons object in this example isn't georeferenced. But the principles shown here can also be applied to georeferenced rasters and SpatialGrid objects.