Skip to content

Instantly share code, notes, and snippets.

@statnmap
Created March 10, 2019 13:59
Show Gist options
  • Save statnmap/4789cb9d0800ddb81db2627b13913d6b to your computer and use it in GitHub Desktop.
Save statnmap/4789cb9d0800ddb81db2627b13913d6b to your computer and use it in GitHub Desktop.
Transform hexbin object into sf polygon
library(hexbin)
library(sf)
# Create an hexbin
set.seed(101)
x <- rnorm(10000)
y <- rnorm(10000)
bin <- hexbin(x, y)
str(bin)
# st_as_sf(bin)
## or
plot(hexbin(x, y + x*(x+1)/4),
main = "(X, X(X+1)/4 + Y) where X,Y ~ rnorm(10000)")
# hexbin_to_sf ----
# hexbin_to_sf <- function(dat, crs = NA_crs_) {
library(grid)
library(sf)
library(hexbin)
library(dplyr)
crs <- NA_crs_
dat <- bin # hexbin object
xbins <- dat@xbins
shape <- dat@shape
tmp <- hcell2xy(dat, check.erosion = FALSE)
mincnt = 1; maxcnt = max(dat@count)
cnt <- dat@count
good <- mincnt <= cnt & cnt <= maxcnt
xnew <- tmp$x[good]
ynew <- tmp$y[good]
cnt <- cnt[good]
sx <- xbins/diff(dat@xbnds)
sy <- (xbins * shape)/diff(dat@ybnds)
inner <- 0.5
outer <- (2 * inner)/sqrt(3)
dx <- inner/sx
dy <- outer/(2 * sy)
rad <- sqrt(dx^2 + dy^2)
hexC <- hexcoords(dx, dy, sep = NULL)
# hexpolygon(hcell2xy(bin)$x, hcell2xy(bin)$y, hexC)
x <- hcell2xy(bin)$x
y <- hcell2xy(bin)$y
area <- pmin(area, maxarea)
radius <- sqrt(area)
n <- length(x)
n6 <- rep.int(6:6, n)
hUnit = "native"
fill = 1; border = 0
# grid.polygon(x = unit(rep.int(hexC$x, n) + rep.int(x, n6), hUnit),
# y = unit(rep.int(hexC$y, n) + rep.int(y, n6), hUnit),
# id.lengths = n6, gp = gpar(col = border, fill = fill))
# def.unit <- "native"
# grid.polygon(pltx, plty, default.units = def.unit, id = NULL,
# id.lengths = n6, gp = gpar(fill = pen, col = border))
pg <- polygonGrob(x = unit(rep.int(hexC$x, n) + rep.int(x, n6), hUnit),
y = unit(rep.int(hexC$y, n) + rep.int(y, n6), hUnit),
id = NULL,
id.lengths = n6,
default.units = def.unit,
# name = name,
gp = gpar(fill = fill, col = border)
# vp = vp
)
# str(pg)
n_coords <- rep(1:length(pg$id.lengths), times = pg$id.lengths)
l.pol <- lapply(1:length(pg$id.lengths), function(i) {
list(matrix(c(pg$x[n_coords == i], pg$x[n_coords == i][1],
pg$y[n_coords == i], pg$y[n_coords == i][1]),
ncol = 2
))
})
sf_pol <- sf::st_multipolygon(l.pol) %>%
sf::st_sfc(crs = crs) %>%
sf::st_cast("POLYGON") %>%
st_sf() %>%
mutate(cnt = cnt)
plot(sf_pol)
# }
@Maschette
Copy link

Thanks heaps for this, running the code I can't find maxarea? is this meant to be max(area)?

@Maschette
Copy link

Actually as far as I can tell area is not used? I go it to work without it :)

Thank you so much!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment