### BASED ON TARAS KADUK'S https://taraskaduk.com/2017/11/26/pixel-maps/

### setup ----------------------------------------------------------------------

# supress warnings/messages
knitr::opts_chunk$set(message = FALSE, warning = FALSE)

# install and/or load packages
pkg <- c("tidyverse", "maps", "ggmap", "here")
sapply(pkg, function(x){
  if(! x %in% installed.packages()) install.packages(x)
  require(x, character.only = TRUE)
})
## tidyverse      maps     ggmap      here 
##      TRUE      TRUE      TRUE      TRUE
### cities data ----------------------------------------------------------------

# cities to geolocate
city_name <- c("breda NL", 
               "utrecht", "rotterdam", "tilburg", "amsterdam",
               "london", "singapore", "kuala lumpur", "zanzibar",
               "antwerp", "middelkerke", "maastricht", "bruges",
               "san fransisco", "vancouver", "willemstad", "hurghada",
               "paris", "rome", "bordeaux", "berlin", 
               "kos", "crete", "kefalonia", "corfu", 
               "dubai", " barcalona", "san sebastian", "dominican republic",
               "porto", "gran canaria", "albufeira", "istanbul", 
               "lake como", "oslo", "riga", "newcastle",
               "dublin")

# initialize dataframe
tibble(
  city = city_name,
  lon = rep(NA, length(city_name)),
  lat = rep(NA, length(city_name))
) ->
  cities

# loop cities through API to overcome SQ limit
for(c in city_name){
  temp <- tibble(lon = NA)
  # geolocate until found or tried 5 times
  attempt <- 0
  while(is.na(temp$lon) & attempt < 5) {
    temp <- geocode(c)
    attempt <- attempt + 1
  } 
  # write to dataframe
  cities[cities$city == c, -1] <- temp
}

head(cities)
## # A tibble: 6 x 3
##   city         lon   lat
##   <chr>      <dbl> <dbl>
## 1 breda NL   4.77   51.6
## 2 utrecht    5.12   52.1
## 3 rotterdam  4.48   51.9
## 4 tilburg    5.06   51.6
## 5 amsterdam  4.90   52.4
## 6 london    -0.128  51.5
### dot data -------------------------------------------------------------------

# generate worldwide dots
lat <- data_frame(lat = seq(-90, 90, by = 1))
lon <- data_frame(lon = seq(-180, 180, by = 1))
dots <- merge(lat, lon, all = TRUE)

# exclude water-based dots
dots %>%
  mutate(country = map.where("world", lon, lat),
         lakes = map.where("lakes", lon, lat)) %>%
  filter(!is.na(country) & is.na(lakes)) %>% 
  select(-lakes) ->
  dots

head(dots)
##   lat  lon                   country
## 1 -83 -173                Antarctica
## 2 -83 -172                Antarctica
## 3 -83 -171                Antarctica
## 4  60 -167 USA:Alaska:Nunivak Island
## 5  60 -166 USA:Alaska:Nunivak Island
## 6  65 -166                USA:Alaska
### visualize ------------------------------------------------------------------

# plot the data
dots %>% ggplot(aes(x = lon, y = lat)) + 
  geom_point(col = "black", size = 0.25) +
  geom_point(data = cities, col = "red", size = 3, alpha = 0.7) + 
  theme_void() +
  theme(
    panel.background = element_rect(fill = "#006994"),
    plot.background = element_rect(fill = "#006994")
  ) -> dot_map

# save plot
ggsave(here("worlmap_dots.png"), dot_map, dpi = 600, width = 8, height = 4.5)

dot_map