# library(rgdal)
# library(sp)
library(sf)
# library(rgeos)

#
# Load data
#

# State boundaries
statefp <- "data/cb_states/cb_2022_us_state_20m.shp"
# statebnds <- readOGR(statefp)
statebnds <- st_read(statefp)

# Conterminous United States
inconterm <- !(statebnds$STUSPS %in% c("AK", "HI", "PR"))
conterm <- statebnds[inconterm,]

# Store locations
stores <- read.csv("data/walmart-openings-geocoded.csv", 
                   sep=",")
walmarts <- stores[stores$store == "walmart",]


#
# Draw boundaries
#

# States and United States territories
plot(st_geometry(statebnds))

# Conterminous United States
plot(st_geometry(conterm))



#
# Draw locations
#

par(mar = c(0, 0, 0, 0))
plot(st_geometry(conterm))
points(walmarts$lng, walmarts$lat, 
       col="blue")



#
# Draw scaled symbols
#

# Count locations by state
wal <- read.csv("data/walmart_addresses.csv", 
                header=FALSE)
statecnts <- table(wal$V4)

# Centroids
# state_centroids <- gCentroid(conterm, byid=TRUE)
state_centroids <- st_centroid(conterm)

# Merge
statesdf <- data.frame(abbrev = conterm$STUSPS)
# statesdf <- cbind(statesdf, state_centroids@coords)
statesdf <- cbind(statesdf, 
                  st_coordinates(state_centroids))
imatch <- match(statesdf$abbrev, dimnames(statecnts)[[1]])
statesdf$walmarts <- as.numeric(statecnts)[imatch]

# Symbols
par(mar=c(0,0,0,0))
plot(st_geometry(conterm), border="#cccccc")
symbols(statesdf$X, statesdf$Y, 
        circles = sqrt(statesdf$walmarts)/7,
        bg = "#b7c1e2",
        inches = FALSE,
        add = TRUE)




#
# Lines
#

par(mar=c(0,0,0,0))
plot(st_geometry(conterm), border="#cccccc")
lines(x = statesdf$X[1:2],
      y = statesdf$Y[1:2],
      lwd = 2, 
      col = "#30437b")
points(x = statesdf$X[1:2],
       y = statesdf$Y[1:2],
       pch = 21,
       bg = "#30437b",
       col = "#000000")
text(x = statesdf$X[1:2],
     y = statesdf$Y[1:2]+.75, 
     labels = statesdf$abbrev[1:2])


# A longer path
roadtrip <- read.csv("data/brewery-road-trip-path.tsv", 
                     sep = "\t")
par(mar=c(0,0,0,0))
plot(st_geometry(conterm), border="#cccccc")
lines(roadtrip$lng, roadtrip$lat, lwd = 3)



#
# Projections
#

# Project map
conterm_p4s <- 
    "+proj=aea +lat_1=29.5 +lat_2=45.5 +lon_0=97.2w"
# conterm_albers <- spTransform(conterm, CRS(conterm_p4s))
conterm_albers <- st_transform(conterm, conterm_p4s)

# Project points
# walsp <- walmarts
# coordinates(walsp) <- c("lng", "lat")
# proj4string(walsp) <- CRS("+proj=longlat")
# walmarts_albers <- spTransform(walsp, CRS(conterm_p4s))
wal_sf <- st_as_sf(walmarts, 
                   coords = c("lng", "lat"), 
                   crs = "+proj=longlat")
wal_sf_albers <- st_transform(wal_sf, conterm_p4s)

# Projected points on map
par(mar = c(0,0,0,0))
plot(st_geometry(conterm_albers), 
     border = "#cccccc", lwd = .5)
points(st_coordinates(wal_sf_albers), 
       col="#30437b", cex = .5)


# Projected scaled symbols
statesdf_sf <- st_as_sf(statesdf,
                        coords = c("X", "Y"),
                        crs = "+proj=longlat")
statesdf_sf_albers <- st_transform(statesdf_sf, conterm_p4s)
par(mar=c(0,0,0,0))
plot(st_geometry(conterm_albers), border="#cccccc")
symbols(st_coordinates(statesdf_sf_albers), 
        circles = sqrt(statesdf$walmarts)*12000,
        bg = "#b7c1e2",
        inches = FALSE,
        add = TRUE)


# Projected road trip
roadtrip_sf <- st_as_sf(roadtrip,
                        coords = c("lng", "lat"),
                        crs = "+proj=longlat")
roadtrip_sf_albers <- st_transform(roadtrip_sf, conterm_p4s)

# Map the projected road trip
par(mar=c(0,0,0,0))
plot(st_geometry(conterm_albers), border="#cccccc")
lines(st_coordinates(roadtrip_sf_albers), 
      lwd = 3)