[R-sig-Geo] Build hexagonal lattice coverage

[Tim Keitt]

Thanks Roger!

I went ahead and wrote a couple of functions:

genHexGrid <- function(dx, ll = c(0, 0), ur = c(1, 1))
{
        dy <- sqrt(3) * dx / 2

        x <- seq(ll[1], ur[1] - dx / 2, dx)
        y <- seq(ll[2], ur[2], dy)

        y <- rep(y, each = length(x))
        x <- rep(c(x, x + dx / 2), length.out = length(y))

        x <- x + (ur[1] - max(x)) / 2
        y <- y + (ur[2] - max(y)) / 2

        list(x = x, y = y)
}

genPolyList <- function(hexGrid)
{
        dx <- hexGrid$x[2] - hexGrid$x[1]
        dy <- dx / sqrt(3)

        x.offset <- c(-dx / 2, 0, dx / 2, dx / 2, 0, -dx / 2)
        y.offset <- c(dy / 2, dy, dy / 2, -dy / 2, -dy, -dy / 2)

        f <- function(i) list(x = hexGrid$x[i] + x.offset,
                              y = hexGrid$y[i] + y.offset)

        lapply(1:length(hexGrid$x), f)
}

xy <- genHexGrid(0.05)
plot(xy, asp = 1, axes = F, pch = 19, xlab = NA, ylab = NA, cex = 0.5)
lapply(genPolyList(xy), polygon, xpd = NA, ypd = NA)

Now to figure out the sp stuff.

Cheers,
Tim

On 3/3/07, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>
>
>


-- 
Timothy H. Keitt, University of Texas at Austin
Contact info and schedule at http://www.keittlab.org/tkeitt/
Reprints at http://www.keittlab.org/tkeitt/papers/
ODF attachment? See http://www.openoffice.org/