[R-sig-Geo] Build hexagonal lattice coverage

[Roger Bivand]

On Fri, 2 Mar 2007, Tim Keitt wrote:

> Hi All,
> 
> I'd like to make a hexagonal lattice (stored in a shapefile, ogr
> coverage or postgis geometry). Is the general approach to make a
> series of independent hexagons and place them in the plane such that
> adjacent edges sit directly on top of each other? Or can one build a
> mesh out of connected line segments? If it is the latter, how does one
> associate a data value with each hexagon in the mesh? I'd like to
> achieve something like the output of the hexbin package, but display
> the result in ArcGIS or QGIS.

This is a bit "cheap", but might help (for points on equilateral 
triangles of side equal to 1):

library(deldir)
x1 <- seq(1,10,1)
x2 <- x1 + 0.5
step <- cos(30*(pi/180))
y1 <- seq(1, 10, step*2)
y2 <- seq(1+step, 10, step*2)
xy1 <- expand.grid(x1, y1)
xy2 <- expand.grid(x2, y2)
crds <- as.matrix(rbind(xy1, xy2))
z <- deldir(crds[,1], crds[,2])
w <- tile.list(z)
polys <- vector(mode="list", length=length(w))
library(sp)
for (i in seq(along=polys)) {
  pcrds <- cbind(w[[i]]$x, w[[i]]$y)
  pcrds <- rbind(pcrds, pcrds[1,])
  polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
}
SP <- SpatialPolygons(polys)
SPDF <- SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1], 
  y=crds[,2], row.names=sapply(slot(SP, "polygons"),
  function(x) slot(x, "ID"))))
Areas <- sapply(slot(SP, "polygons"), function(x) slot(x, "area"))
print(fivenum(Areas), digits=8)
print(sort(Areas), digits=8)
SPDF1 <- SPDF[Areas < 0.866027,]
plot(SPDF1, axes=TRUE, xlim=c(1,10), ylim=c(1,10))
points(crds)
library(maptools)
writePolyShape(SPDF1, "hex")
# or writeOGR() to other formats.

The areas are not exactly equal, varying in the sixth digit after the 
decimal point by about 1. But having deldir build the tesselation saves 
the mess of doing it directly, and the hex grid of points is easier to 
build. The default bounding box of tile.list() seems to be OK.

The external points do not get hexagons, so the initial setup will need to 
take this into account, by adding external buffer points.

Using tile.list() from deldir seems easier than alternatives from tripack, 
the original tesselation question about using it came from Simon Goring, 
who had done most of the work apart from building the Polygons objects.

Hope this helps,

Roger

> 
> THK
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no