# Code written by Marisa Szubryt 2024/01/07, modified from links below

# https://cran.r-project.org/web/packages/fipio/readme/README.html
# https://rdrr.io/cran/fipio/f/README.md
# https://www.reddit.com/r/Rlanguage/comments/se5w7b/how_to_replace_comma_that_is_not_followed_by_a/ 
# https://search.r-project.org/CRAN/refmans/totalcensus/html/convert_fips_to_names.html 
# https://search.r-project.org/CRAN/refmans/usdata/html/abbr2state.html
# https://stackoverflow.com/questions/13093931/remove-last-word-from-string 
# https://stackoverflow.com/questions/72408598/r-extract-first-two-characters-from-a-column-in-a-dataframe 
# https://stackoverflow.com/questions/18115550/combine-two-or-more-columns-in-a-dataframe-into-a-new-column-with-a-new-name 

# Set working directory and clear environment
setwd("/Users/mbs/Desktop/Heterotheca/Heterotheca-RCode/Coordinate_county_maps")
rm(list=ls()); graphics.off(); gc()

# Load in map-making packages for coordinate to county conversion
install.packages('dyplr')
install.packages('fipio')
install.packages('stringr')
install.packages('totalcensus')
install.packages('usdata')
library(dplyr)
library(fipio)
library(stringr)
library(totalcensus)
library(usdata)

# Read in CSV with coordinates
# The below CSV works
points <- read.csv('example_coordinates_simple.csv'); attach(points); points; names(points)

# The below CSV doesn't work 
# points <- read.csv('coordinates.csv'); attach(points); points; names(points)

# Add in necessary columns
points$fips<-'';points$state_fips<-'';points$county_fips<-'';points$state_name<-'';points$county_name<-''; points$SC_name<-''

# Extract latitude and longitude
points_longitude <- points$longitude; points_longitude <- as.numeric(points_longitude); points_longitude
points_latitude <- points$latitude; points_latitude <- as.numeric(points_latitude); points_latitude

# Add FIPS codes (2 numbers for state, 3 numbers for county)
fips_codes <- c()
for (i in c(1:length(points_longitude))) {
  fips_iteration <- fipio::coords_to_fips(x = data.frame(points_longitude[i], points_latitude[i], coords = c("X", "Y")))
  fips_codes <- c(fips_codes, fips_iteration)}
fips_codes <- fips_codes [!is.na(fips_codes)]; fips_codes
points$fips <- fips_codes; points$fips <- as.character(points$fips); points$fips; points

# Split 5 number FIPS codes into state and county fip codes
points$state_fips<- substr(points$fips,1,2); fips_states <- points$state_fips; fips_states <- as.character(fips_states); fips_states
points$county_fips<- substr(points$fips,3,5); fips_counties <- points$county_fips; fips_counties <- as.character(fips_counties); fips_counties

# Generate state and county names from FIPS codes 
states_names <- convert_fips_to_names(fips_states); states_names
counties_names <- convert_fips_to_names(fips_counties, states=states_names, geo_header='COUNTY'); counties_names

# Add state/county names into dataframe
points$state_name <- states_names; points$county_name <- counties_names; points

# Convert state abbreviations to full name
points$state_name <- abbr2state(points$state_name); points

# Remove the 'County' portion of the counties_names vector
points <- points %>% mutate(county_name=word(county_name, 1, -2)); points

# Convert state and county to lower case
points$state_name <- tolower(points$state_name); points$county_name <- tolower(points$county_name); points

# Convert the state and county names into a format usable by the county-highlighting code
points$SC_name <- paste(points$state_name, ",", points$county_name); points
points$SC_name <- gsub(" , (?! )", ",", points$SC_name, perl=TRUE); points

# Extract out the list of "state,county" to insert into the county-highlighting code
plantCounties <- points$SC_name; plantCounties

# Export the amended dataframe for reference
write.csv(points, "example_coordinates_simple_amended.csv", row.names=FALSE)

###########################################################################

# Code written by Patrick Alexander (2024/01/02), modified by Marisa Szubryt

# Load in map-making packages
install.packages('ggplot2')
install.packages('maps')
install.packages('sf')
install.packages('tidyverse')
library(ggplot2)
library(sf)
library(maps)
library(tidyverse)

# Run code to automatically set up map-making
states <- st_as_sf(maps::map("state", plot = FALSE, fill = TRUE))
counties <- st_as_sf(maps::map("county", plot = FALSE, fill = TRUE))
plantStates <- substring(plantCounties, 1, regexpr(",", plantCounties) - 1) %>% unique()
counties <- counties %>% 
  mutate(state=substring(ID, 1, regexpr(',', ID)-1)) %>%
  filter(state %in% plantStates)
plantStates <- states %>% filter(ID %in% plantStates)
plantCounties <- counties %>% filter(ID %in% plantCounties)
bbox <- st_bbox(plantStates)
zone <- floor(((mean(c(bbox$xmin,bbox$xmax))+180)/6) +1)
projBbox <- st_bbox(st_transform(plantStates, crs = paste('+proj=utm +zone=',zone,' +datum=WGS84', sep='')))

# Can change thickness/color of state/county bounaries and fill color of selected counties
ggplot(NULL) + geom_sf(data=counties, lwd=0.33, color='black', fill='white') +
  geom_sf(data=plantCounties, lwd=0.33, fill='steelblue') +
  geom_sf(data = states, lwd=0.75, color='black', fill='transparent') +
  coord_sf(crs=paste('+proj=utm +zone=',zone,' +datum=WGS84',sep=''),
           xlim=c(projBbox$xmin, projBbox$xmax), ylim=c(projBbox$ymin, projBbox$ymax), expand=TRUE) +
  theme_void()

fipio::as_fips(state = "FAKE")                                                                                                                                                            
#> Error in x[is.na(x)] <- y[!is.na(y)] : replacement has length zero