[GIS] Hex grid on EPSG:4326 globe (for use in Leaflet)

dggsdistancehexagonal grid

I'm looking for a way to evenly distribute a hexagon grid (x,y,z) on the EPSG:4326 globe. Knowing this is not mathematically possible I'm accepting the following limitations:

  • Poles don't need coverage – from lat -75 to lat 75 is sufficient.
  • International date line – 0/180 lng is OK to be messed up.
  • Distortion is OK as long as its not noticed by eye when viewing the grid at zoom ~100*100 hexes.

My hexagonal grids topology is defined according to this picture:
hex grid
(source: keekerdc.com)

Some further conditions

  • The hexagons are ~50 meter height.
  • I'm assuming hex(0,0,0) to be located at LatLng(0,0)

I need a function hexToLatLng(x,y,z) and a function LatLngToHex(lat,lng)
How can these lookup functions be created in JavaScript?

Below are some of my attempts.

By using the estimates

  • M meter in north direction = M *(1/111111) degrees lat
  • M meter in east direction = M *(1/(111111*cos(lat)) degrees lng

I implemented hexToLatLng and LatLngToHex. However increasing errors seems to sum up and give big shifts in longitude accuracy when moving away from equator. I guess its because the estimates assuming a sphere and EPSG:4326 is ellipse?

Because the hex system is planar I first tried EPSG:900913 and proj4js until I learned EPSG:900913 is not very suitable for distance calculations. I would also prefer sticking to EPSG:4326 as its easy to work with in my Leaflet map. But any thoughts on this track also most welcome.

I did find this thread but I'm hoping there are other solutions given I'm accepting the drawbacks at the poles
https://stackoverflow.com/questions/749264/covering-earth-with-hexagonal-map-tiles

Best Answer

It's true that you cannot evenly distribute hexagons across the 4326 globe, but you can get close using a discrete global grid. The introduction of 12 pentagons allows you to tile the rest of the globe quite nicely with hexagons. In particular, my R package dggridR allows you to easily construct such a grid:

library(dggridR)
#Generate a global grid whose cells are ~100,000 km^2
dggs         <- dgconstruct(area=100000, metric=FALSE, resround='nearest')
#Save the cells to a KML file for use in other software
gridfilename <- dgearthgrid(dggs,savegrid=TRUE)

The resulting KML file will have each cell labeled with a unique id.

We can also plot the grid to see that it is behaving as expected:

#Retrieve the same grid in a format useful for plotting
grid      <- dgearthgrid(dggs,frame=TRUE,wrapcells=TRUE,savegrid=FALSE)
#Get country outlines
countries <- map_data("world")

#Plot on a spherical map
ggplot() + 
    geom_polygon(data=countries, aes(x=long, y=lat, group=group), fill=NA, color="black")   +
    geom_polygon(data=grid,      aes(x=long, y=lat, group=group), fill="blue", alpha=0.4)   +
    geom_path   (data=grid,      aes(x=long, y=lat, group=group), alpha=0.4, color="white") +
    coord_map("ortho", orientation = c(-38.49831, -179.9223, 0))+
    xlab('')+ylab('')+
    theme(axis.ticks.x=element_blank())+
    theme(axis.ticks.y=element_blank())+
    theme(axis.text.x=element_blank())+
    theme(axis.text.y=element_blank())

A discrete global grid generated by dggridR

Note that one of the 12 pentagons is visible. These are, by default, orientated to be in out-of-the-way places. Future version of the package will allow you to orientate them.

The package is built on top of Kevin Sahr's DGGRID software. I haven't yet looked into how the cells are addressed for identifying neighbours (a future version of the package should address this), but it may be that the following publication would shed light on it:

Sahr, Kevin. “Location Coding on Icosahedral Aperture 3 Hexagon Discrete Global Grids.” Computers, Environment and Urban Systems 32.3 (2008): 174–187. CrossRef. Web.

Good luck!

Related Question