Sampling the surface of a sphere with a grid

2016-07-13

An answer to this question on Stack Overflow.

Question

I want to sample the surface of a sphere using a grid. The requirement is that each cell on the surface should have equal area. I saw that there are solutions like geodesic grids, but are there simpler solutions to that?

Answer

You don't specify if you have a language preference. I'm going to suggest a solution in R. If you need a different language, many of them have interfaces to R that you could use. This solution also produces an output file in a standard format that any language should be able to read.

You're right that geodesic grids are tricky. You're also right that they're the only way you're going to get cells of equal area.

Fortunately, I've built an R package called dggridR which makes it easy to work with them.

The following code accomplishes what you want by distributing hexagonal cells of equal area uniformly across the surface of the Earth. It also saves the vertices of these hexagonal cells to a file for use in other contexts.

#Include libraries
library(dggridR)
library(rgdal)
library(dplyr)
N <- 100    #How many cells to sample
#Distribute the points uniformly on a sphere using equations from
#http://mathworld.wolfram.com/SpherePointPicking.html
u     <- runif(N)
v     <- runif(N)
theta <- 2*pi*u      * 180/pi
phi   <- acos(2*v-1) * 180/pi
lon   <- theta-180
lat   <- phi-90
df    <- data.frame(lat=lat,lon=lon)
#Construct a global grid in which every hexagonal cell has an area of
#100,000 km^2. You could, of course, choose a much smaller value, but these
#will show up when I map them later.
#Note: Cells can only have certain areas, the `dgconstruct()` function below
#will tell you which area is closest to the one you want. You can also round
#up or down.
#Note: 12 cells are actually pentagons with an area 5/6 that of the hexagons
#But, with millions and millions of hexes, you are unlikely to choose one
#Future versions of the package will make it easier to reject the pentagons
dggs    <- dgconstruct(area=100000, metric=FALSE, resround='nearest')
#Get the corresponding grid cells for each randomly chosen lat-long
df$cell <- dgtransform(dggs,df$lat,df$lon)
#Get the hexes for each of these cells
grid    <- dgcellstogrid(dggs,df$cell,frame=FALSE)
#Save the hex boundaries to a handy KML file for use in your project
writeOGR(grid, "dgg_sample_cells.kml", "cells", "KML")

At this point you're done! You can parse weed out any highly unlikely pentagons by counting vertices in the KML, if need be.

For kicks, let's plot which polygons were chosen. Be wary of doing this if you have too many samples, though.

#Get the grid in a more convenient format
grid           <- dgcellstogrid(dggs,df$cell,frame=TRUE,wrapcells=TRUE)
#Get polygons for each country of the world
countries      <- map_data("world")
#Plot everything on a flat map
p<- 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="green", alpha=0.4)    +
    geom_path   (data=grid,      aes(x=long, y=lat, group=group), alpha=0.4, color="white")
p

[![Randomly selected cells on a Mercator projection of a discrete global grid][1]][1]

Note that the cells are distorted by the projection of this flat map. However, if we plot the cells on a sphere, we see that they are of equal area, as promised:

#Replot on a spherical projection
p+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())

[![Randomly selected cells on a spherically-projected discrete global grid][2]][2]

You could also get randomly selected cells like so:

#Choose a set of cells randomly distributed over the Earth
library(dggridR)
dggs    <- dgconstruct(spacing=1000, metric=FALSE, resround='down')
N       <- 100                                 #Number of cells
maxcell <- dgmaxcell(dggs)                     #Get maximum cell id
cells   <- sample(1:maxcell, N, replace=FALSE) #Choose random cells
grid    <- dgcellstogrid(dggs,cells,frame=TRUE,wrapcells=TRUE) #Get grid

[1]: https://i.sstatic.net/tVVxf.png [2]: https://i.sstatic.net/483S1.png