Display mean abundance (instead of count) in a hexagon

Question

In leaflet I'd like to show the mean abundance as a hexagon. I found "addHexbin" in leaflet.extras2 package but it appears to only count the number of records in each hexbin. What I really want is either sum or ideally the mean from another column "abun". I considered transforming the dataset so I have one row for each individual abundance record but some of the records have abundances of 20k+ which would result in a very large dataset and it still wouldn't calculate mean abundance in each hexagon. Calculating within the hexagon seems like a basic necessary function of hexbins so maybe I'm missing something obvious.

library(leaflet)
library(leaflet.extras2)
  
n <- 100000
df <- data.frame(lat = rnorm(n, 42.0285, .01),
                 lng = rnorm(n, -93.65, .01),
                 abun = rnorm(n, 20000, 1000))
  
leaflet(df)  %>%
  addTiles() %>%
  addHexbin(lng = ~lng, lat = ~lat,
            options = hexbinOptions(
              colorRange = c( "purple","blue","green", "yellow"),
              radiusRange = c(1, 20)
            )) 

I also tried using another unofficial R package Leaflethex which has the desired sum/mean function but doesn't integrate well into my Shiny App and caused some strange results (sidebar menu items disappearing, not working with leaflet proxy)

Answer 1

The solution I found was from a video by Jared Lander where he discusses using h3jsr h3 hexagon library as the base hexagon shape. The data is converted into a spatial object and matched with the H3 address. Any analysis can be done within this hexbin (mean/sum) and then can be plotted in leaflet.

library(dplyr)
library(h3jsr)
library(leafgl) #for faster plotting of hexbins
library(sf)

n <- 100000
df <- data.frame(lat = rnorm(n, 42.0285, 1),
                 lng = rnorm(n, -93.65, 1),
                 abun = runif(n, 0, 10000))


#dataframe to spatial object
df2 <- st_as_sf(x = df,                         
                coords = c("lng", "lat"),
                crs = "+datum=WGS84",
                na.fail=FALSE)


#Creates the hexbins for the coverage area and then calculate mean abun in each hexbin
# res=4 controls the size of the hexbin; smaller # = larger hexbins.  
hex5 <-df2 |>
  mutate(hex=point_to_cell(geometry, res=4)) |> #geometry column matches to H3 address
  st_drop_geometry()  |> #drop the geometry for faster processing
  group_by(hex) %>% #mean abundance by hex 
  summarise(mean_abun=mean(abun),
            .groups = 'drop')%>% 
  cell_to_polygon(simple=FALSE)#back to spatial 



#mapping the hexbins
leaflet("map") %>%
  addProviderTiles('CartoDB.Voyager')%>%addScaleBar(position='bottomleft') %>%
  addGlPolygons(data=hex5, fillColor="mean_abun",popup = ~mean_abun)