DataProcessingScriptCFD.R

#Big Thanks To Rebecca!!!!
#Credit her Github https://github.com/HobanLab/QUBR_GenGeoEcoDemoCorr/blob/main/analyses/Data_Processing_Script.R and https://github.com/HobanLab/QUBR_GenGeoEcoDemoCorr/blob/main/analyses/hypothesis_1_clumping_CLEANED.R
#Catherine's land management edits to Rebecca's work.  This code uses the framework Rebecca laid out to analyze Quercus brandegii trees along three riverbed sites.  I am adapting it for use of butternuts in two sites that have been cleared for trail building.

# There are 5 main types of dataframes created in this script:

# 1) "fixed_field_data_processed_sf_trans_coordinates" (the processed spatial/size/shape data for individual trees)
# this dataframe was then filtered for each population, i.e. LM_fixed_field_data_processed
# 2) "all_points_fixed_field_data_processed_terrain" (the processed spatial, tree size/shape, ELEVATION, SLOPE, and ASPECT for individual trees)
# this dataframe was then filtered for each population, i.e. LM_fixed_field_data_processed_terrain
# 3) "LM_fixed_field_data_processed_terrain_dist" (the processed spatial, tree size/shape, elevation, slope, aspect, and DISTANCE TO RIVER for individual trees)
# there is only the three dataframes filtered for each population, i.e. LM_fixed_field_data_processed_terrain_dist
# 4) "LM_fixed_field_data_processed_soils" (the processed spatial, tree size/shape, elevation, slope, aspect, distance to river, and SOIL METRICS for individual trees)
# there is only the three dataframes filtered for each population, i.e. LM_fixed_field_data_processed_soils
# 5) "all_known_pop_soils" (the processed spatial and soil metrics for the 20 known QUBR population points)

# NOTE: Uncomment and run line 273-349, generating the 15m elevation rasters, if they are no longer available, ex: "./data/15 m Elevation Raster/CEM_15_utm_LM.tif" is no longer available or needs to be updated

#### Loading libraries and relevant data ####

library(tidyverse) 
library(moments) # for calculating the moments of each variable
library(sf) # for plotting spatial objects
library(smatr)
library(ggpmisc)
library(PMCMRplus) # for Dunn test
library(geomtextpath) # for PCA graphing
library(spatstat) # to run the Ripley's K function: Kest
library(stars) # for sf_rasterize function
library(raster) #to use crop
library(starsExtra) #to use dist_to_nearest
library(geostatsp) # To successfully use as.im
library(tmaptools)
library(mgcv) #to use GAM function 
library(plotly) #to 3d plot variables
library(MuMIn) #to be able to use dredge
library(visreg) # to be able to plot Aspect/categorical variables with GAM
library(spatialEco) # for heat load index function
library(geosphere) # for finding the distance of populations to the coast

#### Loading and processing relevant data ####

# loading in the tree data (size, elevation, lat/lon, ID, size/shape)

fixed_field_data_processed <- data_ILM #imports the csv created from the 2025 data merged with the data processed in analyzing_morpho_data_cleaned.R

#adding a sequential column, "X," to number each tree and editing column names to make sure they align with the code 
anyDuplicated(fixed_field_data_processed[,c('gps_n','gps_w')])
duplicates <- which(duplicated(fixed_field_data_processed[,c('gps_n','gps_w')]) == FALSE)
fixed_field_data_processed <- fixed_field_data_processed[c(duplicates),]
anyDuplicated(fixed_field_data_processed[,c('gps_n','gps_w')])

# creating the point shapefiles of the tree locations for each population in UTM 18 N

#creating a point shapefile of all points with lat lon coordinates and other attributes in WGS 1984
#sf objects are dataframes with rows representing simple features with attributes and a simple feature geometry list-column (sfc)
fixed_field_data_processed_sf <- st_as_sf(fixed_field_data_processed, 
                                          coords = c("gps_w", "gps_n"), crs = 4326)

#creating a transformed point shapefile with UTM 18 N an equal area projection
fixed_field_data_processed_sf_transformed <- st_transform(fixed_field_data_processed_sf, crs = 26918)

#storing point shapefiles for the trees by population

ILM_fixed_field_data_processed_sf <- fixed_field_data_processed_sf_transformed %>%
  st_as_sf()

#create dataframe with X and Y UTM coordinates

fixed_field_data_processed_sf_trans_coords <- st_coordinates(fixed_field_data_processed_sf_transformed) #creates a dataframe with separate x and y columns from the UTM 18N transformation
fixed_field_data_processed_sf_trans_coordinates <- fixed_field_data_processed_sf_transformed %>%
  cbind(fixed_field_data_processed_sf_trans_coords) #combines the x and y coordinate data frame with the transformed sf dataframe


#### Transformations of tree size/shape variables (log and square root) for linear models ####

#creating columns with transformations: logged DBH
fixed_field_data_processed_sf_trans_coordinates <- fixed_field_data_processed_sf_trans_coordinates %>%
  mutate(dbh_lg = log(dbh_cm))

#creating columns with transformations: square root all of the size/shape variables
fixed_field_data_processed_sf_trans_coordinates <- fixed_field_data_processed_sf_trans_coordinates %>%
  mutate(dbh_sqrt = sqrt(dbh_cm))

#creating columns with transformations: inverse all of the variables
fixed_field_data_processed_sf_trans_coordinates <- fixed_field_data_processed_sf_trans_coordinates %>%
  mutate(dbh_inv = (1/dbh_cm))

# Creating fixed_field_data_processed dataframes for each population with the nearest neighbor columns 

ILM_fixed_field_data_processed <- fixed_field_data_processed_sf_trans_coordinates %>%
  filter(site == "ILM") 



#### Loading in ArcGIS river and Baja California Sur shapefile and storing out polygons for each population ####

setwd("~/land_management_analysis_cfd/trails_all_ILM")
trails_ILM <- st_read("BufferedFeatures.shp")
trails_ILM <- trails_ILM$geometry[1]
plot(trails_ILM)

#changing the coordinate reference system of the river polygons to be equal area projection (UTM 12N), uses meters as distance measurement 
#uses meters as distance measurement
trails_ILM_trans <- st_transform(trails_ILM, crs = 26918)

#ensuring the river outlines are shapefiles for the distance measurements
trails_ILM_trans <- st_sf(geometry = trails_ILM_trans)

# Importing BCS and River Shapefiles 

#turning the Baja California Sur polygon into a shapefile, to be able to visualize the point locations
setwd("~/land_management_analysis_cfd/ILM_boundaries")
ILM_polygon <- read_sf("ILM_boundaries.shp")
ILM_polygon <- st_as_sf(ILM_polygon) #ensures foreign will be an sf object (in this case a multipolygon)
plot(ILM_polygon)

#### Creating buffers and boundaries around the rivers ####

#creating buffers around the rivers

#ILM
trails_buffer_ILM <- st_buffer(trails_ILM, 10) #80 m buffer
ggplot()+  #plotting the river shapefile, the buffer, and the tree points
  geom_sf(data = trails_buffer_ILM)+
  #geom_sf(data = trails_ILM_trans)+
  geom_sf(data = ILM_fixed_field_data_processed_sf)

library(tmap)
tm_shape(ILM_polygon) + 
  tm_polygons(fill = "bisque", col = "bisque") +
  tm_shape(trails_ILM) +
  tm_polygons(fill = "burlywood4", col = "burlywood4") +
  tm_shape(ILM_fixed_field_data_processed_sf) +
  tm_dots(fill = "green4") +
  tm_title("ILM butternuts (green) and ILM trails (brown)") +
  tm_grid() +
  tm_graticules()

#creating bounding boxes for each population

#creating a boundary box for LM with the UTM 12 N min and max lat lon values and then turning it into a simple feature geometry
ILM_fixed_field_data_processed_box <- fixed_field_data_processed_sf_transformed %>%
  filter(site == "ILM") %>%
  st_bbox %>%
  st_as_sfc()

#creating bounding boxes for all of the river shapefiles for each population
ILM_box <- st_bbox(trails_ILM_trans)

#i deleted everything w elevation and aspect rasters.  if i needed that to rasterize we will get there


#### Descriptive Summary ####

#histograms
ggplot(fixed_field_data_processed_sf_trans_coordinates) + #DBH
  geom_histogram(aes(x = dbh_cm))+
  xlab("Diameter at breast height")+
  ylab("Frequency")

#LM
ggplot(ILM_fixed_field_data_processed) + #DBH
  geom_histogram(aes(x = dbh_cm))+
  xlab("Diameter at breast height")+
  ylab("Frequency")


#Summary Statistics
# Create a df which contains the "classical" univariate dist'n stats of each of the important variables
field_data_summarized <- fixed_field_data_processed %>%
  dplyr::select(dbh_cm, height_ft, live_adult_canker_base, live_adult_girdle, live_adult_canker_1stlivebranch, adult_percent_live_canopy, seedling_germ_yr, seedling_girdle, seedling_previous_stem_count) %>%  # Keep only the columns we are interested in getting summary values of
  summarise(across(everything(), list(mean = mean, median = median, var = var, sd = sd), na.rm=TRUE)) # Create columns which summarize the mean, median, variance, and standard deviation of each of the selected columns --> these will be used on the histogram plots
# View(field_data_summarized)

#### Comparing the hand-held collected GPS elevation data to the raster derived elevation data ####
#no rasters!!!
#so i deleted a big fat section here too

#### Creating the distance to river columns ####

#LM

#turning the river polygon into a linestring object and then into a raster, to be able to later calculate the distances
trails_ILM_trans_points <- st_cast(trails_ILM_trans, "MULTILINESTRING") #turning the polyline of the river into a multilinestring object
trails_ILM_trans_point_raster <- st_rasterize(trails_ILM_trans_points) #creating a raster out of the river linestring object
plot(trails_ILM_trans_point_raster,, breaks = "equal") #plotting the river linestring object

#turning the river buffer polygon into a linestring object and then into a raster to be able to later calculate the distances
trails_buffer_ILM_trans <- st_transform(trails_buffer_ILM, crs = 26918)
trails_buffer_ILM_trans <- st_sf(geometry = trails_buffer_ILM_trans)
trails_ILM_buffer_trans_outline <- st_cast(trails_buffer_ILM_trans, "MULTILINESTRING") #turning the polygon of the river buffer into a linestring object
trails_buffer_ILM_point_raster <- st_rasterize(trails_ILM_buffer_trans_outline) #creating a raster of river buffer linestring object
plot(trails_buffer_ILM_point_raster) #plotting the river buffer linestring object

#generating a distance to river raster with the distances of each cell in the buffer raster from the river edge points, whereby the river raster cells are set to a distance of 0 m
trails_buffer_ILM_point_raster[is.na(trails_buffer_ILM_point_raster[])] <- 0  #making sure the cell that are not the river buffer linestring raster have a 0 value
dist_near_trails_buffer_ILM <- dist_to_nearest(trails_buffer_ILM_point_raster, trails_ILM_trans_points, progress = T) #creating a raster of the distances of each cell in the buffer raster to the linestring object of the river polygon, this can take a while to run
plot(dist_near_trails_buffer_ILM) #plotting the distance to river raster


## Making it so the cells in the distance raster within or overlapping with the river raster are assigned 1 

#LM

#Assigning points within and overlapping with the river to be "true"
ILM_points_intersects_trails <- st_intersects(ILM_fixed_field_data_processed, trails_ILM_trans, sparse = F) #creating a list of true or falses for whether points intersect the river shapefiles
ILM_fixed_field_data_processed_intersects_trails <- cbind(ILM_fixed_field_data_processed, ILM_points_intersects_trails) #binding the list of true or falses with the point data
#printing the river polygon and the tree points, colored by whether are or aren't within or overlapping with the river
ggplot()+
  geom_sf(data=trails_ILM_trans)+
  geom_sf(data=ILM_fixed_field_data_processed)+
  geom_sf(data=ILM_fixed_field_data_processed_intersects_trails, aes(color = ILM_points_intersects_trails))


## Extracting distance to river for each tree using the distance to river raster

#LM
ILM_distance_data_pts <- st_extract(dist_near_trails_buffer_ILM, ILM_fixed_field_data_processed) #extracting distance to river for each tree
ILM_fixed_field_data_processed_distance  <- cbind(ILM_fixed_field_data_processed, ILM_distance_data_pts) #binding the distance to river data for each point to the LM point dataframe


## Assigning all points within/overlapping river to distances of 0

#LM
ILM_fixed_field_data_processed_distance <- ILM_fixed_field_data_processed_distance %>% 
  mutate(d = case_when((ILM_fixed_field_data_processed_intersects_trails$ILM_points_intersects_trails == T) ~ 0,  #assigns 0 to points within river
                       (ILM_fixed_field_data_processed_intersects_trails$ILM_points_intersects_trails == F) ~ d)) #to points outside of river, it leaves the original distance value



#################################################################################################
################################################################################################
##################################################################################################
#NOTHING HERE HAS BEEN FIXED YET
#### Load in Soil Metric Data ####

# loading in soil textures raster tifs from CONABIO
clay_05 <- raster(paste0("./data/Soil Grid/clay content/clay content 0-5.tif"))
clay_200 <- raster(paste0("./data/Soil Grid/clay content/clay content 100-200.tif"))
silt_05 <- raster(paste0("./data/Soil Grid/silt/silt 0-5.tif"))
silt_200 <-raster(paste0("./data/Soil Grid/silt/silt 100-200.tif"))
sand_05 <- raster(paste0("./data/Soil Grid/sand/sand 0-5.tif"))
sand_200 <- raster(paste0("./data/Soil Grid/sand/sand 100-200.tif"))
ph_05 <- raster(paste0("./data/Soil Grid/pH/ph_0-5.tif")) #0-5 cm ph
ph_200 <- raster(paste0("./data/Soil Grid/pH/ph_100-200.tif")) #100-200 ph
ocd_05 <- raster(paste0("./data/Soil Grid/organic carbon density/OCD_0-5.tif")) #0-5cm organic carbon density
ocd_200 <- raster(paste0("./data/Soil Grid/organic carbon density/OCR_100-200.tif")) #100-200cm organic carbon density
coarse_frag_05 <- raster(paste0("./data/Soil Grid/coarse fragments/coarse_fragments_0-5.tif")) #0-5 cm coarse fragments
coarse_frag_200 <- raster(paste0("./data/Soil Grid/coarse fragments/coarse_fragments_100-200.tif")) #100-200 cm coarse fragments
cat_ex_cap_05 <-raster(paste0("./data/Soil Grid/cation exchange capacity/Cat_exc_0-5.tif")) #0-5 cm cation exchange capacity
cat_ex_cap_200 <- raster(paste0("./data/Soil Grid/cation exchange capacity/Cat_exc_100-200.tif")) #100-200 cm cation exchange capacity
bulk_dens_05 <- raster(paste0("./data/Soil Grid/bulk density/bulk_density_0-5.tif")) #0-5 cm bulk density
bulk_dens_200 <- raster(paste0("./data/Soil Grid/bulk density/bulk_density_100-200.tif")) #100-200 cm bulk density
vol_wat_10kpa_05 <- raster(paste0("./data/Soil Grid/vol. water content at -10 kPa/vol_water_-10_0-5.tif"))  #0-5 cm -10 kpa volumn water content
vol_wat_10kpa_200 <- raster(paste0("./data/Soil Grid/vol. water content at -10 kPa/vol_water_-10_100-200.tif"))  #100-200 cm -10 kpa volumn water content
vol_wat_33kpa_05 <- raster(paste0("./data/Soil Grid/vol. water content at -33 kPa /vol_water_0-5.tif")) #0-5 cm -33 kpa volumn water content
vol_wat_33kpa_200 <- raster(paste0("./data/Soil Grid/vol. water content at -33 kPa /vol_water_100-200.tif")) #100-200 cm -33 kpa volumn water content
vol_wat_1500kpa_05 <- raster(paste0("./data/Soil Grid/vol. water content at -1500 kPa/vol_water_-1500kPa_0-5.tif"))  #0-5 cm -1500 kpa volumn water content
vol_wat_1500kpa_200 <- raster(paste0("./data/Soil Grid/vol. water content at -1500 kPa/vol_water_-1500_100-200.tif")) #100-200 cm -1500 kpa volumn water content
nitrogen_05 <- raster(paste0("./data/Soil Grid/Nitrogen/nitrogen 0-5.tif")) #0-5 cm nitrogen content 
nitrogen_200 <- raster(paste0("./data/Soil Grid/Nitrogen/nitrogen 100-200.tif")) #100-200 cm nitrogen content 
Soil_Organic_Carbon_05 <- raster(paste0("./data/Soil Grid/Soil Organic Carbon/SOC 0-5.tif")) #0-5 cm soil organic carbon
Soil_Organic_Carbon_200 <- raster(paste0("./data/Soil Grid/Soil Organic Carbon/SOC 100-200.tif")) #100-200 cm soil organic carbon

#project the soil metric rasters to equal area projection (UTM 12N) which uses meters as distance measurement 

clay_05_utm <- projectRaster(clay_05, crs=26912) #converting the 0-5 cm clay raster to utm 12
clay_200_utm <- projectRaster(clay_200, crs=26912) #converting the 90-200 cm clay raster to utm 12
silt_05_utm <- projectRaster(silt_05, crs=26912)
silt_200_utm <- projectRaster(silt_200, crs=26912)
sand_05_utm <- projectRaster(sand_05, crs=26912)
sand_200_utm <- projectRaster(sand_200, crs=26912)
ph_05_utm <- projectRaster(ph_05, crs=26912) 
ph_200_utm <- projectRaster(ph_200, crs=26912) 
ocd_05_utm <- projectRaster(ocd_05, crs=26912)
ocd_200_utm <- projectRaster(ocd_200, crs=26912)
coarse_frag_05_utm <- projectRaster(coarse_frag_05, crs=26912)
coarse_frag_200_utm <- projectRaster(coarse_frag_200, crs=26912)
cat_ex_cap_05_utm <- projectRaster(cat_ex_cap_05, crs=26912)
cat_ex_cap_200_utm <- projectRaster(cat_ex_cap_200, crs=26912)
bulk_dens_05_utm <- projectRaster(bulk_dens_05, crs=26912)
bulk_dens_200_utm <- projectRaster(bulk_dens_200, crs=26912)
vol_wat_10kpa_05_utm <- projectRaster(vol_wat_10kpa_05, crs=26912)
vol_wat_10kpa_200_utm <- projectRaster(vol_wat_10kpa_200, crs=26912)
vol_wat_33kpa_05_utm <- projectRaster(vol_wat_33kpa_05, crs=26912)
vol_wat_33kpa_200_utm <- projectRaster(vol_wat_33kpa_200, crs=26912)
vol_wat_1500kpa_05_utm <- projectRaster(vol_wat_1500kpa_05, crs=26912)
vol_wat_1500kpa_200_utm <- projectRaster(vol_wat_1500kpa_200, crs=26912)
nitrogen_05_utm <- projectRaster(nitrogen_05, crs=26912)
nitrogen_200_utm <- projectRaster(nitrogen_200, crs=26912)
Soil_Organic_Carbon_05_utm <- projectRaster(Soil_Organic_Carbon_05, crs=26912)
Soil_Organic_Carbon_200_utm <- projectRaster(Soil_Organic_Carbon_200, crs=26912)
sandy_avail_water_0.5_utm <- vol_wat_33kpa_05_utm - vol_wat_1500kpa_05_utm # Sand Available Water 0-5 cm
sandy_avail_water_100.200_utm <- vol_wat_33kpa_200_utm - vol_wat_1500kpa_200_utm # Sand Available Water 100-200 cm
clay_loam_avail_water_0.5_utm <- vol_wat_10kpa_05_utm - vol_wat_1500kpa_05_utm # Clay/Loam Available Water 0-5 cm
clay_loam_avail_water_100.200_utm <- vol_wat_10kpa_200_utm - vol_wat_1500kpa_200_utm # Clay/Loam Available Water 100-200 cm

nitrogen_200_utm <- projectRaster(nitrogen_200, crs=26912)

"Elevation", "Slope", "Eastness", "Northness", "TWI", "HLI", "Distance to River"

## cropping the soil metrics to bounding boxes around each population ##

# LM

# using the extent of the box around the rivers to crop the raster for each soil texture layer
clay_05_LM <- crop(clay_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100)) 
clay_200_LM <- crop(clay_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
silt_05_LM <- crop(silt_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
silt_200_LM <- crop(silt_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
sand_05_LM <- crop(sand_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
sand_200_LM <- crop(sand_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
ph_05_LM <- crop(ph_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
ph_200_LM <- crop(ph_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
ocd_05_LM <- crop(ocd_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
ocd_200_LM <- crop(ocd_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
coarse_frag_05_LM <- crop(coarse_frag_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
coarse_frag_200_LM <- crop(coarse_frag_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
cat_ex_cap_05_LM <- crop(cat_ex_cap_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
cat_ex_cap_200_LM <- crop(cat_ex_cap_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
bulk_dens_05_LM <- crop(bulk_dens_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
bulk_dens_200_LM <- crop(bulk_dens_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
vol_wat_10kpa_05_LM <- crop(vol_wat_10kpa_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
vol_wat_10kpa_200_LM <- crop(vol_wat_10kpa_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
vol_wat_33kpa_05_LM <- crop(vol_wat_33kpa_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
vol_wat_33kpa_200_LM <- crop(vol_wat_33kpa_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
vol_wat_1500kpa_05_LM <- crop(vol_wat_1500kpa_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
vol_wat_1500kpa_200_LM <- crop(vol_wat_1500kpa_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
nitrogen_05_LM <-  crop(nitrogen_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
nitrogen_200_LM <- crop(nitrogen_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
Soil_Organic_Carbon_05_LM <- crop(Soil_Organic_Carbon_05_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))
Soil_Organic_Carbon_200_LM <- crop(Soil_Organic_Carbon_200_utm, extent(LM_box[1]-100, LM_box[3]+100, LM_box[2]-100, LM_box[4]+100))


# plotting an example soil metric raster with river shapefile and tree points  
ggplot()+
  geom_raster(data = as.data.frame(Soil_Organic_Carbon_05_LM, xy=T), aes(x=x, y=y, fill = SOC.0.5))+
  geom_sf(data = river_LM_trans)+
  geom_sf(data = LM_fixed_field_data_processed)

#creating stacks of the soil raster layers (using multiple for visualization purposes)
soil_stack_LM_soil_text <- stack(clay_05_LM, clay_200_LM, silt_05_LM, silt_200_LM, sand_05_LM, sand_200_LM) #the stack of all of the soil texture rasters
soil_stack_LM_other <- stack(ph_05_LM, ph_200_LM, ocd_05_LM, ocd_200_LM, coarse_frag_05_LM, coarse_frag_200_LM, #the stack of all of the other soil variables, with different extents than the soil texture rasters
                             cat_ex_cap_05_LM, cat_ex_cap_200_LM, bulk_dens_05_LM, bulk_dens_200_LM, vol_wat_10kpa_05_LM,
                             vol_wat_10kpa_200_LM, vol_wat_33kpa_05_LM, vol_wat_33kpa_200_LM, vol_wat_1500kpa_05_LM, 
                             vol_wat_1500kpa_200_LM) 
soil_stack_LM_extra <- stack(nitrogen_05_LM, nitrogen_200_LM, Soil_Organic_Carbon_05_LM, Soil_Organic_Carbon_200_LM)



#plotting the stacked rasters
plot(soil_stack_LM_soil_text) #version with soil textures
plot(soil_stack_LM_soil_text, zlim = c(100, 710)) #version where the plots have the same scale
plot(soil_stack_LM_other) #version with other variables
plot(soil_stack_LM_other, zlim = c(30, 360)) #version where the plots have the same scale
plot(soil_stack_LM_extra) #version with other variables
plot(soil_stack_LM_extra, zlim = c(30, 360)) #version where the plots have the same scale


#LC

# using the extent of the box around the rivers to crop the raster for each soil texture layer
clay_05_LC <- crop(clay_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100)) 
clay_200_LC <- crop(clay_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
silt_05_LC <- crop(silt_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
silt_200_LC <- crop(silt_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
sand_05_LC <- crop(sand_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
sand_200_LC <- crop(sand_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
ph_05_LC <- crop(ph_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
ph_200_LC <- crop(ph_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
ocd_05_LC <- crop(ocd_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
ocd_200_LC <- crop(ocd_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
coarse_frag_05_LC <- crop(coarse_frag_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
coarse_frag_200_LC <- crop(coarse_frag_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
cat_ex_cap_05_LC <- crop(cat_ex_cap_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
cat_ex_cap_200_LC <- crop(cat_ex_cap_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
bulk_dens_05_LC <- crop(bulk_dens_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
bulk_dens_200_LC <- crop(bulk_dens_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
vol_wat_10kpa_05_LC <- crop(vol_wat_10kpa_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
vol_wat_10kpa_200_LC <- crop(vol_wat_10kpa_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
vol_wat_33kpa_05_LC <- crop(vol_wat_33kpa_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
vol_wat_33kpa_200_LC <- crop(vol_wat_33kpa_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
vol_wat_1500kpa_05_LC <- crop(vol_wat_1500kpa_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
vol_wat_1500kpa_200_LC <- crop(vol_wat_1500kpa_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
nitrogen_05_LC <-  crop(nitrogen_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
nitrogen_200_LC <- crop(nitrogen_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
Soil_Organic_Carbon_05_LC <- crop(Soil_Organic_Carbon_05_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))
Soil_Organic_Carbon_200_LC <- crop(Soil_Organic_Carbon_200_utm, extent(LC_box[1]-100, LC_box[3]+100, LC_box[2]-100, LC_box[4]+100))

#creating stacks of the soil raster layers (using multiple for visualization purposes)
soil_stack_LC_soil_text <- stack(clay_05_LC, clay_200_LC, silt_05_LC, silt_200_LC, sand_05_LC, sand_200_LC) #the stack of all of the soil texture rasters
soil_stack_LC_other <- stack(ph_05_LC, ph_200_LC, ocd_05_LC, ocd_200_LC, coarse_frag_05_LC, coarse_frag_200_LC, #the stack of all of the other soil variables, with different extents than the soil texture rasters
                             cat_ex_cap_05_LC, cat_ex_cap_200_LC, bulk_dens_05_LC, bulk_dens_200_LC, vol_wat_10kpa_05_LC,
                             vol_wat_10kpa_200_LC, vol_wat_33kpa_05_LC, vol_wat_33kpa_200_LC, vol_wat_1500kpa_05_LC, 
                             vol_wat_1500kpa_200_LC) 
soil_stack_LC_extra <- stack(nitrogen_05_LC, nitrogen_200_LC, Soil_Organic_Carbon_05_LC, Soil_Organic_Carbon_200_LC)

#plotting the stacked rasters
plot(soil_stack_LC_soil_text) #version with soil textures
plot(soil_stack_LC_soil_text, zlim = c(100, 710)) #version where the plots have the same scale
plot(soil_stack_LC_other) #version with other variables
plot(soil_stack_LC_other, zlim = c(30, 360)) #version where the plots have the same scale
plot(soil_stack_LC_extra) #version with other variables
plot(soil_stack_LC_extra, zlim = c(30, 180)) #version where the plots have the same scale


#SD

#using the extent of the box around the rivers to crop the raster for each soil texture layer
clay_05_SD <- crop(clay_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100)) 
clay_200_SD <- crop(clay_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
silt_05_SD <- crop(silt_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
silt_200_SD <- crop(silt_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
sand_05_SD <- crop(sand_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
sand_200_SD <- crop(sand_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
ph_05_SD <- crop(ph_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
ph_200_SD <- crop(ph_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
ocd_05_SD <- crop(ocd_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
ocd_200_SD <- crop(ocd_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
coarse_frag_05_SD <- crop(coarse_frag_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
coarse_frag_200_SD <- crop(coarse_frag_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
cat_ex_cap_05_SD <- crop(cat_ex_cap_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
cat_ex_cap_200_SD <- crop(cat_ex_cap_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
bulk_dens_05_SD <- crop(bulk_dens_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
bulk_dens_200_SD <- crop(bulk_dens_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
vol_wat_10kpa_05_SD <- crop(vol_wat_10kpa_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
vol_wat_10kpa_200_SD <- crop(vol_wat_10kpa_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
vol_wat_33kpa_05_SD <- crop(vol_wat_33kpa_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
vol_wat_33kpa_200_SD <- crop(vol_wat_33kpa_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
vol_wat_1500kpa_05_SD <- crop(vol_wat_1500kpa_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
vol_wat_1500kpa_200_SD <- crop(vol_wat_1500kpa_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
nitrogen_05_SD <-  crop(nitrogen_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
nitrogen_200_SD <- crop(nitrogen_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
Soil_Organic_Carbon_05_SD <- crop(Soil_Organic_Carbon_05_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))
Soil_Organic_Carbon_200_SD <- crop(Soil_Organic_Carbon_200_utm, extent(SD_box[1]-100, SD_box[3]+100, SD_box[2]-100, SD_box[4]+100))

#creating stacks of the soil raster layers (using multiple for visualization purposes)
soil_stack_SD_soil_text <- stack(clay_05_SD, clay_200_SD, silt_05_SD, silt_200_SD, sand_05_SD, sand_200_SD) #the stack of all of the soil texture rasters
soil_stack_SD_other <- stack(ph_05_SD, ph_200_SD, ocd_05_SD, ocd_200_SD, coarse_frag_05_SD, coarse_frag_200_SD, #the stack of all of the other soil variables, with different extents than the soil texture rasters
                             cat_ex_cap_05_SD, cat_ex_cap_200_SD, bulk_dens_05_SD, bulk_dens_200_SD, vol_wat_10kpa_05_SD,
                             vol_wat_10kpa_200_SD, vol_wat_33kpa_05_SD, vol_wat_33kpa_200_SD, vol_wat_1500kpa_05_SD, 
                             vol_wat_1500kpa_200_SD) 
soil_stack_SD_extra <- stack(nitrogen_05_SD, nitrogen_200_SD,Soil_Organic_Carbon_05_SD,  Soil_Organic_Carbon_200_SD)

#plotting the stacked rasters
plot(soil_stack_SD_soil_text)
plot(soil_stack_SD_soil_text, zlim = c(130, 710)) #version where the plots have the same scale
plot(soil_stack_SD_other)
plot(soil_stack_SD_other, zlim = c(45, 360)) #version where the plots have the same scale
plot(soil_stack_SD_extra)
plot(soil_stack_SD_extra, zlim = c(25, 340)) #version where the plots have the same scale

#creating X sequential columns in point data which will make it easier to select random points from each grid later

#creating an x_sequential column that is 1 through the number of LM points, which will make it easier to randomly choose one point
LM_fixed_field_data_processed <- LM_fixed_field_data_processed %>%
  mutate(X_sequential = 1:nrow(LM_fixed_field_data_processed))

#creating an x_sequential column that is 1 through the number of LC points, which will make it easier to randomly choose one point
LC_fixed_field_data_processed <- LC_fixed_field_data_processed %>%
  mutate(X_sequential = 1:nrow(LC_fixed_field_data_processed))

#creating an x_sequential column that is 1 through the number of SD points, which will make it easier to randomly choose one point
SD_fixed_field_data_processed <- SD_fixed_field_data_processed %>%
  mutate(X_sequential = 1:nrow(SD_fixed_field_data_processed))

## Extracting the soil data to the tree points 

#LM
LM_soil_text_raster_250_data_pts <- raster::extract(soil_stack_LM_soil_text, LM_fixed_field_data_processed) #extracting soil textures for each point value
LM_soil_other_raster_250_data_pts <- raster::extract(soil_stack_LM_other, LM_fixed_field_data_processed) #extracting the other soil variables for each point value
LM_soil_extra_raster_250_data_pts <- raster::extract(soil_stack_LM_extra, LM_fixed_field_data_processed) #extracting the extra soil variables for each point value
LM_fixed_field_data_processed_soils <- cbind(LM_fixed_field_data_processed, LM_soil_text_raster_250_data_pts) #bind the soil textures data for each point to the LM point dataframe
LM_fixed_field_data_processed_soils <- cbind(LM_fixed_field_data_processed_soils, LM_soil_other_raster_250_data_pts) #bind the other soil variable data for each point to the LM point dataframe
LM_fixed_field_data_processed_soils <- cbind(LM_fixed_field_data_processed_soils, LM_soil_extra_raster_250_data_pts) #bind the extra soil variable data for each point to the LM point dataframe

#LC
LC_soil_text_raster_250_data_pts <- raster::extract(soil_stack_LC_soil_text, LC_fixed_field_data_processed) #extracting soil textures for each point value
LC_soil_other_raster_250_data_pts <- raster::extract(soil_stack_LC_other, LC_fixed_field_data_processed) #extracting the other soil variables for each point value
LC_soil_extra_raster_250_data_pts <- raster::extract(soil_stack_LC_extra, LC_fixed_field_data_processed) #extracting the extra soil variables for each point value
LC_fixed_field_data_processed_soils <- cbind(LC_fixed_field_data_processed, LC_soil_text_raster_250_data_pts) #bind the soil textures data for each point to the LC point dataframe
LC_fixed_field_data_processed_soils <- cbind(LC_fixed_field_data_processed_soils, LC_soil_other_raster_250_data_pts) #bind the other soil variable data for each point to the LC point dataframe
LC_fixed_field_data_processed_soils <- cbind(LC_fixed_field_data_processed_soils, LC_soil_extra_raster_250_data_pts) #bind the extra soil variable data for each point to the LC point dataframe

#SD
SD_soil_text_raster_250_data_pts <- raster::extract(soil_stack_SD_soil_text, SD_fixed_field_data_processed) #extracting soil textures for each point value
SD_soil_other_raster_250_data_pts <- raster::extract(soil_stack_SD_other, SD_fixed_field_data_processed) #extracting the other soil variables for each point value
SD_soil_extra_raster_250_data_pts <- raster::extract(soil_stack_SD_extra, SD_fixed_field_data_processed) #extracting the extra soil variables for each point value
SD_fixed_field_data_processed_soils <- cbind(SD_fixed_field_data_processed, SD_soil_text_raster_250_data_pts) #bind the soil textures data for each point to the LC point dataframe
SD_fixed_field_data_processed_soils <- cbind(SD_fixed_field_data_processed_soils, SD_soil_other_raster_250_data_pts) #bind the other soil variable data for each point to the LC point dataframe
SD_fixed_field_data_processed_soils <- cbind(SD_fixed_field_data_processed_soils, SD_soil_extra_raster_250_data_pts) #bind the extra soil variable data for each point to the LC point dataframe


#### Creating Sandy and Clay/Loamy Available Water Columns ####

# Making four new soil metric variables for soil available water, which equals the field capacity - permanent wilting point

# LM

LM_fixed_field_data_processed_soils <- LM_fixed_field_data_processed_soils %>%
  mutate(sandy_avail_water_0.5 = vol_water_0.5 - vol_water_.1500kPa_0.5) %>% # Sand Available Water 0-5 cm
  mutate(sandy_avail_water_100.200 = vol_water_100.200 - vol_water_.1500_100.200) %>% # Sand Available Water 100-200 cm
  mutate(clay_loam_avail_water_0.5 = vol_water_.10_0.5 - vol_water_.1500kPa_0.5) %>% # Clay/Loam Available Water 0-5 cm
  mutate(clay_loam_avail_water_100.200 = vol_water_.10_100.200 - vol_water_.1500_100.200) # Clay/Loam Available Water 100-200 cm

# LC

LC_fixed_field_data_processed_soils <- LC_fixed_field_data_processed_soils %>%
  mutate(sandy_avail_water_0.5 = vol_water_0.5 - vol_water_.1500kPa_0.5) %>%
  mutate(sandy_avail_water_100.200 = vol_water_100.200 - vol_water_.1500_100.200) %>% # Sand Available Water 100-200 cm
  mutate(clay_loam_avail_water_0.5 = vol_water_.10_0.5 - vol_water_.1500kPa_0.5) %>% # Clay/Loam Available Water 0-5 cm
  mutate(clay_loam_avail_water_100.200 = vol_water_.10_100.200 - vol_water_.1500_100.200)  # Clay/Loam Available Water 100-200 cm

# SD

SD_fixed_field_data_processed_soils <- SD_fixed_field_data_processed_soils %>%
  mutate(sandy_avail_water_0.5 = vol_water_0.5 - vol_water_.1500kPa_0.5) %>% # Sand Available Water 0-5 cm
  mutate(sandy_avail_water_100.200 = vol_water_100.200 - vol_water_.1500_100.200) %>% # Sand Available Water 100-200 cm
  mutate(clay_loam_avail_water_0.5 = vol_water_.10_0.5 - vol_water_.1500kPa_0.5) %>% # Clay/Loam Available Water 0-5 cm
  mutate(clay_loam_avail_water_100.200 = vol_water_.10_100.200 - vol_water_.1500_100.200) # Clay/Loam Available Water 100-200 cm

#### Creating All Population Variables Extraction, Storing, and Uploading ####

#extracting the slope in degrees, using the queens method (neighbor = 8)
all_populations_slope_raster_15 <- terra::terrain(CEM_15_utm_all_populations, unit = 'degrees', neighbors = 8, 'slope')

#extracting the aspect in degrees, using the queens method (neighbor = 8)
all_populations_aspect_raster_15 <- terra::terrain(CEM_15_utm_all_populations, unit = 'degrees', neighbors = 8, 'aspect')


## TWI

# 1. we fill in the depressions in the DEM
dem_path <- normalizePath("./data/15m_Elevation_Raster/CEM_15_utm_all_populations.tif")
out_path <- normalizePath("./data/topographic_wetness_index", mustWork = TRUE)
wbt_fill_depressions(dem = dem_path, output = file.path(out_path, "dem_filled_all_populations.tif"), )

# 2. we want to create a slope raster using the filled DEM raster even though we already have slope rasters
wbt_slope(dem = "./data/topographic_wetness_index/dem_filled_all_populations.tif", output = "./data/topographic_wetness_index/slope_all_populations.tif")

# 3. we create the specific contributing area
# we are using the deterministic 8 which factors in the eight immediate neighbors 
wbt_d8_flow_accumulation(input = "./data/topographic_wetness_index/dem_filled_all_populations.tif", output = "./data/topographic_wetness_index/spec_catch_area_all_populations.tif", out_type = 'specific contributing area')

# 4. we create the topographic wetness index raster
wbt_wetness_index(sca = "./data/topographic_wetness_index/spec_catch_area_all_populations.tif", slope = "./data/topographic_wetness_index/slope_all_populations.tif", output = "./data/topographic_wetness_index/twi_all_populations.tif")

#importing the rasters and setting the crs

#filled DEM
dem_filled_all_populations <- raster(paste0("./data/topographic_wetness_index/dem_filled_all_populations.tif"))
terra::crs(dem_filled_all_populations) <- CRS("EPSG:26912")

#slope
slope_all_populations <- raster(paste0("./data/topographic_wetness_index/slope_all_populations.tif"))
terra::crs(slope_all_populations) <- CRS("EPSG:26912")

#specific catchment area
spec_catch_area_all_populations <- raster(paste0("./data/topographic_wetness_index/spec_catch_area_all_populations.tif"))
terra::crs(spec_catch_area_all_populations) <- CRS("EPSG:26912")

#topographic wetness index
twi_all_populations <- raster(paste0("./data/topographic_wetness_index/twi_all_populations.tif"))
terra::crs(twi_all_populations) <- CRS("EPSG:26912")


# HLI 

#need to make sure the DEM is a spatraster for the heat load index function
CEM_15_utm_all_populations <- as(CEM_15_utm_all_populations, "SpatRaster")

#creating the heat load index raster
heat.load.raster.all.pops <- hli(CEM_15_utm_all_populations, #our LM DEM 
                                 check = TRUE, #check for projection integrity and calculate central latitude for non-geographic projections
                                 force.hemisphere = "northern") #calculated based on our rasters being in the northern hemisphere


#### Generating the 20 QUBR Population Soil and Spatial Dataframe ####

#downloading the data containing the locations (lat/lon) of the 20 known populations
all_pop_locations.df <- read.csv(file = "./data/Known QUBR populations.xlsx - More accurate GPS coords for pops (12_2025).csv")

#creating a point shapefile of all points with lat lon coordinates and other attributes in WGS 1984
#sf objects are dataframes with rows representing simple features with attributes and a simple feature geometry list-column (sfc)
all_pop_locations.df_sf <- st_as_sf(all_pop_locations.df, 
                                    coords = c("Classic.Longitude", "Classic.Latitude"), crs = 4326)

#creating a transformed point shapefile with UTM 12 N an equal area projection
all_pop_locations.df_sf_transformed <- st_transform(all_pop_locations.df_sf, crs = 26912)

#create dataframe with X and Y UTM coordinates
all_pop_locations.df_sf_trans_coords <- st_coordinates(all_pop_locations.df_sf_transformed) #creates a dataframe with seperate x and y columns from the UTM 12N transformation
all_pop_locations.df_sf_trans_coordinates <- all_pop_locations.df_sf_transformed %>%
  cbind(all_pop_locations.df_sf_trans_coords) #combines the x and y coordinate data frame with the transformed sf dataframe

# importing the Baja California Sur (BCS) shapefile data

#turning the BCS polygon into a shapefile and visualize its outline
BCS_polygon <- read_sf("./data/Shapefiles/BCS_Shapefile/bcs_entidad.shp")
BCS_polygon <- st_as_sf(BCS_polygon) #making sure the BCS polygon is converted to a simple feature 
BCS_polygon_UTM <- st_transform(BCS_polygon, crs = 26912) # transforming the BCS polygon to a UTM 12 N an equal area projection
BCS_polygon_UTM <- st_as_sf(BCS_polygon_UTM) #making sure the transformed BCS polygon is converted to a simple feature 


#cropping the BCS polygon to just be the southern region of where the 20 known populations are with a 7 km radius 

# Creating a bounding box of the 20 known population points
bbox_points <- st_bbox(all_pop_locations.df_sf_trans_coordinates)

# Convert the bounding box to simple feature geometry list-column polygon
bbox_poly <- st_as_sfc(bbox_points)

# Buffering the BCS polygon by 7,000 meters (7 km)
bbox_poly_buffered <- st_buffer(bbox_poly, dist = 7000)

# Making the coordinate reference system matches of the buffered polygon matches the larger polygon
bbox_poly_buffered <- st_transform(bbox_poly_buffered, st_crs(BCS_polygon_UTM))

# Cropping (using intersection) the big polygon with buffered bounding box polygon
BCS_polygon_box_sf_cropped <- st_intersection(BCS_polygon_UTM, bbox_poly_buffered)

# Checking the cropped polygon boundaries changed (should be smaller than original bounding box of points)
print(st_bbox(BCS_polygon_box_sf_cropped))

#plotting the original extent of the BCS polygon, the buffered polygon, the cropped extent, and the 20 population points
ggplot()+
  geom_sf(data=BCS_polygon_UTM)+
  geom_sf(data=bbox_poly_buffered)+
  geom_sf(data=BCS_polygon_box_sf_cropped, color = "red")+
  geom_sf(data=all_pop_locations.df_sf_trans_coordinates)

#plotting the original extent of the cropped BCS polygon and the 20 population points
ggplot()+
  geom_sf(data=BCS_polygon_box_sf_cropped, color = "red")+
  geom_sf(data=all_pop_locations.df_sf_trans_coordinates)

#using the extent of the cropped river polygon to crop the raster for each soil texture layer
#clay content 0-5 cm
clay_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(clay_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
clay_05_clipped <- crop(clay_05_utm, clay_05_bbox_poly_buffered) # Crop the raster by the polygon
clay_05_all_pop <- mask(clay_05_clipped, clay_05_bbox_poly_buffered) # Mask the raster to the polygon
#clay content 100-200 cm
clay_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(clay_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
clay_200_clipped <- crop(clay_200_utm, clay_200_bbox_poly_buffered) # Crop the raster by the polygon
clay_200_all_pop <- mask(clay_200_clipped, clay_200_bbox_poly_buffered) # Mask the raster by the polygon
#silt content 0-5 cm
silt_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(silt_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
silt_05_clipped <- crop(silt_05_utm, silt_05_bbox_poly_buffered) # Crop the raster by the polygon
silt_05_all_pop <- mask(silt_05_clipped, silt_05_bbox_poly_buffered) # Mask the raster by the polygon
#silt content 100-200 cm
silt_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(silt_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
silt_200_clipped <- crop(silt_200_utm, silt_200_bbox_poly_buffered) # Crop the raster by the polygon
silt_200_all_pop <- mask(silt_200_clipped, silt_200_bbox_poly_buffered) # Mask the raster by the polygon
#sand content 0-5 cm
sand_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(sand_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
sand_05_clipped <- crop(sand_05_utm, sand_05_bbox_poly_buffered) # Crop the raster by the polygon
sand_05_all_pop <- mask(sand_05_clipped, sand_05_bbox_poly_buffered) # Mask the raster by the polygon
#sand content 100-200 cm
sand_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(sand_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
sand_200_clipped <- crop(sand_200_utm, sand_200_bbox_poly_buffered) # Crop the raster by the polygon
sand_200_all_pop <- mask(sand_200_clipped, sand_200_bbox_poly_buffered) # Mask the raster by the polygon
#pH 0-5 cm
ph_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(ph_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
ph_05_clipped <- crop(ph_05_utm, ph_05_bbox_poly_buffered) # Crop the raster by the polygon
ph_05_all_pop <- mask(ph_05_clipped, ph_05_bbox_poly_buffered) # Mask the raster by the polygon
#pH 100-200 cm
ph_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(ph_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
ph_200_clipped <- crop(ph_200_utm, ph_200_bbox_poly_buffered) # Crop the raster by the polygon
ph_200_all_pop <- mask(ph_200_clipped, ph_200_bbox_poly_buffered) # Mask the raster by the polygon
#Organic Carbon Content (OCD) 0-5 cm
ocd_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(ocd_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
ocd_05_clipped <- crop(ocd_05_utm, ocd_05_bbox_poly_buffered) # Crop the raster by the polygon
ocd_05_all_pop <- mask(ocd_05_clipped, ocd_05_bbox_poly_buffered) # Mask the raster by the polygon
#Organic Carbon Content (OCD) 100-200 cm
ocd_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(ocd_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
ocd_200_clipped <- crop(ocd_200_utm, ocd_200_bbox_poly_buffered) # Crop the raster by the polygon
ocd_200_all_pop <- mask(ocd_200_clipped, ocd_200_bbox_poly_buffered) # Mask the raster by the polygon
#Coarse Fragment Content 0-5 cm
coarse_frag_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(coarse_frag_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
coarse_frag_05_clipped <- crop(coarse_frag_05_utm, coarse_frag_05_bbox_poly_buffered) # Crop the raster by the polygon
coarse_frag_05_all_pop <- mask(coarse_frag_05_clipped, coarse_frag_05_bbox_poly_buffered) # Mask the raster by the polygon
#Coarse Fragment Content 100-200 cm
coarse_frag_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(coarse_frag_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
coarse_frag_200_clipped <- crop(coarse_frag_200_utm, coarse_frag_200_bbox_poly_buffered) # Crop the raster by the polygon
coarse_frag_200_all_pop <- mask(coarse_frag_200_clipped, coarse_frag_200_bbox_poly_buffered) # Mask the raster by the polygon
#Cation Exchange Capacity 0-5 cm
cat_ex_cap_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(cat_ex_cap_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
cat_ex_cap_05_clipped <- crop(cat_ex_cap_05_utm, cat_ex_cap_05_bbox_poly_buffered) # Crop the raster by the polygon
cat_ex_cap_05_all_pop <- mask(cat_ex_cap_05_clipped, cat_ex_cap_05_bbox_poly_buffered) # Mask the raster by the polygon
#Cation Exchange Capacity 100-200 cm
cat_ex_cap_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(cat_ex_cap_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
cat_ex_cap_200_clipped <- crop(cat_ex_cap_200_utm, cat_ex_cap_200_bbox_poly_buffered) # Crop the raster by the polygon
cat_ex_cap_200_all_pop <- mask(cat_ex_cap_200_clipped, cat_ex_cap_200_bbox_poly_buffered) # Mask the raster by the polygon
#Bulk Density 0-5 cm
bulk_dens_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(bulk_dens_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
bulk_dens_05_clipped <- crop(bulk_dens_05_utm, bulk_dens_05_bbox_poly_buffered) # Crop the raster by the polygon
bulk_dens_05_all_pop <- mask(bulk_dens_05_clipped, bulk_dens_05_bbox_poly_buffered) # Mask the raster by the polygon
#Bulk Density 100-200 cm
bulk_dens_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(bulk_dens_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
bulk_dens_200_clipped <- crop(bulk_dens_200_utm, bulk_dens_200_bbox_poly_buffered) # Crop the raster by the polygon
bulk_dens_200_all_pop <- mask(bulk_dens_200_clipped, bulk_dens_200_bbox_poly_buffered) # Mask the raster by the polygon
#Volume of Water Content at -10 kpa 0-5 cm
vol_wat_10kpa_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(vol_wat_10kpa_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
vol_wat_10kpa_05_clipped <- crop(vol_wat_10kpa_05_utm, vol_wat_10kpa_05_bbox_poly_buffered) # Crop the raster by the polygon
vol_wat_10kpa_05_all_pop <- mask(vol_wat_10kpa_05_clipped, vol_wat_10kpa_05_bbox_poly_buffered) # Mask the raster by the polygon
#Volume of Water Content at -10 kpa 100-200 cm
vol_wat_10kpa_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(vol_wat_10kpa_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
vol_wat_10kpa_200_clipped <- crop(vol_wat_10kpa_200_utm, vol_wat_10kpa_200_bbox_poly_buffered) # Crop the raster by the polygon
vol_wat_10kpa_200_all_pop <- mask(vol_wat_10kpa_200_clipped, vol_wat_10kpa_200_bbox_poly_buffered) # Mask the raster by the polygon
#Volume of Water Content at -33 kpa 0-5 cm
vol_wat_33kpa_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(vol_wat_33kpa_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
vol_wat_33kpa_05_clipped <- crop(vol_wat_33kpa_05_utm, vol_wat_33kpa_05_bbox_poly_buffered) # Crop the raster by the polygon
vol_wat_33kpa_05_all_pop <- mask(vol_wat_33kpa_05_clipped, vol_wat_33kpa_05_bbox_poly_buffered) # Mask the raster by the polygon
#Volume of Water Content at -33 kpa 100-200 cm
vol_wat_33kpa_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(vol_wat_33kpa_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
vol_wat_33kpa_200_clipped <- crop(vol_wat_33kpa_200_utm, vol_wat_33kpa_200_bbox_poly_buffered) # Crop the raster by the polygon
vol_wat_33kpa_200_all_pop <- mask(vol_wat_33kpa_200_clipped, vol_wat_33kpa_200_bbox_poly_buffered) # Mask the raster by the polygon
#Volume of Water Content at -1500 kpa 0-5 cm
vol_wat_1500kpa_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(vol_wat_1500kpa_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
vol_wat_1500kpa_05_clipped <- crop(vol_wat_1500kpa_05_utm, vol_wat_1500kpa_05_bbox_poly_buffered) # Crop the raster by the polygon
vol_wat_1500kpa_05_all_pop <- mask(vol_wat_1500kpa_05_clipped, vol_wat_1500kpa_05_bbox_poly_buffered) # Mask the raster by the polygon
#Volume of Water Content at -1500 kpa 100-200 cm
vol_wat_1500kpa_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(vol_wat_1500kpa_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
vol_wat_1500kpa_200_clipped <- crop(vol_wat_1500kpa_200_utm, vol_wat_1500kpa_200_bbox_poly_buffered) # Crop the raster by the polygon
vol_wat_1500kpa_200_all_pop <- mask(vol_wat_1500kpa_200_clipped, vol_wat_1500kpa_200_bbox_poly_buffered) # Mask the raster by the polygon
#Nitrogen Content 0-5 cm
nitrogen_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(nitrogen_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
nitrogen_05_clipped <- crop(nitrogen_05_utm, nitrogen_05_bbox_poly_buffered) # Crop the raster by the polygon
nitrogen_05_all_pop <- mask(nitrogen_05_clipped, nitrogen_05_bbox_poly_buffered) # Mask the raster by the polygon
#Nitrogen Content 100-200 cm
nitrogen_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(nitrogen_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
nitrogen_200_clipped <- crop(nitrogen_200_utm, nitrogen_200_bbox_poly_buffered) # Crop the raster by the polygon
nitrogen_200_all_pop <- mask(nitrogen_200_clipped, nitrogen_200_bbox_poly_buffered) # Mask the raster by the polygon
#Soil Organic Carbon 0-5 cm
Soil_Organic_Carbon_05_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(Soil_Organic_Carbon_05_utm)) # Make sure your sfc polygon is in the same CRS as the raster
Soil_Organic_Carbon_05_clipped <- crop(Soil_Organic_Carbon_05_utm, Soil_Organic_Carbon_05_bbox_poly_buffered) # Crop the raster by the polygon
Soil_Organic_Carbon_05_all_pop <- mask(Soil_Organic_Carbon_05_clipped, Soil_Organic_Carbon_05_bbox_poly_buffered) # Mask the raster by the polygon
#Soil Organic Carbon 100-200 cm
Soil_Organic_Carbon_200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(Soil_Organic_Carbon_200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
Soil_Organic_Carbon_200_clipped <- crop(Soil_Organic_Carbon_200_utm, Soil_Organic_Carbon_200_bbox_poly_buffered) # Crop the raster by the polygon
Soil_Organic_Carbon_200_all_pop <- mask(Soil_Organic_Carbon_200_clipped, Soil_Organic_Carbon_200_bbox_poly_buffered) # Mask the raster by the polygon
#Sand Available Water Content 0-5 cm
sandy_avail_water_0.5_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(sandy_avail_water_0.5_utm)) # Make sure your sfc polygon is in the same CRS as the raster
sandy_avail_water_0.5_clipped <- crop(sandy_avail_water_0.5_utm, sandy_avail_water_0.5_bbox_poly_buffered) # Crop the raster by the polygon
sandy_avail_water_0.5_all_pop <- mask(sandy_avail_water_0.5_clipped, sandy_avail_water_0.5_bbox_poly_buffered) # Mask the raster by the polygon
#Sand Available Water Content 100-200 cm
sandy_avail_water_100.200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(sandy_avail_water_100.200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
sandy_avail_water_100.200_clipped <- crop(sandy_avail_water_100.200_utm, sandy_avail_water_100.200_bbox_poly_buffered) # Crop the raster by the polygon
sandy_avail_water_100.200_all_pop <- mask(sandy_avail_water_100.200_clipped, sandy_avail_water_100.200_bbox_poly_buffered) # Mask the raster by the polygon
# Clay/Loam Available Water Content 0-5 cm
clay_loam_avail_water_0.5_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(clay_loam_avail_water_0.5_utm)) # Make sure your sfc polygon is in the same CRS as the raster
clay_loam_avail_water_0.5_clipped <- crop(clay_loam_avail_water_0.5_utm, clay_loam_avail_water_0.5_bbox_poly_buffered) # Crop the raster by the polygon
clay_loam_avail_water_0.5_all_pop <- mask(clay_loam_avail_water_0.5_clipped, clay_loam_avail_water_0.5_bbox_poly_buffered) # Mask the raster by the polygon
# Clay/Loam Available Water Content 100-200 cm
clay_loam_avail_water_100.200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(clay_loam_avail_water_100.200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
clay_loam_avail_water_100.200_clipped <- crop(clay_loam_avail_water_100.200_utm, clay_loam_avail_water_100.200_bbox_poly_buffered) # Crop the raster by the polygon
clay_loam_avail_water_100.200_all_pop <- mask(clay_loam_avail_water_100.200_clipped, clay_loam_avail_water_100.200_bbox_poly_buffered) # Mask the raster by the polygon
# Elevation
clay_loam_avail_water_100.200_bbox_poly_buffered <- st_transform(BCS_polygon_box_sf_cropped, crs(clay_loam_avail_water_100.200_utm)) # Make sure your sfc polygon is in the same CRS as the raster
clay_loam_avail_water_100.200_clipped <- crop(clay_loam_avail_water_100.200_utm, clay_loam_avail_water_100.200_bbox_poly_buffered) # Crop the raster by the polygon
clay_loam_avail_water_100.200_all_pop <- mask(clay_loam_avail_water_100.200_clipped, clay_loam_avail_water_100.200_bbox_poly_buffered) # Mask the raster by the polygon



#Plotting a confirmation that we properly cropped the rasters by plotting the clay rasters (blue) with the cropped polygon (green) around it
ggplot()+
  geom_raster(data = as.data.frame(clay_05_utm, xy=T), aes(x=x, y=y, fill = clay.content.0.5)) +
  scale_fill_gradientn(colours=c("yellow","red"), name = "clay_05_utm")+
  # Add new fill scale
  ggnewscale::new_scale_fill() +
  geom_raster(data = as.data.frame(clay_05_all_pop, xy=T), aes(x=x, y=y, fill = clay.content.0.5)) +
  scale_fill_gradientn(colours=c("lightblue","darkblue"), name = "clay_05_all_pop") +
  geom_sf(data = BCS_polygon_box_sf_cropped, fill = NA, color = "green")+
  geom_sf(data=all_pop_locations.df_sf_trans_coordinates)

#creating raster stacks of the 0-5 and 100-200 cm raster layers for each soil metric for later cell data extraction
soil_stack_clay <- stack(clay_05_all_pop, clay_200_all_pop) #stacked clay
soil_stack_silt <- stack(silt_05_all_pop, silt_200_all_pop) #stacked silt
soil_stack_sand <- stack(sand_05_all_pop, sand_200_all_pop) #stacked sand
soil_stack_ph <- stack(ph_05_all_pop, ph_200_all_pop) #stacked ph
soil_stack_ocd <- stack(ocd_05_all_pop, ocd_200_all_pop)  #stacked ocd
soil_stack_coarse_frag <- stack(coarse_frag_05_all_pop, coarse_frag_200_all_pop)  #stacked coarse fragment
soil_stack_cat_ex <- stack(cat_ex_cap_05_all_pop, cat_ex_cap_200_all_pop)  #stacked cation exchange capacity
soil_stack_bulk_dens <- stack(bulk_dens_05_all_pop, bulk_dens_200_all_pop)  #stacked bulk density
soil_stack_vol_wat_10kpa <- stack(vol_wat_10kpa_05_all_pop, vol_wat_10kpa_200_all_pop)  #stacked volume water content at -10 kpa
soil_stack_vol_wat_33kpa <- stack(vol_wat_33kpa_05_all_pop, vol_wat_33kpa_200_all_pop)  #stacked volume water content at -33 kpa
soil_stack_vol_wat_1500kpa <- stack(vol_wat_1500kpa_05_all_pop, vol_wat_1500kpa_200_all_pop) #stacked volume water content at -1500 kpa
soil_stack_nitrogen <- stack(nitrogen_05_all_pop, nitrogen_200_all_pop)  #stacked nitrogen
soil_stack_soc <- stack(Soil_Organic_Carbon_05_all_pop, Soil_Organic_Carbon_200_all_pop)  #stacked soil organic carbon
soil_stack_sandy_water <- stack(sandy_avail_water_0.5_all_pop, sandy_avail_water_100.200_all_pop) #stacked sandy available water
soil_stack_clay_loam_water <- stack(clay_loam_avail_water_0.5_all_pop, clay_loam_avail_water_100.200_all_pop) #stacked clay/loam available water

#plotting a stacked raster, example with clay content at 0-5 and 100-200 cm
plot(soil_stack_clay) #version with soil textures
plot(soil_stack_clay, zlim = c(0, 350)) #version where the plots have the same scale

#extracting and storing the soil data for each known population point for each soil metric
all_known_pop_soil_clay <- raster::extract(soil_stack_clay, all_pop_locations.df_sf_trans_coordinates) #extracting clay
all_known_pop_soil_silt <- raster::extract(soil_stack_silt, all_pop_locations.df_sf_trans_coordinates) #extracting silt
all_known_pop_soil_sand <- raster::extract(soil_stack_sand, all_pop_locations.df_sf_trans_coordinates) #extracting sand
all_known_pop_soil_ph <- raster::extract(soil_stack_ph, all_pop_locations.df_sf_trans_coordinates) #extracting ph
all_known_pop_soil_ocd <- raster::extract(soil_stack_ocd, all_pop_locations.df_sf_trans_coordinates) #extracting ocd
all_known_pop_soil_coarse_frag <- raster::extract(soil_stack_coarse_frag, all_pop_locations.df_sf_trans_coordinates) #extracting coarse fragment
all_known_pop_soil_cat_ex <- raster::extract(soil_stack_cat_ex, all_pop_locations.df_sf_trans_coordinates) #extracting cation exchange capacity
all_known_pop_soil_bulk_dens <- raster::extract(soil_stack_bulk_dens, all_pop_locations.df_sf_trans_coordinates) #extracting bulk density
all_known_pop_soil_vol_wat_10kpa <- raster::extract(soil_stack_vol_wat_10kpa, all_pop_locations.df_sf_trans_coordinates) # extracting volume water content at -10 kpa
all_known_pop_soil_vol_wat_33kpa <- raster::extract(soil_stack_vol_wat_33kpa, all_pop_locations.df_sf_trans_coordinates) # extracting volume water content at -33 kpa
all_known_pop_soil_vol_wat_1500kpa <- raster::extract(soil_stack_vol_wat_1500kpa, all_pop_locations.df_sf_trans_coordinates) # extracting volume water content at -1500 kpa
all_known_pop_soil_nitrogen <- raster::extract(soil_stack_nitrogen, all_pop_locations.df_sf_trans_coordinates) # extracting nitrogen
all_known_pop_soil_soc <- raster::extract(soil_stack_soc, all_pop_locations.df_sf_trans_coordinates) # extracting soil organic carbon
all_known_pop_soil_sandy_water <- raster::extract(soil_stack_sandy_water, all_pop_locations.df_sf_trans_coordinates) # extracting
all_known_pop_soil_clay_loam_water <- raster::extract(soil_stack_clay_loam_water, all_pop_locations.df_sf_trans_coordinates) # extracting
all_known_pop_soil_elevation <- raster::extract(CEM_15_utm_all_populations, all_pop_locations.df_sf_trans_coordinates) # extracting
all_known_pop_soil_slope <- raster::extract(all_populations_slope_raster_15, all_pop_locations.df_sf_trans_coordinates) # extracting
all_known_pop_soil_aspect <- raster::extract(all_populations_aspect_raster_15, all_pop_locations.df_sf_trans_coordinates) # extracting
all_known_pop_soil_TWI <- raster::extract(twi_all_populations, all_pop_locations.df_sf_trans_coordinates) # extracting
all_known_pop_soil_HLI <- raster::extract(heat.load.raster.all.pops, all_pop_locations.df_sf_trans_coordinates) # extracting


# All Points

#creating the eastness and northness values
all_populations_aspect_raster_15
#first, converting aspect values to radian
all_known_pop_soil_aspect_radian = ((all_known_pop_soil_aspect * pi) / 180) # creating a column that is the radians

#creating eastness
all_known_pop_soil_eastness = sin(all_known_pop_soil_aspect_radian) # creating the eastness column

#creating northness
all_known_pop_soil_northness = cos(all_known_pop_soil_aspect_radian) # creating the northness column



#combine all of the extracted population soil textures data across all soil metrics
all_known_pop_soils <- cbind(all_pop_locations.df_sf_trans_coordinates, all_known_pop_soil_clay) #adding clay to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_silt) #adding silt to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_sand) #adding sand to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_ph) #adding ph to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_ocd) #adding organic carbon density content to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_cat_ex) #adding cation exchange capacity to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_bulk_dens) #adding bulk density to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_vol_wat_10kpa) #adding volume water content at -10 kpa to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_vol_wat_33kpa) #adding volume water content at -33 kpa to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_vol_wat_1500kpa) #adding volume water content at -1500 kpa to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_nitrogen) #adding nitrogen to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_soc) #adding soil organic carbon to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_sandy_water) #adding sandy available water to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_clay_loam_water) #adding clay/loam available waterr to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_elevation) #adding elevation to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_slope) #adding slope to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_eastness) #adding eastness to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_northness) #adding northness to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_aspect) #adding aspect to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_TWI) #adding topographic wetness index to the dataframe
all_known_pop_soils <- cbind(all_known_pop_soils, all_known_pop_soil_HLI) #adding heat load index to the dataframe

#correcting column names for clarity
all_known_pop_soils <- all_known_pop_soils %>%
  mutate(sandy_avail_water_0.5 = layer.1) %>%
  mutate(sandy_avail_water_100.200 = layer.2) %>% 
  mutate(clay_loam_avail_water_0.5 = layer.1.1) %>%
  mutate(clay_loam_avail_water_100.200 = layer.2.1)  %>%
  mutate(all_known_pop_soil_HLI = lyr.1)


## Creating the distance to coast column ##

#plotting the original extent of the BCS polygon, the buffered polygon, the cropped extent, and the 20 population points
ggplot()+
  geom_sf(data=BCS_polygon_UTM)+
  geom_sf(data=bbox_poly_buffered)+
  geom_sf(data=BCS_polygon_box_sf_cropped, color = "red")+
  geom_sf(data=all_pop_locations.df_sf_trans_coordinates)

#  converting the all populations locations to spatial data with longitudes and latitudes to be able to calculate distances
all_pop_locations.df_sf_trans_coordinates_sp <- st_transform(all_pop_locations.df_sf_trans_coordinates, crs = 4326)
all_pop_locations.df_sf_trans_coordinates_sp <- as(st_geometry(all_pop_locations.df_sf_trans_coordinates_sp$geometry), "Spatial")

#  converting the BCS polygon to spatial data with longitudes and latitudes to be able to calculate distances
BCS_polygon_UTM_sp_coast <- st_transform(BCS_polygon_UTM, crs = 4326)
BCS_polygon_UTM_sp_coast <- as(st_geometry(BCS_polygon_UTM_sp_coast), "Spatial")

# calculates the shortest distance (meters) of each population to the coast
Distance <- dist2Line(p = all_pop_locations.df_sf_trans_coordinates_sp, 
                      line = BCS_polygon_UTM_sp_coast)

# add the distance to coast column 
all_known_pop_soils <- cbind(all_known_pop_soils, Distance)