Grid2Polygons

June 25, 2012
By

(This article was first published on 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)

center

Figure 1: Simple spatial grid data frame.

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)

center

Figure 2: Simple gridded data represented with spatial polygons.

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)

center

Figure 3: Distance from river represented with spatial polygons.

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

center

Figure 4: Topographic information represented with a raster image.

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, at)
cols <- terrain.colors(length(at) - 1)[col.idxs]
plot(dem.plys, border = "transparent", col = cols)

center

Figure 5: Topographic information represented with spatial polygons.

Finally, transform the projection of the spatial polygons object (1.5 seconds) and replot (fig. 6).

dem.plys.trans <- rgdal::spTransform(dem.plys, CRS = CRS("+proj=eqc"))
plot(dem.plys.trans, border = "transparent", col = cols)
par(op)

center

Figure 6: Topographic information represented with spatial polygons and new projection.

To leave a comment for the author, please follow the link and comment on their blog: jfisher-usgs.

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



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Sponsors

Mango solutions



plotly webpage

dominolab webpage



Zero Inflated Models and Generalized Linear Mixed Models with R

Quantide: statistical consulting and training

datasociety

http://www.eoda.de





ODSC

ODSC

CRC R books series





Six Sigma Online Training









Contact us if you wish to help support R-bloggers, and place your banner here.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)