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.