**jfisher-usgs**, and kindly contributed to R-bloggers)

I’d like to introduce you to the *Grid2Polygons* function; an R function for converting **sp** spatial objects from class *SpatialGridDataFrame* to *SpatialPolygonsDataFrame*. The significance of this conversion is that spatial polygons can be transformed to a different projection or datum with the *spTransform* function in package **rgdal**. Postscript files created with spatial polygons are reduced in size and result in a much “cleaner” version of your image. Disadvantages of the conversion include long computational times and irreversible leveling, partitioning the range of *z* values. A general explanation of the algorithm is provided here; inspiration provided here.

To access the function install the Grid2Polygons package (source):

```
install.packages("Grid2Polygons")
library(Grid2Polygons)
```

See help documentation for argument descriptions:

```
?Grid2Polygons
```

The following examples highlight the functions usefulness:

## Example 1

Construct a simple spatial grid data frame.

```
m <- 5
n <- 6
z <- c(1.1, 1.5, 4.2, 4.1, 4.3, 4.7,
1.2, 1.4, 4.8, 4.8, NA, 4.1,
1.7, 4.2, 1.4, 4.8, 4.0, 4.4,
1.1, 1.3, 1.2, 4.8, 1.6, NA,
3.3, 2.9, NA, 4.1, 1.0, 4.0)
x <- rep(0:6, m + 1)
y <- rep(0:5, each = n + 1)
xc <- c(rep(seq(0.5, 5.5, by = 1), m))
yc <- rep(rev(seq(0.5, 4.5, by = 1)), each = n)
grd <- data.frame(z = z, xc = xc, yc = yc)
coordinates(grd) <- ~ xc + yc
gridded(grd) <- TRUE
grd <- as(grd, "SpatialGridDataFrame")
```

Plot the grid using a gray scale to indicate values of *z* (**fig. 1**).

```
image(grd, col = gray.colors(30), axes = TRUE)
grid(col = "black", lty = 1)
points(x = x, y = y, pch = 16)
text(cbind(xc, yc), labels = z)
text(cbind(x = x + 0.1, y = rev(y + 0.1)), labels = 1:42, cex = 0.6)
```

Convert the grid to spatial polygons and overlay in plot (**fig. 2**). Leveling is specified with cut locations at 1, 2, 3, 4, and 5, and *z*-values set equal to the midpoint between breakpoints. A “winding rule” is used to determine if a polygon ring is filled (island) or is a hole in another polygon.

```
plys <- Grid2Polygons(grd, level = TRUE, at = 1:5)
cols <- rainbow(4, alpha = 0.3)
plot(plys, col = cols, add = TRUE)
x <- rep(0:6, m + 1)
y <- rep(0:5, each = n + 1)
legend("top", legend = plys[[1]], fill = cols, bty = "n", xpd = TRUE, inset = c(0, -0.1), ncol = 4)
```

## Example 2

Apply the conversion function to the *meuse* data set, included in the **sp** package. The effect of leveling is shown in **figure 3**.

```
data(meuse.grid)
coordinates(meuse.grid) <- ~ x + y
gridded(meuse.grid) <- TRUE
meuse.grid <- as(meuse.grid, "SpatialGridDataFrame")
meuse.plys <- Grid2Polygons(meuse.grid, "dist", level = FALSE)
op <- par(mfrow = c(1, 2), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0))
z <- meuse.plys[[1]]
col.idxs <- findInterval(z, sort(unique(na.omit(z))))
cols <- heat.colors(max(col.idxs))[col.idxs]
plot(meuse.plys, col = cols)
title("level = FALSE", line = -7)
meuse.plys.lev <- Grid2Polygons(meuse.grid, "dist", level = TRUE)
z <- meuse.plys.lev[[1]]
col.idxs <- findInterval(z, sort(unique(na.omit(z))))
cols <- heat.colors(max(col.idxs))[col.idxs]
plot(meuse.plys.lev, col = cols)
title("level = TRUE", line = -7)
par(op)
```

## Example 3

A real-world example using topographic information on the eastern Snake River Plain; raster plot shown in **figure 4**.

```
library(rgdal)
data(DEM)
at <- seq(500, 4000, by = 250)
op <- par(oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0))
image(DEM, breaks = at, col = terrain.colors(length(at) - 1))
```

Convert the grid to spatial polygons (**fig. 5**). The conversion took 22.7 seconds on my machine. The size of the postscript file created from the polygon image (1.89 MB) is 82.6 percent smaller than the size of the raster image file (10.86 MB).

```
dem.plys <- Grid2Polygons(DEM, level = TRUE, at = at)
z <- dem.plys[[1]]
col.idxs <- findInterval(z, sort(unique(na.omit(z))))
cols <- terrain.colors(max(col.idxs))[col.idxs]
plot(dem.plys, border = "transparent", col = cols)
```

Finally, remove the Albers projection from the spatial polygons object (1.5 seconds) and replot (**fig. 6**).

```
dem.plys.ll <- rgdal::spTransform(dem.plys, CRS = CRS("+proj=longlat +datum=WGS84"))
plot(dem.plys.ll, border = "transparent", col = cols)
par(op)
```

**leave a comment**for the author, please follow the link and comment on his blog:

**jfisher-usgs**.

R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...