Transform hexbin object into sf polygon

hexbin_to_sf

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

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

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

Thank you so much!