Both SpatialPolygons and SpatialPolygons consist of list of Polygons. Each of these Polygons in turn consist of vertices, i.e., the points where the edges of the polygon meet. When you transform a SpatialPolygons(DataFrame) to a different projection (e.g., from long-lat coordinates to UTM) with the spTransform-function, actually only the vertices are transformed (and not the lines connecting them). This means that also in the new projection straight lines are drawn between the vertices, whereas in reality they should be curved (due to the quasi-spherical shape of our world). For most shapes this isn't much of an issue. But when your shape covers a large area, you could run into serious problems.
This problem can be solved by adding more vertices to the polygons. In other words, add new vertices by interpolating between the existing vertices. If you do this before re-projecting, the new projection will be a lot more accurate. I tried to find an existing package in R that would do this, but I could not find any (please leave a comment if you did find an alternative). Therefore, I've written a small function to take care of this myself, which I will demonstrate in this post.
The source code of this function and the accompanying example is given below. This code is also used to produced the graphical illustration below. Note how the red line crosses the islands Northeast of Scotland in the original projection. Compare this with the new projection of the red lines with and without segmentation. Also note that the red lines in the new projection are straight when we don't apply segmentation, and they are curved (as they are supposed to be) when we do apply segmentation.
The segmentation is achieved relatively simple by looping all the lists of Polygons with the lapply function and then looping each line segment with the apply function (see R script below). The new vertices are derived via linear interpolation using the approx function. You could go for something more fancy like spline interpolation, but linear will do just fine in this example. As I've implemented all this in the function named ‘SPSegmentation’, you can simply call it on any SpatialPolygons or SpatialPolygonsDataFrame object you would like. With some modifications, you could also use this function of SpatialLines and SpatialLinesDataFrame objects.
In case you want to achieve the opposite (i.e., reduce the number of vertices in a SpatialPolygons(DataFrame)), you can use the gSimplify function from the rgeos package. This uses a Ramer–Douglas–Peucker algorithm to simplify the shape, which could save you some computing time, but your shape will also become less detailed as a result.
| ## package for handling spatial data: | |
| library(sp) | |
| ## package that contains simple world map data | |
| library(maptools) | |
| ## package for modifying spatial data: | |
| library(rgeos) | |
| ## package used for setting up clipping areas | |
| library(raster) | |
| ## function that segments the polygons of SpatialPolygons(DataFrame) object | |
| ## 'object' by a factor of 'segments', where 'segments' should be an | |
| ## integer value >= 1 | |
| SPSegmentation <- function(object, segments) { | |
| ## apply to each element of the SpatialPolygons object: | |
| object@polygons <- lapply(object@polygons, function(pol) { | |
| ## apply to each polygon of each element of each SpatialPolygons object: | |
| pol@Polygons <- lapply(pol@Polygons, function(p) { | |
| ## apply linear interpolation to the coordinates: | |
| p@coords <- apply(p@coords, 2, function(x) { | |
| approx(x, n = (nrow(p@coords) - 1)*segments + 1)$y | |
| }) | |
| return(p) | |
| }) | |
| return(pol) | |
| }) | |
| return(object) | |
| } | |
| data(wrld_simpl) | |
| ## first clip the world data to an area useful for our example... | |
| ## start be defining the clipping area: | |
| clip <- as(extent(-10, 15, 50, 60), "SpatialPolygons") | |
| ## add the proper projection to the data: | |
| proj4string(clip) <- proj4string(wrld_simpl) | |
| ## and clip it: | |
| wrld_clip <- gIntersection(wrld_simpl, clip, byid=TRUE) | |
| ## setup a simple rectangle: | |
| simple_rect <- as(extent(-8, 13, 51, 59), "SpatialPolygons") | |
| ## add the proper projection to the data: | |
| proj4string(simple_rect) <- proj4string(wrld_simpl) | |
| ## transform the clipped world map to UTM: | |
| wrld_clip_trans <- spTransform(wrld_clip, CRS("+proj=utm +zone=31 +datum=WGS84")) | |
| ## transform the simple rectangle without segmentation: | |
| simple_rect_trans <- spTransform(simple_rect, CRS("+proj=utm +zone=31 +datum=WGS84")) | |
| ## now first segment the simple rectangle (increasing the | |
| ## number of vertices by a factor of 10) | |
| segmented_rect <- SPSegmentation(simple_rect, 10) | |
| ## transform the segmented rectangle to UTM: | |
| segmented_rect_trans <- spTransform(segmented_rect, CRS("+proj=utm +zone=31 +datum=WGS84")) | |
| ## plot everything: | |
| png("segmentation-demo.png", 1200, 400, res = 300, type = "cairo") | |
| par(mfcol=c(1,3), mar = c(0, 0, 2, 0)) | |
| plot(wrld_clip, col = "lightgreen", main = "original projection") | |
| plot(simple_rect, add = T, border = "red") | |
| plot(wrld_clip_trans, col = "lightgreen", main = "new projection\nno segmentation") | |
| plot(simple_rect_trans, add = T, border = "red") | |
| plot(wrld_clip_trans, col = "lightgreen", main = "new projection\nsegmentation") | |
| plot(segmented_rect_trans, add = T, border = "red") | |
| dev.off() |
