geohashTools
In geohashTools: Tools for Working with Geohashes

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = '#>',
  warning = FALSE,
  message = FALSE
)

What is a Geohash?

Developed by Gustavo Niemeyer, the geohash is system of nestable, compact global coordinates based on Z-order curves. The system consists of carving the earth into equally-sized rectangles (when projected into latitude/longitude space) and nesting this process recursively. Geohashes are a grid-like hierarchical spatial indexing system. The precision of a geohash is dictated by the length of the character string that encodes it.

The geohash system is useful for indexing or aggregating point data with latitude and longitude coordinates. Each geohash of a given length uniquely identifies a section of the globe.

Precision of geohashes

data.frame(
  `Geohash Length` = 1:8,
  `KM error`= c(2500.0, 630.0, 78.0, 20.0, 2.4, 0.61, 0.076, 0.019)
)

Encoding geohashes

Encoding is the process of turning latitude/longitude coordinates into geohash strings. For example, Parque Nacional Tayrona in Colombia is located at roughly 11.3113917 degrees of latitude, -74.0779006 degrees of longitude. This can be expressed more compactly as:

library(geohashTools)
gh_encode(11.3113917, -74.0779006)

These 6 characters identify this point on the globe to within 1.2 kilometers (east-west) and .6 kilometers (north-south).

The park is quite large, and this is too precise to cover the park; we can "zoom out" by reducing the precision (which is the number of characters in the output, 6 by default):

gh_encode(11.3113917, -74.0779006, precision = 5L)

Example: Encoding many points

We can use this as a simple, regular level of spatial aggregation for spatial points data. Here with randomly-selected coordinates:

coords = data.frame(
  x=rnorm(20L),
  y=rnorm(20L)
)
gh <- gh_encode(coords$x, coords$y)
gh

Decoding geohashes

The reverse of encoding geohashes is of course decoding them -- taking a given geohash string and converting it into global coordinates. For example, the Ethiopian coffee growing region of Yirgacheffe is roughly at sc54v:

gh_decode('sc54v')

It can also be helpful to know just how precisely we've identified these coordinates; the include_delta argument gives the cell half-widths in both directions in addition to the cell centroid:

gh_decode('sc54v', include_delta = TRUE)

For more detail on precision, see the table earlier on this vignette which shows the approximate level potential delta at different precision levels.

In terms of latitude and longitude, all geohashes with the same precision have the same dimensions (though the physical size of the "rectangle" changes depending on the latitude); as such it's easy to figure out thecell half-widths from the precision alone using gh_delta:

gh_delta(5L)

You can also pass entire vectors into gh_decode to decode multiple geohashes at once.

gh_decode(gh)

Geohash neighborhoods

One unfortunate consequence of the geohash system is that, while geohashes that are lexicographically similar (e.g. wxyz01 and wxyz12) are certainly close to one another, the converse is not true -- for example, 7gxyru and k58n2h are neighbors! Put another way, small movements on the globe occasionally have visually huge jumps in the geohash-encoded output.

The gh_neighbors function is designed to address this. Calling gh_neighbours will return all of the geohashes adjacent to a given geohash (or vector of geohashes) at the same level of precision.

For example, the Merlion statue in Singapore is roughly at w21z74nz, but this level of precision zooms in a bit too far. The geohash neighborhood thereof can be found with:

gh_neighbors('w21z74nz')

Working with spatial formats

There are several geohashTools helper functions for converting geohashes into sp and sf class objects.

gh_to_*

The gh_to_sp, gh_to_spdf and gh_to_sf functions convert geohash or vector of geohashes into spatial objects of sp and sf class respectively. Consider the previous example with the Singapore Merlion and some random data.

library(sf)
merlion_ghs <- gh_neighbors('w21z74')
merlion_nbhd <- gh_to_sf(merlion_ghs)

# Example plot of geohashes neighbouring w21z74
plot(merlion_nbhd, col = NA, reset = FALSE, key.pos = NULL)
text(
  st_coordinates(st_centroid(merlion_nbhd)),
  labels = row.names(merlion_nbhd)
)

gh_covering

Sometimes we have a set of spatial datapoints that we want to aggregate or index using geohashes. The gh_covering function produces a grid of geohashes that overlap with the spatial points. For this function to work, the spatial object has to be in the WGS84 (EPSG 4326) coordinate reference system that the geohashes use.

Let's use the included meuse dataset in the sp package which shows point coordinates for metal deposits. The first element in this object shows values for cadmium deposits.

if (!requireNamespace('ggplot2', quietly = TRUE)) {
  install.packages('ggplot2')
}

library(ggplot2)

data(meuse, package = 'sp')
meuse_sf = st_as_sf(meuse, coords = c('x', 'y'), crs = 28992L, agr = 'constant')
meuse_sf <- st_transform(meuse_sf, crs = 4326L)

ggplot() +
  geom_sf(data = meuse_sf, aes(colour = cadmium)) +
  theme_void()

By default, gh_covering creates a grid that covers the extent of the bounding box for the spatial object.

meuse_gh <- gh_covering(meuse_sf)
meuse_gh

We can visualize what this looks like with geohashes overlayed on top.

ggplot() +
  geom_sf(data = meuse_sf) +
  geom_sf(data = meuse_gh, fill  = NA) +
  geom_sf_text(data = meuse_gh, aes(label = rownames(meuse_gh))) +
  theme_void()

Alternatively, using gh_covering with the parameter minimal = TRUE will only create geohashes for intersecting objects.

meuse_gh <- gh_covering(meuse_sf, minimal = TRUE)

ggplot() +
  geom_sf(data = meuse_sf) +
  geom_sf(data = meuse_gh, fill  = NA) +
  geom_sf_text(data = meuse_gh, aes(label = rownames(meuse_gh))) +
  theme_void()