library(tidyverse)

geomid <- function(lat, lon)
{
  # http://www.geomidpoint.com/calculation.html
  lat <- lat * pi / 180
  lon <- lon * pi / 180
  x <- mean(sum(cos(lat) * cos(lon)))
  y <- mean(sum(cos(lat) * sin(lon)))
  z <- mean(sum(sin(lat)))
  lat <- atan2(y, x)
  hyp <- sqrt(x * x + y * y)
  lon <- atan2(z, hyp)
  c(lat = lat * 180 / pi, lon = lon * 180 / pi)
}

# calculate the geographic midpoint for each county in the US
map_data("county") %>% 
  mutate(polyname = str_c(region, subregion, sep = ",")) %>%
  group_by(polyname) %>%
  group_modify(~ {
    as.data.frame(t(geomid(.x$lat, .x$long)))
  })