Skip to content

Commit

Permalink
Fixed bug in aggregation 'dailyPrecipitationFreeEventDistribution'
Browse files Browse the repository at this point in the history 
- old version of aggregation ‘dailyPrecipitationFreeEventDistribution’
only increased index ‘nv’ of result object ‘resMeans’ if there was at
least one precipitation free spell; if it rained on every day, then
‘nv’ wasn’t increased appropriately, and following aggregation results
were written to incorrect positions
- old version of aggregation ‘dailyPrecipitationEventSizeDistribution’
failed if there was no precipitation at all

- aggregations ‘dailyPrecipitationFreeEventDistribution’ and
‘dailyPrecipitationEventSizeDistribution’ performed similar tasks, but
used independent code; this code is now unified in a new function
‘tabulate_values_in_bins’
  • Loading branch information
dschlaep committed Sep 14, 2016
1 parent ba4c169 commit 2829718
Show file tree
Hide file tree
Showing 2 changed files with 63 additions and 48 deletions.
71 changes: 23 additions & 48 deletions 2_SWSF_p4of5_Code_v51.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3144,30 +3144,19 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer
if(print.debug) print("Aggregation of dailyPrecipitationEventSizeDistribution")
if(!exists("prcp.dy")) prcp.dy <- get_PPT_dy(sc, runData, simTime)

#prcp-event sizes in bins
events <- lapply(simTime$useyrs, FUN=function(y) floor((temp <- prcp.dy$ppt[simTime2$year_ForEachUsedDay == y])[temp>0]/bin.prcpSizes)*bin.prcpSizes)
bins.summary <- (0:6) * bin.prcpSizes #aggregate to maximal 7 bins
bins.available <- sort(unique(unlist(events))) #bins present
counts.available <- as.matrix(sapply(simTime$useyrs - startyr + 1, FUN=function(y) sapply(bins.available, FUN=function(b) sum(events[[y]] == b))))
counts.summary <- matrix(data=0, ncol=simTime$no.useyr, nrow=length(bins.summary))
counts.summary[bins.summary %in% bins.available, ] <- counts.available[bins.available %in% bins.summary, ]
if(max(bins.summary) < max(bins.available)){ #if bins available that are outside the summary bins, then sum them up into the highest summary bin
if(length(bins.suming <- which(bins.available >= max(bins.summary))) > 1){
counts.summary[length(bins.summary), ] <- apply(counts.available[bins.suming, ], MARGIN=2, FUN=sum)
} else {
counts.summary[length(bins.summary), ] <- counts.available[bins.suming, ]
}
}
eventsPerYear <- apply(counts.summary, MARGIN=2, FUN=sum)
freq.summary <- sweep(counts.summary, MARGIN=2, STATS=eventsPerYear, FUN="/")

resMeans[nv] <- mean(eventsPerYear)
resSDs[nv] <- sd(eventsPerYear)
resMeans[(nv+1):(nv+7)] <- apply(freq.summary, MARGIN=1, FUN=mean)
resSDs[(nv+1):(nv+7)] <- apply(freq.summary, MARGIN=1, FUN=sd)
#prcp-event sizes in bins
ppt_sizes <- tabulate_values_in_bins(
x = prcp.dy$ppt, method = "values", vcrit = 0,
bins = bin.prcpSizes, nbins = 7,
simTime = simTime, simTime2 = simTime2)

resMeans[nv] <- mean(ppt_sizes$eventsPerYear)
resSDs[nv] <- sd(ppt_sizes$eventsPerYear)
resMeans[(nv+1):(nv+7)] <- apply(ppt_sizes$freq.summary, 1, mean)
resSDs[(nv+1):(nv+7)] <- apply(ppt_sizes$freq.summary, 1, sd)
nv <- nv+8

rm(events, counts.available, counts.summary, freq.summary, eventsPerYear)
rm(ppt_sizes)
}
#15
if(any(simulation_timescales=="yearly") & aon$yearlyPET){
Expand Down Expand Up @@ -3224,8 +3213,8 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer
monthlyTemp = meanmonthlyTemp)

resMeans[nv:(nv+2)] <- unlist(di.normals)
resMeans[(nv+3):(nv+5)] <- apply(temp <- cbind(di.ts$ai, di.ts$TD, di.ts$criteria12), MARGIN=2, FUN=mean, na.rm=TRUE)
resSDs[(nv+3):(nv+5)] <- apply(temp, MARGIN=2, FUN=sd, na.rm=TRUE)
resMeans[(nv+3):(nv+5)] <- sapply(di.ts, mean, na.rm = TRUE)
resSDs[(nv+3):(nv+5)] <- sapply(di.ts, sd, na.rm = TRUE)
nv <- nv+6

rm(di.ts, di.normals)
Expand Down Expand Up @@ -3273,33 +3262,19 @@ do_OneSite <- function(i_sim, i_labels, i_SWRunInformation, i_sw_input_soillayer
if(print.debug) print("Aggregation of dailyPrecipitationFreeEventDistribution")
if(!exists("prcp.dy")) prcp.dy <- get_PPT_dy(sc, runData, simTime)

#duration of prcp-free days in bins
durations <- lapply(simTime$useyrs, FUN=function(y) floor(((temp <- rle(prcp.dy$ppt[simTime2$year_ForEachUsedDay == y] == 0))$lengths[temp$values]-1)/bin.prcpfreeDurations)*bin.prcpfreeDurations )
if(length(unlist(durations))==0) {if(!be.quiet) print("Unable to complete aggregation of dailyPrecipitationFreeEventDistribution")
}else{
bins.summary <- (0:3) * bin.prcpfreeDurations #aggregate to maximal 4 bins
bins.available <- sort(unique(unlist(durations))) #bins present
counts.available <- as.matrix(sapply(simTime$useyrs - startyr + 1, FUN=function(y) sapply(bins.available, FUN=function(b) sum(durations[[y]] == b))))
counts.summary <- matrix(data=0, ncol=simTime$no.useyr, nrow=length(bins.summary))
counts.summary[bins.summary %in% bins.available, ] <- counts.available[bins.available %in% bins.summary, ]
if(max(bins.summary) < max(bins.available)){ #if bins available that are outside the summary bins, then sum them up into the highest summary bin
if(length(bins.suming <- which(bins.available >= max(bins.summary))) > 1){
counts.summary[length(bins.summary), ] <- apply(counts.available[bins.suming, ], MARGIN=2, FUN=sum)
} else {
counts.summary[length(bins.summary), ] <- counts.available[bins.suming, ]
}
}
eventsPerYear <- apply(counts.summary, MARGIN=2, FUN=sum)
freq.summary <- sweep(counts.summary, MARGIN=2, STATS=eventsPerYear, FUN="/")
#duration of prcp-free days in bins
ppt_free <- tabulate_values_in_bins(
x = prcp.dy$ppt <= tol, method = "duration",
bins = bin.prcpfreeDurations, nbins = 4,
simTime = simTime, simTime2 = simTime2)

resMeans[nv] <- mean(eventsPerYear)
resSDs[nv] <- sd(eventsPerYear)
resMeans[(nv+1):(nv+4)] <- apply(freq.summary, MARGIN=1, FUN=mean)
resSDs[(nv+1):(nv+4)] <- apply(freq.summary, MARGIN=1, FUN=sd)
resMeans[nv] <- mean(ppt_free$eventsPerYear)
resSDs[nv] <- sd(ppt_free$eventsPerYear)
resMeans[(nv+1):(nv+4)] <- apply(ppt_free$freq.summary, 1, mean)
resSDs[(nv+1):(nv+4)] <- apply(ppt_free$freq.summary, 1, sd)
nv <- nv+5

rm(durations, counts.available, counts.summary, freq.summary, eventsPerYear)
}
rm(ppt_free)
}
#21
if(any(simulation_timescales=="monthly") & aon$monthlySPEIEvents){
Expand Down
40 changes: 40 additions & 0 deletions 2_SWSF_p5of5_Functions_v51.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1787,6 +1787,46 @@ daily_spells_permonth <- compiler::cmpfun(function(x, simTime2) {
matrix(temp, nrow = 12)
})

tabulate_values_in_bins <- compiler::cmpfun(function(x, method = c("duration", "values"),
vcrit = NULL, bins, nbins, simTime, simTime2) {
method <- match.arg(method)

bins.summary <- (seq_len(nbins) - 1) * bins

dat <- if (method == "duration" && is.logical(x)) {
# bin duration of TRUE runs
lapply(simTime$useyrs, function(y) {
temp <- rle(x[simTime2$year_ForEachUsedDay == y])
temp <- floor((temp$lengths[temp$values] - 1) / bins) * bins
findInterval(temp, vec = bins.summary)
})
} else if (method == "values") {
# bin values
lapply(simTime$useyrs, function(y) {
temp <- x[simTime2$year_ForEachUsedDay == y]
if (!is.null(vcrit)) temp <- temp[temp > vcrit]
floor(temp / bins) * bins
findInterval(temp, vec = bins.summary)
})
} else {
print("'tabulate_values_in_bins' cannot be calculated")
NULL
}

if (length(unlist(dat)) > 0) {
counts.summary <- sapply(dat, function(x)
tabulate(x, nbins = length(bins.summary)))
eventsPerYear <- apply(counts.summary, 2, sum)
freq.summary <- sweep(counts.summary, 2, STATS = eventsPerYear, FUN = "/")
rm(counts.summary)

} else {
freq.summary <- matrix(0, nrow = length(bins.summary), ncol = simTime$no.useyr)
eventsPerYear <- rep(0, simTime$no.useyr)
}

list(eventsPerYear = eventsPerYear, freq.summary = freq.summary)
})

#------------------------DAILY WEATHER
#TODO replace with Rsoilwat31::getWeatherData_folders
Expand Down

0 comments on commit 2829718

Please sign in to comment.