| Previous changeset 84:1863d847749b (2017-11-15) Next changeset 86:5e1f9fde280c (2017-11-22) |
|
Commit message:
Uploaded |
|
modified:
insect_phenology_model.R |
| b |
| diff -r 1863d847749b -r 8f5c05be4efe insect_phenology_model.R --- a/insect_phenology_model.R Wed Nov 15 07:36:38 2017 -0500 +++ b/insect_phenology_model.R Wed Nov 22 13:02:16 2017 -0500 |
| [ |
| b'@@ -1,661 +1,646 @@\n-#!/usr/bin/env Rscript\r\n-\r\n-suppressPackageStartupMessages(library("optparse"))\r\n-\r\n-option_list <- list(\r\n- make_option(c("-a", "--adult_mort"), action="store", dest="adult_mort", type="integer", help="Adjustment rate for adult mortality"),\r\n- make_option(c("-b", "--adult_accum"), action="store", dest="adult_accum", type="integer", help="Adjustment of DD accumulation (old nymph->adult)"),\r\n- make_option(c("-c", "--egg_mort"), action="store", dest="egg_mort", type="integer", help="Adjustment rate for egg mortality"),\r\n- make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"),\r\n- make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),\r\n- make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),\r\n- make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"),\r\n- make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"),\r\n- make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),\r\n- make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),\r\n- make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),\r\n- make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),\r\n- make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"),\r\n- make_option(c("-t", "--se_plot"), action="store", dest="se_plot", help="Plot SE"),\r\n- make_option(c("-v", "--input"), action="store", dest="input", help="Temperature data for selected location"),\r\n- make_option(c("-y", "--young_nymph_accum"), action="store", dest="young_nymph_accum", type="integer", help="Adjustment of DD accumulation (egg->young nymph)")\r\n-)\r\n-\r\n-parser <- OptionParser(usage="%prog [options] file", option_list=option_list)\r\n-args <- parse_args(parser, positional_arguments=TRUE)\r\n-opt <- args$options\r\n-\r\n-parse_input_data = function(input_file, num_rows) {\r\n- # Read in the input temperature datafile into a data frame.\r\n- temperature_data_frame <- read.csv(file=input_file, header=T, strip.white=TRUE, sep=",")\r\n- if (dim(temperature_data_frame)[2] == 6) {\r\n- # The input data has the following 6 columns:\r\n- # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX\r\n- # Add a column containing the daylight length for each day.\r\n- temperature_data_frame <- add_daylight_length(temperature_data_frame, num_rows)\r\n- }\r\n- # Return the temperature_data_frame.\r\n- temperature_data_frame\r\n-}\r\n-\r\n-add_daylight_length = function(temperature_data_frame, num_rows) {\r\n- # Return a vector of daylight length (photoperido profile) for\r\n- # the number of days specified in the input temperature data\r\n- # (from Forsythe 1995).\r\n- p = 0.8333\r\n- latitude <- temperature_data_frame[1, 1]\r\n- daylight_length_vector <- NULL\r\n- for (i in 1:num_rows) {\r\n- # Get the day of the year from the current row\r\n- # of the temperature data for computation.\r\n- doy <- temperature_data_frame[i, 4]\r\n- theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)))\r\n- phi <- asin(0.39795 * cos(theta))\r\n- # Compute the length of daylight for the day of the year.\r\n- daylight_length_vector[i] <- 24 - (24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi))))\r\n- }\r\n- # Append daylight_length_vec'..b'lications), 1, sd)\n+# SE for eggs\n+se.se <- apply(S0.rep, 1, sd) / sqrt(opt$replications)\n+# SE value for P\n+g0.se <- apply(g0.rep, 1, sd) / sqrt(opt$replications)\n+# SE for F1\n+g1.se <- apply(g1.rep, 1, sd) / sqrt(opt$replications)\n+# SE for F2\n+g2.se <- apply(g2.rep, 1, sd) / sqrt(opt$replications)\n+# SE for P adult\n+g0a.se <- apply(g0a.rep, 1, sd) / sqrt(opt$replications)\n+# SE for F1 adult\n+g1a.se <- apply(g1a.rep, 1, sd) / sqrt(opt$replications)\n+# SE for F2 adult\n+g2a.se <- apply(g2a.rep, 1, sd) / sqrt(opt$replications)\n+\n+dev.new(width=20, height=30)\n+\n+# Start PDF device driver to save charts to output.\n+pdf(file=opt$output, width=20, height=30, bg="white")\n+\n+par(mar = c(5, 6, 4, 4), mfrow=c(3, 1))\n+\n+# Subfigure 1: population size by life stage\n+title <- paste("BSMB total population by life stage :", opt$location, ": Lat:", latitude, ":", start_date, "to", end_date, sep=" ")\n+plot(day.all, sa, main=title, type="l", ylim=c(0, max(se + se.se, sn + sn.se, sa + sa.se)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)\n+# Young and old nymphs.\n+lines(day.all, sn, lwd=2, lty=1, col=2)\n+# Eggs\n+lines(day.all, se, lwd=2, lty=1, col=4)\n+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"))\n+axis(2, cex.axis=3)\n+leg.text <- c("Egg", "Nymph", "Adult")\n+legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)\n+if (opt$se_plot == 1) {\n+ # Add SE lines to plot\n+ # SE for adults\n+ lines (day.all, sa + sa.se, lty=2)\n+ lines (day.all, sa - sa.se, lty=2) \n+ # SE for nymphs\n+ lines (day.all, sn + sn.se, col=2, lty=2)\n+ lines (day.all, sn - sn.se, col=2, lty=2)\n+ # SE for eggs\n+ lines (day.all, se + se.se, col=4, lty=2)\n+ lines (day.all, se - se.se, col=4, lty=2)\n+}\n+\n+# Subfigure 2: population size by generation\n+title <- paste("BSMB total population by generation :", opt$location, ": Lat:", latitude, ":", start_date, "to", end_date, sep=" ")\n+plot(day.all, g0, main=title, type="l", ylim=c(0, max(g2)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)\n+lines(day.all, g1, lwd = 2, lty = 1, col=2)\n+lines(day.all, g2, lwd = 2, lty = 1, col=4)\n+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"))\n+axis(2, cex.axis=3)\n+leg.text <- c("P", "F1", "F2")\n+legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)\n+if (opt$se_plot == 1) {\n+ # Add SE lines to plot\n+ # SE for adults\n+ lines (day.all, g0+g0.se, lty=2)\n+ lines (day.all, g0-g0.se, lty=2) \n+ # SE for nymphs\n+ lines (day.all, g1+g1.se, col=2, lty=2)\n+ lines (day.all, g1-g1.se, col=2, lty=2)\n+ # SE for eggs\n+ lines (day.all, g2+g2.se, col=4, lty=2)\n+ lines (day.all, g2-g2.se, col=4, lty=2)\n+}\n+\n+# Subfigure 3: adult population size by generation\n+title <- paste("BSMB adult population by generation :", opt$location, ": Lat:", latitude, ":", start_date, "to", end_date, sep=" ")\n+plot(day.all, g0a, ylim=c(0, max(g2a) + 100), main=title, type="l", axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)\n+lines(day.all, g1a, lwd = 2, lty = 1, col=2)\n+lines(day.all, g2a, lwd = 2, lty = 1, col=4)\n+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"))\n+axis(2, cex.axis=3)\n+leg.text <- c("P", "F1", "F2")\n+legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)\n+if (opt$se_plot == 1) {\n+ # Add SE lines to plot\n+ # SE for adults\n+ lines (day.all, g0a+g0a.se, lty=2)\n+ lines (day.all, g0a-g0a.se, lty=2) \n+ # SE for nymphs\n+ lines (day.all, g1a+g1a.se, col=2, lty=2)\n+ lines (day.all, g1a-g1a.se, col=2, lty=2)\n+ # SE for eggs\n+ lines (day.all, g2a+g2a.se, col=4, lty=2)\n+ lines (day.all, g2a-g2a.se, col=4, lty=2)\n+}\n+\n+# Turn off device driver to flush output.\n+dev.off()\n' |