Hexagon_state_map.R

library(tidyverse)
library(tigris)
library(maptools)
library(rgeos)
library(leaflet)

# Number of Hexagons
n_hex <- 3300

# Get state map
#mn_map  <- tigris::block_groups(state = c("MN"), cb = T)
#mn_data <- mn_map@data

# Remove Lake Superior block groups
#mn_map  <- subset(mn_map, !GEOID %in% c("271379901000", "270759901000", "270319900000"))
#plot(subset(mn_map, GEOID == "270314801002"))
#plot(mn_map)

# County map
mn_map <- tigris::counties(state = c("MN"), cb = T)

plot(mn_map)

mn_data <- mn_map@data

# Convert to UTM
utm_proj <-  CRS("+init=epsg:3857") # UTM: CRS("+init=epsg:26915")

lat_long_proj <- proj4string(mn_map) #"+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0"

mn_map <- spTransform(mn_map, utm_proj)
plot(mn_map)

# Extract all MN coordinates
mn_df <- fortify(mn_map)

# Find state boundaries
mn_bounds <- summarize(mn_df, 
                       min_lat  = min(lat),
                       max_lat  = max(lat),
                       min_long = min(long),
                       max_long = max(long))

# MN length to width ratio
mn_ratio <- with(mn_df, diff(range(lat)) / diff(range(long)))

# Hex size
hex_size <- sqrt(with(mn_df, diff(range(lat)) * diff(range(long))) / n_hex)

# Hex longitude
hex_long <- ceiling(diff(range(mn_df$long)) / hex_size) + 1

# Hex latitude
hex_lat <- ceiling(diff(range(mn_df$lat)) / hex_size) + 2
  
# Create hex grid
hex_grid <- data_frame(lat  = as.numeric(sapply(0:hex_lat, function(i) rep(max(mn_df$lat) - ((i-1) * hex_size ), hex_long + 1))),
                       long = rep(as.numeric(sapply(0:hex_long, function(i) min(mn_df$long) + ((i-1) * hex_size))), hex_lat + 1))


# Adjust every other row for hexagon fitting
hex_grid$adjust <- rep(sapply(1:(hex_lat+1), function(i) i %% 2 == 0),  each = hex_long + 1)

hex_grid$long <- hex_grid$long + hex_size * 0.5 * hex_grid$adjust

hex_grid$point_id <- 1:nrow(hex_grid)

coordinates(hex_grid) <- ~long + lat

proj4string(hex_grid) <- utm_proj     # CRS("+init=epsg:26915")

plot(mn_map)

plot(hex_grid, col = "red", add = T)


# Remove points further than 1 hexagon away from MN
#hex_buffer <- gBuffer(hex_grid, width = hex_size/2, byid = T)
#hex_buffer <- crop(hex_buffer, mn_map)
#plot(mn_map)
#plot(hex_buffer, col = "blue", add = T)
#plot(hex_grid, add =T)
#hex_grid <- subset(hex_grid, point_id %in% hex_buffer$point_id)
#plot(hex_grid, col = "red", add =T)


# Create Hexagon grid
hex_pts <- spsample(hex_grid, type = "hexagonal", cellsize = hex_size, n = length(hex_grid)) 


# Convert to lat/long
hex_grid <- spTransform(hex_grid,  CRS("+init=epsg:4326")) # CRS("+init=epsg:4326"))

# Get coordinates
hex_df <- data_frame(lat = hex_grid@coords[,'lat'], long = hex_grid@coords[,'long'])

# Save map
write_csv(hex_df, "C://Users//dkvale//Desktop//mn_hex_map.csv")


# ------------ #
# Leaflet test #
# ------------ #

hex_pts <- HexPoints2SpatialPolygons(hex_pts)

hex_pts <- hex_pts[mn_map, ]

a <- over(hex_pts, mn_map)

hex_pts$county <- a$NAME

hex_pts$aqi_color <- cut(as.numeric(as.factor(hex_pts$county)), breaks = 3, labels = c("#9BF59B", "#ffff00", "#ff7e00")) 

plot(hex_pts, col = "blue")
plot(mn_map, add = T)

hex_pts <- spTransform(hex_pts,  CRS("+init=epsg:4326")) # CRS("+init=epsg:4326"))

h_map <- leaflet(hex_pts) %>% 
         addProviderTiles("CartoDB.Positron", options = providerTileOptions(opacity = 1)) %>% 
         addProviderTiles("Stamen.TonerLines", options = providerTileOptions(opacity = 0.2)) %>%
         addPolygons(stroke       = T,
                     color        = "#d3d3d3",
                     weight       = 0.2,
                     smoothFactor = 0.3,
                     fillOpacity  = 0.5,
                     fillColor    = ~aqi_color, # sample(c("#9BF59B", "#ffff00", "#ff7e00"), 1085, replace = T),
                     popup        = paste0("<b><span style='font-size:1.5em;'>", signif(5, 2), "</span> </b><br>",
                                           gsub("_", " ", "pancake"), "<br><br>",
                                           hex_pts$county, ": ", "hex_pts$point_id", " "))

h_map