comparison insect_phenology_model.R @ 106:aa1890120743 draft

Uploaded
author greg
date Tue, 05 Dec 2017 11:21:30 -0500
parents d6e37828b094
children 2c2903e7ce68
comparison
equal deleted inserted replaced
105:375a2b8e7207 106:aa1890120743
200 legend_text <- c("Egg", "Nymph", "Adult") 200 legend_text <- c("Egg", "Nymph", "Adult")
201 columns <- c(4, 2, 1) 201 columns <- c(4, 2, 1)
202 } else if (chart_type == "pop_size_by_generation") { 202 } else if (chart_type == "pop_size_by_generation") {
203 title <- paste(insect, ": Total pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") 203 title <- paste(insect, ": Total pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
204 legend_text <- c("P", "F1", "F2") 204 legend_text <- c("P", "F1", "F2")
205 columns <- col=c(1, 2, 4) 205 columns <- c(1, 2, 4)
206 } else if (chart_type == "adult_pop_size_by_generation") { 206 } else if (chart_type == "adult_pop_size_by_generation") {
207 title <- paste(insect, ": Adult pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ") 207 title <- paste(insect, ": Adult pop. by generation :", location, ": Lat:", latitude, ":", start_date, "-", end_date, sep=" ")
208 legend_text <- c("P", "F1", "F2") 208 legend_text <- c("P", "F1", "F2")
209 columns <- col=c(1, 2, 4) 209 columns <- c(1, 2, 4)
210 } 210 }
211 plot(days, group1, main=title, type="l", ylim=c(0, maxval), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3) 211 plot(days, group1, main=title, type="l", ylim=c(0, maxval), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
212 legend("topleft", legend_text, lty=c(1, 1, 1), columns, cex=3) 212 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3)
213 lines(days, group2, lwd=2, lty=1, col=2) 213 lines(days, group2, lwd=2, lty=1, col=2)
214 lines(days, group3, lwd=2, lty=1, col=4) 214 lines(days, group3, lwd=2, lty=1, col=4)
215 axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) 215 axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
216 axis(2, cex.axis=3) 216 axis(2, cex.axis=3)
217 if (plot_std_error==1) { 217 if (plot_std_error==1) {
229 229
230 temperature_data_frame <- parse_input_data(opt$input, opt$num_days) 230 temperature_data_frame <- parse_input_data(opt$input, opt$num_days)
231 # All latitude values are the same, 231 # All latitude values are the same,
232 # so get the value from the first row. 232 # so get the value from the first row.
233 latitude <- temperature_data_frame$LATITUDE[1] 233 latitude <- temperature_data_frame$LATITUDE[1]
234
235 cat("Number of days: ", opt$num_days, "\n")
236 234
237 # Initialize matrices. 235 # Initialize matrices.
238 Eggs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 236 Eggs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
239 YoungNymphs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 237 YoungNymphs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
240 OldNymphs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications) 238 OldNymphs.replications <- matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications)
619 # Mean value for F2 adults. 617 # Mean value for F2 adults.
620 F2_adults <- apply(F2_adults.replications, 1, mean) 618 F2_adults <- apply(F2_adults.replications, 1, mean)
621 # Standard error for F2 adult. 619 # Standard error for F2 adult.
622 F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications) 620 F2_adults.std_error <- apply(F2_adults.replications, 1, sd) / sqrt(opt$replications)
623 621
622 # Display the total number of days in the Galaxy history item blurb.
623 cat("Number of days: ", opt$num_days, "\n")
624
624 dev.new(width=20, height=40) 625 dev.new(width=20, height=40)
625 626
626 # Start PDF device driver to save charts to output. 627 # Start PDF device driver to save charts to output.
627 pdf(file=opt$output, width=20, height=40, bg="white") 628 pdf(file=opt$output, width=20, height=40, bg="white")
628 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)) 629 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1))