Skip to content
This repository has been archived by the owner on Oct 21, 2021. It is now read-only.

Commit

Permalink
Added pdsi timeseries
Browse files Browse the repository at this point in the history 
Estimated soil water with thornthwaite and calculated PDSI from timeseries incorporating field AWC
  • Loading branch information
@flakesw
flakesw committed Dec 23, 2020
1 parent b668fac commit 30fd193
Show file tree
Hide file tree
Showing 7 changed files with 166,211 additions and 13,029 deletions.
730 changes: 365 additions & 365 deletions .Rhistory
View file

Large diffs are not rendered by default.

2 changes: 2 additions & 0 deletions .Rproj.user/shared/notebooks/paths
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
C:/Users/Sam/Documents/Research/MS Thesis/Understory/Raw data/ALL_climate_variables.csv="FA323732"
C:/Users/Sam/Documents/Research/MS Thesis/Understory/Raw data/soils_missing_imputed_with_new_ptf_2020-12-7.csv="4DD99E8C"
C:/Users/Sam/Documents/Research/MS Thesis/Understory/plot_level_understory.R="3F501BF8"
C:/Users/Sam/Documents/Research/MS Thesis/Understory/plot_level_understory_tc.R="B649CF02"
C:/Users/Sam/Documents/Research/MS Thesis/Understory/prism_data_extract_timeseries.R="F70A5D84"
Expand Down
103 changes: 103 additions & 0 deletions Raw data/soil_new_awc.csv
View file

Large diffs are not rendered by default.

165,650 changes: 153,001 additions & 12,649 deletions clean data/climate_data_monthly.csv
View file

Large diffs are not rendered by default.

12,649 changes: 12,649 additions & 0 deletions clean data/climate_data_yearly.csv
View file

Large diffs are not rendered by default.

75 changes: 65 additions & 10 deletions prism_data_extract_timeseries.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,12 +104,6 @@ clim_data[clim_data$vpd < 0, "vpd"]<- 0

write.csv(clim_data, file="./clean data/climate_data_monthly.csv")
clim_data <- read.csv("./clean data/climate_data_monthly.csv")
#some code to plot points on top of raster to make sure stuff looks right
# e <- drawExtent()
# plot(e)
# plot(tdrast, add = TRUE)
# points(sites)


#initialize data frame
clim_ann <- data.frame(plot = rep(unique(clim_data$site), each = nyear),
Expand All @@ -125,7 +119,6 @@ clim_ann <- data.frame(plot = rep(unique(clim_data$site), each = nyear),
ppt_tot = NA)



#extract winter precip and last-fall + current-season VPD
for (i in 1:npoints){
plot <- clim_ann$plot[i * nyear]
Expand Down Expand Up @@ -190,10 +183,72 @@ for (i in 1:nrow(clim_ann)){
clim_ann$fdsi <- .44*clim_ann$stdP - .56*clim_ann$stdVPD


write.csv(clim_ann, file = "./clean data/climate_data_monthly.csv")

#-------------------------------------------------------------------------------
# Thornthwaite water balance
#-------------------------------------------------------------------------------

library("ClimClass")
library("scPDSI")
library("dplyr")
library("purrr")
library("tidyr")
library("tibble")
library("TSstudio")

clim_month <- read.csv("./clean data/climate_data_monthly.csv")
clim_normals <- clim_month[between(clim_month$year, 1980, 2010), ] %>%
stats::aggregate(by = list(.$site), FUN = mean)
soil <- read.csv("./Raw data/soil_new_awc.csv")


latitudes <- as.data.frame(sites@coords)$YCoord

months <- as.factor(month.abb)

tt_vars <- clim_month %>%
dplyr::select(site, year, month, ppt, tmean) %>%
rename(P = ppt, Tm = tmean)%>%
split(.$site) %>%
{purrr::pmap(list(series = ., latitude = latitudes, TAW = soil$mean_fc),
~thornthwaite(series = ..1, latitude = ..2, TAW = ..3))}


test_site <- "DES1014"

ppt <- tt_vars %>%
{purrr::map(., ~pluck(..1, "W_balance", "Precipitation"))} %>%
{purrr::map(., ~rownames_to_column(., var = "month"))} %>%
{purrr::map(., ~pivot_longer(..1, cols = -month,
names_to = "year",
values_to = "ppt"))} %>%
{purrr::map(., ~mutate(..1, month = factor(month, levels = months)))} %>%
{purrr::map(., ~arrange(..1, year, month))}

pet <- tt_vars %>%
{purrr::map(., ~pluck(..1, "W_balance", "Et0"))} %>%
{purrr::map(., ~rownames_to_column(., var = "month"))} %>%
{purrr::map(., ~pivot_longer(..1, cols = -month,
names_to = "year",
values_to = "Et0"))} %>%
{purrr::map(., ~mutate(..1, month = factor(month, levels = months)))} %>%
{purrr::map(., ~arrange(..1, year, month))}

pdsi <- pmap(list(P = ppt, PE = pet),
~pdsi(P = ..1$ppt, PE = ..2$Et0)) %>%
map(., ~pluck(..1, "X")) %>%
map(., ~ts_reshape(..1)) %>%
map(., ~pivot_longer(..1, cols = -month, names_to = "year")) %>%
map(., ~arrange(..1, year, month)) %>%
bind_rows(. , .id = "column_label") %>%
mutate(year = as.numeric(year)) %>%
mutate(year = year + 1894)


pdsi_ann <- stats::aggregate(pdsi$value,
by = list(pdsi$column_label, pdsi$year),
FUN = mean) %>%
`colnames<-`(c("plot", "year", "pdsi"))

clim_ann <- merge(clim_ann, pdsi_ann, by = c("plot", "year"))

write.csv(clim_ann, file = "./clean data/climate_data_yearly.csv")

31 changes: 26 additions & 5 deletions soils.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,8 @@ vangenuchten <- function(qr, qs, a, n, m, h){
return(1/((1 + (a * h)^n)^m))
}

soil2 <- soil

soil2$awc <- (soil2$Field_cap - soil2$Perm_wilt)/100
soil2$awc1 <- soil2$X330 - soil2$X15000
soil2$awc2 <- soil2$X330.1 - soil2$X15000.1
Expand All @@ -115,28 +117,39 @@ for(i in 1:102){
brookscorey(ptf_params[i+1, 10], ptf_params[i+1, 11],
ptf_params[i+1, 12], ptf_params[i+1, 13],
15000)

soil2$X330.4[i] <- brookscorey(ptf_params[i+1, 10], ptf_params[i+1, 11],
ptf_params[i+1, 12], ptf_params[i+1, 13],
330)
#campbell shiosawa
soil2$awc6[i] <- brookscorey(ptf_params[i+1, 15], ptf_params[i+1, 16],
ptf_params[i+1, 17], ptf_params[i+1, 18],
330) -
brookscorey(ptf_params[i+1, 15], ptf_params[i+1, 16],
ptf_params[i+1, 17], ptf_params[i+1, 18],
15000)
soil2$x330.5[i] <- brookscorey(ptf_params[i+1, 15], ptf_params[i+1, 16],
ptf_params[i+1, 17], ptf_params[i+1, 18],
330)
#rawls brackensiek
soil2$awc7[i] <- brookscorey(ptf_params[i+1, 20], ptf_params[i+1, 21],
ptf_params[i+1, 22], ptf_params[i+1, 23],
330) -
brookscorey(ptf_params[i+1, 20], ptf_params[i+1, 21],
ptf_params[i+1, 22], ptf_params[i+1, 23],
15000)
soil2$x330.6[i] <- brookscorey(ptf_params[i+1, 20], ptf_params[i+1, 21],
ptf_params[i+1, 22], ptf_params[i+1, 23],
330)
#williams
soil2$awc8[i] <- brookscorey(ptf_params[i+1, 25], ptf_params[i+1, 26],
ptf_params[i+1, 27], ptf_params[i+1, 28],
330) -
brookscorey(ptf_params[i+1, 25], ptf_params[i+1, 26],
ptf_params[i+1, 27], ptf_params[i+1, 28],
15000)
soil2$x330.7[i] <- brookscorey(ptf_params[i+1, 25], ptf_params[i+1, 26],
ptf_params[i+1, 27], ptf_params[i+1, 28],
330)

#oosterveld_chang
soil2$awc9[i] <- brookscorey(ptf_params[i+1, 35], ptf_params[i+1, 36],
Expand All @@ -145,13 +158,19 @@ for(i in 1:102){
brookscorey(ptf_params[i+1, 35], ptf_params[i+1, 36],
ptf_params[i+1, 37], ptf_params[i+1, 38],
15000)
soil2$x330.8[i] <- brookscorey(ptf_params[i+1, 35], ptf_params[i+1, 36],
ptf_params[i+1, 37], ptf_params[i+1, 38],
330)
#wosten_chang
soil2$awc10[i] <- vangenuchten(ptf_params[i+1, 45], ptf_params[i+1, 46],
ptf_params[i+1, 47], ptf_params[i+1, 48],
ptf_params[i+1, 49], 330) -
vangenuchten(ptf_params[i+1, 45], ptf_params[i+1, 46],
ptf_params[i+1, 47], ptf_params[i+1, 48],
ptf_params[i+1, 49], 15000)
soil2$x330.9[i] <- vangenuchten(ptf_params[i+1, 45], ptf_params[i+1, 46],
ptf_params[i+1, 47], ptf_params[i+1, 48],
ptf_params[i+1, 49], 330)

#varallyay
soil2$awc11[i] <- vangenuchten(ptf_params[i+1, 51], ptf_params[i+1, 52],
Expand All @@ -160,10 +179,12 @@ for(i in 1:102){
vangenuchten(ptf_params[i+1, 51], ptf_params[i+1, 52],
ptf_params[i+1, 53], ptf_params[i+1, 54],
ptf_params[i+1, 55], 15000)

soil2$x330.10[i] <- vangenuchten(ptf_params[i+1, 51], ptf_params[i+1, 52],
ptf_params[i+1, 53], ptf_params[i+1, 54],
ptf_params[i+1, 55], 330)
}

soil2$mean_awc <- apply(soil2[, c(27:38)], 1, mean)

soil2$mean_awc <- apply(soil2[, grep("awc", names(soil2))], 1, mean)
soil2$mean_fc <- apply(soil2[, grep("330", names(soil2))], 1, mean)

write.csv(soil2, "./Raw data/soil_")
write.csv(soil2, "./Raw data/soil_new_awc.csv")

0 comments on commit 30fd193

Please sign in to comment.