library(raster); library(rgdal); library(dplyr); library(yarrr)

# saving projections/datums for later use
LongLatCRS <-'+init=epsg:4326'
albersCRS <- '+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80'

################################### Prepare shapefile of US states

# read in US states shapefile
b0 <- rgdal::readOGR("D:\\Desktop Files November 2019\\fireflies\\us state bounds\\us_states\\us_states.shp")
class(b0) # sp
crs(b1) # CRS arguments: +proj=longlat +datum=WGS84 +no_defs 

# plot original file
plot(b0)

# subsetting so we just have the continental US
b1 <- b0[b0$name != 'Alaska',]
b1 <- b1[b1$name != 'Hawaii',]
b1 <- b1[b1$name != 'Guam',]
b1 <- b1[b1$name != 'Northern Mariana Islands',]
b1 <- b1[b1$name != 'American Samoa',]
b1 <- b1[b1$name != 'United States Virgin Islands',]
b1 <- b1[b1$name != 'Puerto Rico',]

# plot
plot(b1) # looks nice!


################################### Prepare shapefile of Great Lakes (bc US states looks dumb without them)

# read in great lakes shapefile
gl0 <- rgdal::readOGR("D:\\Working Files\\Main Storage\\ArcGIS\\Vector Shapefiles\\Great_Lakes (1)\\Great_Lakes.shp")
class(gl0) # sp
gl1 <- spTransform(gl0, CRS(LongLatCRS)) # transform
plot(gl1, add = TRUE, col = "lightblue")

################################### Shapefiles for all trees!!!

# https://www.fs.fed.us/nrs/atlas/littlefia/species_table.html#

################################### Prepare shapefile of American Sycamore

# figuring out the damn projection 

### https://www.fs.fed.us/nrs/atlas/littlefia/albers_prj.txt

#COORDINATE SYSTEM DESCRIPTION
#
#Projection               ALBERS
#Units                    METERS             
#Spheroid             CLARKE1866
#Parameters:
#1st standard parallel                                   38  0  0.000
#2nd standard parallel                                   42  0  0.000
#central meridian                                       -82  0  0.000
#latitude of projection's origin                         40  0  0.000
#false easting (meters)                                       0.00000
#false northing (meters)                                      0.00000

# this site was helpful in interpreting the +proj4 string
# https://proj.org/operations/projections/aea.html

LittleProjection <- '+proj=aea +lat_1=38 +lat_2=42 +lat_0=40 +lon_0=-82 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD27 +units=m +no_defs'

# read in shapefile
amsyc0 <- rgdal::readOGR("C:\\Users\\obses\\Desktop\\sycamore\\litt731av.shp")
class(amsyc0) # sp
plot(amsyc0) # looks fine

# define CRS
crs(amsyc0) <- LittleProjection 
plot(amsyc0) # still looks fine

# reproject
amsyc1 <- spTransform(amsyc0, LongLatCRS) # transform

extent(amsyc1)
extent(b1)

plot(b1)
plot(amsyc1, add = TRUE, col = yarrr::transparent(orig.col = "cornflowerblue", trans.val = 0.65, maxColorValue = 255))

################################### Prepare shapefile of Loblolly Pine

# read in shapefile
lobpin0 <- rgdal::readOGR("C:\\Users\\obses\\Desktop\\loblolly\\litt131av.shp") # read shapefile
crs(lobpin0) <- LittleProjection # define projection
lobpin1 <- spTransform(lobpin0, LongLatCRS) # transform
plot(lobpin1, add = TRUE, col = yarrr::transparent(orig.col = "yellow", trans.val = 0.65, maxColorValue = 255))
plot(gl1, add = TRUE, col = "lightblue")

#

clip1 <- raster::intersect(amsyc1, lobpin1)
plot(clip1, add = TRUE, col = yarrr::transparent(orig.col = "darkolivegreen3", trans.val = 0.4, maxColorValue = 255))

################################### export shapefiles

writeOGR(b1, "C:\\Users\\obses\\Desktop\\forcody\\us_states.shp", driver="ESRI Shapefile", layer='test') # us_states
writeOGR(gl1, "C:\\Users\\obses\\Desktop\\forcody\\greatlakes.shp", driver="ESRI Shapefile", layer='test') # great lakes
writeOGR(amsyc1, "C:\\Users\\obses\\Desktop\\forcody\\sycamore.shp", driver="ESRI Shapefile", layer='test') # american sycamore range
writeOGR(lobpin1, "C:\\Users\\obses\\Desktop\\forcody\\loblolly.shp", driver="ESRI Shapefile", layer='test') # loblolly pine range
writeOGR(clip1, "C:\\Users\\obses\\Desktop\\forcody\\overlap.shp", driver="ESRI Shapefile", layer='test') # overlap between the two species

################################### export ORIGINAL shapefiles

writeOGR(b0, "C:\\Users\\obses\\Desktop\\spatial_data_weebly\\us_states.shp", driver="ESRI Shapefile", layer='test') # us_states

writeOGR(gl1, "C:\\Users\\obses\\Desktop\\forcody\\greatlakes.shp", driver="ESRI Shapefile", layer='test') # great lakes
writeOGR(amsyc1, "C:\\Users\\obses\\Desktop\\forcody\\sycamore.shp", driver="ESRI Shapefile", layer='test') # american sycamore range
writeOGR(lobpin1, "C:\\Users\\obses\\Desktop\\forcody\\loblolly.shp", driver="ESRI Shapefile", layer='test') # loblolly pine range
writeOGR(clip1, "C:\\Users\\obses\\Desktop\\forcody\\overlap.shp", driver="ESRI Shapefile", layer='test') # overlap between the two species