I've been playing with st_nearest_feature in R/sf to calculate nearest features on a map.

ct <- counties(cb=T, resolution = '20m') %>% st_transform(5070)
st <- states(cb=T, resolution = '20m') %>% filter(NAME %in% c('Texas','New York')) %>% st_transform(5070)

nr <- st_nearest_feature(ct, st)

ct <- cbind(ct, st_drop_geometry(st)[nr,])