Hexagon binning using Google Maps

An answer to this question on the GIS Stack Exchange.

Question

I have a bunch of geolocated points that I am looking to cluster on a google map using hexagons. These clusters will be stored in a db for quick retrieval and will be reduced to a single point and radius. I am using the http://www.geocodezip.com/v3_MW_example_eshapes.html library to draw the hexagons but because of the mercator projection the hexs overlap at higher latitudes.

Here is a jsfiddle attempt with the overlap http://jsfiddle.net/kworthin/Lkuqykog/2

overlapping hexes

Answer

If you're pregenerating the hex information, as your question suggests, you could use my R package dggridR to generate the hexagons' boundaries.

dggridR uses some rather complicated maths to tile the entire Earth with non-overlapping, equally-sized hexagons.

Here's an example of how you might use dggridR to create the grid you need:

#Include libraries
library(dggridR)
library(dplyr)
library(OpenStreetMap)
library(rgdal)
#Construct a global grid with cells approximately 10 miles across
dggs          <- dgconstruct(spacing=10, metric=FALSE, resround='down')
#Load included test data set
nylat = 40.7127 
nylon = -74.0059
#Generate a bunch of random points
lats <- runif(100,nylat-1,nylat+1)
lons <- runif(100,nylon-1,nylon+1)
df   <- data.frame(lat=lats, lon=lons)
#Get the corresponding hexagonal grid cells for each point
df$cell <- dgtransform(dggs,df$lat,df$lon)
#Get the number of points in each cell
pcounts <- df %>% group_by(cell) %>% summarise(count=n())
#Get the grid cell boundaries for the the points in a form suitable for plotting
grid    <- dgcellstogrid(dggs,df$cell,frame=TRUE)
#Update the grid cells' properties to include the number of points
grid    <- merge(grid,pcounts,by.x="Name",by.y="cell")
#Get a street map of the area
map   <- openmap(c(nylat+1,nylon-1),c(nylat-1,nylon+1),type='osm')
mapLL <- openproj(map)
plot(mapLL)
#Plot everything on a flat map
p<-autoplot(mapLL)
p<- p+    geom_polygon(data=grid,      aes(x=long, y=lat, group=group, fill=count), alpha=0.4)    +
    geom_path   (data=grid,      aes(x=long, y=lat, group=group), alpha=0.4, color="white") +
    scale_fill_gradient(low="blue", high="red")
p

Discrete global grid on New York City

Okay, it works! Now let's save the grid, with the data, to somewhere it can be read by other programs:

#Get the grid cell boundaries in a form suitable for printing to a KML file
grid <- dgcellstogrid(dggs,df$cell,frame=FALSE)
#Update the grid cells' properties to include the number of earthquakes
#in each cell
grid@data$count <- merge(grid@data, pcounts, by.x="Name", by.y="cell", all.x=TRUE)$count
#Write out the grid
writeOGR(grid, "ny_grid.kml", "cells", "KML")

1: https://i.sstatic.net/nh119.png