# load libraries
library(sf)

colors <- list(darkblue = "#062e57",
               lightblue = "#b1dcf3")


# save figure
tiff(filename = "figures/map.tiff",
     width = 11,
     height = 17,
     units = "in",
     res = 600,
     compression = "lzw")

png(filename = "figures/map.png",
    width = 11,
    height = 17,
    units = "in",
    res = 600)

# plot map

plot(sf::st_geometry(extent_poly), 
     col = colors$darkblue, 
     extent = extent_bbox,
     xlim = c(extent$longitude_min, extent$longitude_max),
     ylim = c(extent$latitude_min, extent$latitude_max),
     border = "black")
plot(sf::st_geometry(data$NHDWaterbody), 
     col = colors$lightblue, 
     border = NA,
     extent = extent_bbox,
     xlim = c(extent$longitude_min, extent$longitude_max),
     ylim = c(extent$latitude_min, extent$latitude_max),
     add = TRUE)
plot(sf::st_geometry(data$NHDArea), 
     col = colors$lightblue, 
     border = NA,
     extent = extent_bbox,
     xlim = c(extent$longitude_min, extent$longitude_max),
     ylim = c(extent$latitude_min, extent$latitude_max),
     add = TRUE)
plot(sf::st_geometry(data$NHDFlowline), 
     col = colors$lightblue,
     lwd = data$NHDFlowline$TotDASqKM^0.3204*0.0446,
     border = NA,
     extent = extent_bbox,
     xlim = c(extent$longitude_min, extent$longitude_max),
     ylim = c(extent$latitude_min, extent$latitude_max),
     add = TRUE)
# finish saving figure
dev.off()