Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 125:24f389a2dd93 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Tue, 07 Aug 2018 12:58:20 -0400
parents 534658644efe
children 397c24c1adc9
comparison
equal deleted inserted replaced
124:534658644efe 125:24f389a2dd93
4 4
5 option_list <- list( 5 option_list <- list(
6 make_option(c("--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"), 6 make_option(c("--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"),
7 make_option(c("--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"), 7 make_option(c("--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"),
8 make_option(c("--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"), 8 make_option(c("--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"),
9 make_option(c("--end_date"), action="store", dest="end_date", default=NULL, help="End date for custom date interval"),
10 make_option(c("--input_norm"), action="store", dest="input_norm", help="30 year normals temperature data for selected station"), 9 make_option(c("--input_norm"), action="store", dest="input_norm", help="30 year normals temperature data for selected station"),
11 make_option(c("--input_ytd"), action="store", dest="input_ytd", default=NULL, help="Year-to-date temperature data for selected location"), 10 make_option(c("--input_ytd"), action="store", dest="input_ytd", default=NULL, help="Year-to-date temperature data for selected location"),
12 make_option(c("--insect"), action="store", dest="insect", help="Insect name"), 11 make_option(c("--insect"), action="store", dest="insect", help="Insect name"),
13 make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"), 12 make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"),
14 make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"), 13 make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"),
23 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), 22 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
24 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), 23 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
25 make_option(c("--plot_generations_separately"), action="store", dest="plot_generations_separately", help="Plot Plot P, F1 and F2 as separate lines or pool across them"), 24 make_option(c("--plot_generations_separately"), action="store", dest="plot_generations_separately", help="Plot Plot P, F1 and F2 as separate lines or pool across them"),
26 make_option(c("--plot_std_error"), action="store", dest="plot_std_error", help="Plot Standard error"), 25 make_option(c("--plot_std_error"), action="store", dest="plot_std_error", help="Plot Standard error"),
27 make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"), 26 make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
28 make_option(c("--start_date"), action="store", dest="start_date", default=NULL, help="Start date for custom date interval"), 27 make_option(c("--script_dir"), action="store", dest="script_dir", help="R script source directory"),
29 make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)") 28 make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)")
30 ) 29 )
31 30
32 parser <- OptionParser(usage="%prog [options] file", option_list=option_list); 31 parser <- OptionParser(usage="%prog [options] file", option_list=option_list);
33 args <- parse_args(parser, positional_arguments=TRUE); 32 args <- parse_args(parser, positional_arguments=TRUE);
34 opt <- args$options; 33 opt <- args$options;
35 34
36 add_daylight_length = function(temperature_data_frame) { 35 add_daylight_length = function(temperature_data_frame) {
37 # Return temperature_data_frame with an added column 36 # Return temperature_data_frame with an added column
38 # of daylight length (photoperido profile). 37 # of daylight length (photoperiod profile).
39 num_rows = dim(temperature_data_frame)[1]; 38 num_rows = dim(temperature_data_frame)[1];
40 # From Forsythe 1995. 39 # From Forsythe 1995.
41 p = 0.8333; 40 p = 0.8333;
42 latitude = temperature_data_frame$LATITUDE[1]; 41 latitude = temperature_data_frame$LATITUDE[1];
43 daylight_length_vector = NULL; 42 daylight_length_vector = NULL;
64 # Reset the column names with the additional column for later access. 63 # Reset the column names with the additional column for later access.
65 colnames(data_frame) = append(current_column_names, new_column_name); 64 colnames(data_frame) = append(current_column_names, new_column_name);
66 return(data_frame); 65 return(data_frame);
67 } 66 }
68 67
69 extract_date_interval_rows = function(df, start_date, end_date) { 68 from_30_year_normals = function(norm_data_frame, start_date_doy, end_date_doy, year) {
70 date_interval_rows = df[df$DATE >= start_date & df$DATE <= end_date];
71 return(date_interval_rows);
72 }
73
74 from_30_year_normals = function(temperature_data_frame, norm_data_frame, start_date_doy, end_date_doy) {
75 # The data we want is fully contained within the 30 year normals data. 69 # The data we want is fully contained within the 30 year normals data.
76 first_norm_row = which(norm_data_frame$DOY==start_date_doy); 70 first_norm_row = which(norm_data_frame$DOY==start_date_doy);
77 last_norm_row = which(norm_data_frame$DOY==end_date_doy); 71 last_norm_row = which(norm_data_frame$DOY==end_date_doy);
78 norm_data_frame_rows = last_norm_row - first_norm_row; 72 # Add 1 to the number of rows to ensure that the end date is included.
73 tmp_data_frame_rows = last_norm_row - first_norm_row + 1;
74 tmp_data_frame = get_new_temperature_data_frame(nrow=tmp_data_frame_rows);
79 j = 0; 75 j = 0;
80 for (i in first_norm_row:last_norm_row) { 76 for (i in first_norm_row:last_norm_row) {
81 j = j + 1; 77 j = j + 1;
82 temperature_data_frame[j,] = get_next_normals_row(norm_data_frame, year, i); 78 tmp_data_frame[j,] = get_next_normals_row(norm_data_frame, year, i);
83 } 79 }
84 return (temperature_data_frame); 80 return (tmp_data_frame);
85 } 81 }
86 82
87 get_file_path = function(life_stage, base_name, life_stage_nymph=NULL, life_stage_adult=NULL) { 83 get_new_norm_data_frame = function(is_leap_year, input_norm=NULL, nrow=0) {
88 if (!is.null(life_stage_nymph)) { 84 # The input_norm data has the following 10 columns:
89 lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph); 85 # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX
90 file_name = paste(lsi, tolower(life_stage_nymph), base_name, sep="_"); 86 column_names = c("STATIONID", "LATITUDE","LONGITUDE", "ELEV_M", "NAME", "ST", "MMDD", "DOY", "TMIN", "TMAX");
91 } else if (!is.null(life_stage_adult)) { 87 if (is.null(input_norm)) {
92 lsi = get_life_stage_index(life_stage, life_stage_adult=life_stage_adult); 88 norm_data_frame = data.frame(matrix(ncol=10, nrow));
93 file_name = paste(lsi, tolower(life_stage_adult), base_name, sep="_"); 89 # Set the norm_data_frame column names for access.
90 colnames(norm_data_frame) = column_names;
94 } else { 91 } else {
95 lsi = get_life_stage_index(life_stage); 92 norm_data_frame = read.csv(file=input_norm, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
96 file_name = paste(lsi, base_name, sep="_"); 93 # Set the norm_data_frame column names for access.
97 } 94 colnames(norm_data_frame) = column_names;
98 file_path = paste("output_plots_dir", file_name, sep="/"); 95 if (!is_leap_year) {
99 return(file_path); 96 # All normals data includes Feb 29 which is row 60 in
100 } 97 # the data, so delete that row if we're not in a leap year.
101 98 norm_data_frame = norm_data_frame[-c(60),];
102 get_life_stage_index = function(life_stage, life_stage_nymph=NULL, life_stage_adult=NULL) { 99 # Since we've removed row 60, we need to subtract 1 from
103 # Name collection elements so that they 100 # each value in the DOY column of the data frame starting
104 # are displayed in logical order. 101 # with the 60th row.
105 if (life_stage=="Egg") { 102 num_rows = dim(norm_data_frame)[1];
106 lsi = "01"; 103 for (i in 60:num_rows) {
107 } else if (life_stage=="Nymph") { 104 leap_year_doy = norm_data_frame$DOY[i];
108 if (life_stage_nymph=="Young") { 105 non_leap_year_doy = leap_year_doy - 1;
109 lsi = "02"; 106 norm_data_frame$DOY[i] = non_leap_year_doy;
110 } else if (life_stage_nymph=="Old") { 107 }
111 lsi = "03"; 108 }
112 } else if (life_stage_nymph=="Total") { 109 }
113 lsi="04"; 110 return (norm_data_frame);
114 } 111 }
115 } else if (life_stage=="Adult") { 112
116 if (life_stage_adult=="Pre-vittelogenic") { 113 get_new_temperature_data_frame = function(input_ytd=NULL, nrow=0) {
117 lsi = "05"; 114 # The input_ytd data has the following 6 columns:
118 } else if (life_stage_adult=="Vittelogenic") { 115 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
119 lsi = "06"; 116 if (is.null(input_ytd)) {
120 } else if (life_stage_adult=="Diapausing") { 117 temperature_data_frame = data.frame(matrix(ncol=6, nrow));
121 lsi = "07"; 118 } else {
122 } else if (life_stage_adult=="Total") { 119 temperature_data_frame = read.csv(file=input_ytd, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
123 lsi = "08"; 120 }
124 } 121 colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
125 } else if (life_stage=="Total") { 122 return(temperature_data_frame);
126 lsi = "09";
127 }
128 return(lsi);
129 }
130
131 get_mean_and_std_error = function(p_replications, f1_replications, f2_replications) {
132 # P mean.
133 p_m = apply(p_replications, 1, mean);
134 # P standard error.
135 p_se = apply(p_replications, 1, sd) / sqrt(opt$replications);
136 # F1 mean.
137 f1_m = apply(f1_replications, 1, mean);
138 # F1 standard error.
139 f1_se = apply(f1_replications, 1, sd) / sqrt(opt$replications);
140 # F2 mean.
141 f2_m = apply(f2_replications, 1, mean);
142 # F2 standard error.
143 f2_se = apply(f2_replications, 1, sd) / sqrt(opt$replications);
144 return(list(p_m, p_se, f1_m, f1_se, f2_m, f2_se))
145 } 123 }
146 124
147 get_next_normals_row = function(norm_data_frame, year, index) { 125 get_next_normals_row = function(norm_data_frame, year, index) {
148 # Return the next 30 year normals row formatted 126 # Return the next 30 year normals row formatted
149 # appropriately for the year-to-date data. 127 # appropriately for the year-to-date data.
227 averages = sum(dh) / 24; 205 averages = sum(dh) / 24;
228 } 206 }
229 return(c(curr_mean_temp, averages)) 207 return(c(curr_mean_temp, averages))
230 } 208 }
231 209
232 get_tick_index = function(index, last_tick, ticks, month_labels) {
233 # The R code tries hard not to draw overlapping tick labels, and so
234 # will omit labels where they would abut or overlap previously drawn
235 # labels. This can result in, for example, every other tick being
236 # labelled. We'll keep track of the last tick to make sure all of
237 # the month labels are displayed, and missing ticks are restricted
238 # to Sundays which have no labels anyway.
239 if (last_tick==0) {
240 return(length(ticks)+1);
241 }
242 last_saved_tick = ticks[[length(ticks)]];
243 if (index-last_saved_tick<3) {
244 last_saved_month = month_labels[[length(month_labels)]];
245 if (last_saved_month=="") {
246 # We're safe overwriting a tick
247 # with no label (i.e., a Sunday tick).
248 return(length(ticks));
249 } else {
250 # Don't eliminate a Month label.
251 return(NULL);
252 }
253 }
254 return(length(ticks)+1);
255 }
256
257 get_total_days = function(is_leap_year) {
258 # Get the total number of days in the current year.
259 if (is_leap_year) {
260 return(366);
261 } else {
262 return(365);
263 }
264 }
265
266 get_x_axis_ticks_and_labels = function(temperature_data_frame, prepend_end_doy_norm, append_start_doy_norm, date_interval) {
267 # Generate a list of ticks and labels for plotting the
268 # x axis. There are several scenarios that affect this.
269 # 1. If date_interval is TRUE:
270 # a.
271 if (prepend_end_doy_norm > 0) {
272 prepend_end_norm_row = which(temperature_data_frame$DOY==prepend_end_doy_norm);
273 } else {
274 prepend_end_norm_row = 0;
275 }
276 if (append_start_doy_norm > 0) {
277 append_start_norm_row = which(temperature_data_frame$DOY==append_start_doy_norm);
278 } else {
279 append_start_norm_row = 0;
280 }
281 num_rows = dim(temperature_data_frame)[1];
282 month_labels = list();
283 ticks = list();
284 current_month_label = NULL;
285 last_tick = 0;
286 for (i in 1:num_rows) {
287 # We're plotting the entire year, so ticks will
288 # occur on Sundays and the first of each month.
289 if (i == prepend_end_norm_row) {
290 # Add a tick for the end of the 30 year normnals data
291 # that was prepended to the year-to-date data.
292 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
293 ticks[tick_index] = i;
294 month_labels[tick_index] = "End prepended 30 year normals";
295 last_tick = i;
296 } else if (i == append_start_doy_norm) {
297 # Add a tick for the start of the 30 year normnals data
298 # that was appended to the year-to-date data.
299 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
300 ticks[tick_index] = i;
301 month_labels[tick_index] = "Start appended 30 year normals";
302 last_tick = i;
303 } else if (i==num_rows) {
304 # Add a tick for the last day of the year.
305 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
306 ticks[tick_index] = i;
307 month_labels[tick_index] = "";
308 last_tick = i;
309 } else {
310 # Get the year and month from the date which
311 # has the format YYYY-MM-DD.
312 date = format(temperature_data_frame$DATE[i]);
313 # Get the month label.
314 items = strsplit(date, "-")[[1]];
315 month = items[2];
316 month_label = month.abb[as.integer(month)];
317 if (!identical(current_month_label, month_label)) {
318 # Add a tick for the month.
319 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
320 ticks[tick_index] = i;
321 month_labels[tick_index] = month_label;
322 current_month_label = month_label;
323 last_tick = i;
324 }
325 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
326 if (!is.null(tick_index)) {
327 if (date_interval) {
328 # Add a tick for every day.
329 ticks[tick_index] = i;
330 # Add a blank month label so it is not displayed.
331 month_labels[tick_index] = "";
332 last_tick = i;
333 } else {
334 # Get the day.
335 day = weekdays(as.Date(date));
336 if (day=="Sunday") {
337 # Add a tick if we're on a Sunday.
338 ticks[tick_index] = i;
339 # Add a blank month label so it is not displayed.
340 month_labels[tick_index] = "";
341 last_tick = i;
342 }
343 }
344 }
345 }
346 }
347 return(list(ticks, month_labels));
348 }
349
350 is_leap_year = function(date_str) { 210 is_leap_year = function(date_str) {
351 # Extract the year from the date_str. 211 # Extract the year from the date_str.
352 date = format(date_str); 212 date = format(date_str);
353 items = strsplit(date, "-")[[1]]; 213 items = strsplit(date, "-")[[1]];
354 year = as.integer(items[1]); 214 year = as.integer(items[1]);
395 parse_input_data = function(input_ytd, input_norm, location, start_date, end_date) { 255 parse_input_data = function(input_ytd, input_norm, location, start_date, end_date) {
396 # The end DOY for norm data prepended to ytd data. 256 # The end DOY for norm data prepended to ytd data.
397 prepend_end_doy_norm = 0; 257 prepend_end_doy_norm = 0;
398 # The start DOY for norm data appended to ytd data. 258 # The start DOY for norm data appended to ytd data.
399 append_start_doy_norm = 0; 259 append_start_doy_norm = 0;
260 if (is.null(start_date) && is.null(end_date)) {
261 # We're not dealing with a date interval.
262 date_interval = FALSE;
263 if (is.null(input_ytd)) {
264 # Base all dates on the current date since 30 year
265 # normals data does not include any dates.
266 year = format(Sys.Date(), "%Y");
267 }
268 } else {
269 date_interval = TRUE;
270 year = get_year_from_date(start_date);
271 # Get the DOY for start_date and end_date.
272 start_date_doy = as.integer(strftime(start_date, format="%j"));
273 end_date_doy = as.integer(strftime(end_date, format="%j"));
274 }
400 if (is.null(input_ytd)) { 275 if (is.null(input_ytd)) {
401 # We're processing only the 30 year normals data. 276 # We're processing only the 30 year normals data.
402 processing_year_to_date_data = FALSE; 277 processing_year_to_date_data = FALSE;
278 if (is.null(start_date) && is.null(end_date)) {
279 # We're processing the entire year, so we can
280 # set the start_date to Jan 1.
281 start_date = paste(year, "01", "01", sep="-");
282 }
403 } else { 283 } else {
404 processing_year_to_date_data = TRUE; 284 processing_year_to_date_data = TRUE;
405 } 285 # Read the input_ytd temperature data file into a data frame.
406 if (is.null(start_date) && is.null(end_date)) { 286 temperature_data_frame = get_new_temperature_data_frame(input_ytd=input_ytd);
407 # We're processing the entire year, possibly merging 287 num_ytd_rows = dim(temperature_data_frame)[1];
408 # data from input_norm with data from input_ytd. 288 if (!date_interval) {
409 date_interval = FALSE; 289 start_date = temperature_data_frame$DATE[1];
410 } else { 290 year = get_year_from_date(start_date);
411 date_interval = TRUE; 291 }
412 # Get the DOY for start_date and end_date. 292 }
413 start_date_doy = strftime(start_date, format="%j"); 293 # See if we're in a leap year.
414 end_date_doy = strftime(end_date, format="%j"); 294 is_leap_year = is_leap_year(start_date);
415 }
416 # Read the input_norm temperature datafile into a data frame. 295 # Read the input_norm temperature datafile into a data frame.
417 # The input_norm data has the following 10 columns: 296 norm_data_frame = get_new_norm_data_frame(is_leap_year, input_norm=input_norm);
418 # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX
419 norm_data_frame = read.csv(file=input_norm, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
420 # Set the norm_data_frame column names for access.
421 colnames(norm_data_frame) = c("STATIONID", "LATITUDE","LONGITUDE", "ELEV_M", "NAME", "ST", "MMDD", "DOY", "TMIN", "TMAX");
422 if (processing_year_to_date_data) { 297 if (processing_year_to_date_data) {
423 # Read the input_ytd temperature data file into a data frame.
424 # The input_ytd data has the following 6 columns:
425 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
426 temperature_data_frame = read.csv(file=input_ytd, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
427 # Set the temperature_data_frame column names for access.
428 colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
429 if (date_interval) { 298 if (date_interval) {
430 # We're plotting a date interval. 299 # We're plotting a date interval.
431 start_date_ytd_row = which(temperature_data_frame$DATE==start_date); 300 start_date_ytd_row = which(temperature_data_frame$DATE==start_date);
432 if (start_date_ytd_row > 0) { 301 if (length(start_date_ytd_row) > 0) {
433 # The start date is contained within the input_ytd data. 302 # The start date is contained within the input_ytd data.
303 start_date_ytd_row = start_date_ytd_row[1];
434 start_doy_ytd = as.integer(temperature_data_frame$DOY[start_date_ytd_row]); 304 start_doy_ytd = as.integer(temperature_data_frame$DOY[start_date_ytd_row]);
435 } else { 305 } else {
436 # The start date is contained within the input_norm data. 306 # The start date is contained within the input_norm data.
437 start_date_norm_row = which(norm_data_frame$DATE==start_date); 307 start_date_ytd_row = 0;
308 start_date_norm_row = which(norm_data_frame$DOY==start_date_doy);
438 } 309 }
439 end_date_ytd_row = which(temperature_data_frame$DATE==end_date); 310 end_date_ytd_row = which(temperature_data_frame$DATE==end_date);
440 if (end_date_ytd_row > 0) { 311 if (length(end_date_ytd_row) > 0) {
312 end_date_ytd_row = end_date_ytd_row[1];
441 # The end date is contained within the input_ytd data. 313 # The end date is contained within the input_ytd data.
442 end_doy_ytd = as.integer(temperature_data_frame$DOY[end_date_ytd_row]); 314 end_doy_ytd = as.integer(temperature_data_frame$DOY[end_date_ytd_row]);
443 } 315 } else {
444 date_str = start_date; 316 end_date_ytd_row = 0;
445 # Extract the year from the start date. 317 }
446 date_str_items = strsplit(date_str, "-")[[1]];
447 year = date_str_items[1];
448 } else { 318 } else {
449 # We're plotting an entire year. 319 # We're plotting an entire year.
450 # Get the number of days contained in temperature_data_frame.
451 num_rows = dim(temperature_data_frame)[1];
452 # Get the start date and end date from temperature_data_frame. 320 # Get the start date and end date from temperature_data_frame.
453 start_date_ytd_row = 1; 321 start_date_ytd_row = 1;
454 # Temporarily set start_date to get the year. 322 # Temporarily set start_date to get the year.
455 start_date = temperature_data_frame$DATE[1]; 323 start_date = temperature_data_frame$DATE[1];
456 end_date_ytd_row = num_rows; 324 end_date_ytd_row = num_ytd_rows;
457 end_date = temperature_data_frame$DATE[num_rows]; 325 end_date = temperature_data_frame$DATE[num_ytd_rows];
458 date_str = format(start_date); 326 date_str = format(start_date);
459 # Extract the year from the start date. 327 # Extract the year from the start date.
460 date_str_items = strsplit(date_str, "-")[[1]]; 328 date_str_items = strsplit(date_str, "-")[[1]];
461 # Get the year. 329 # Get the year.
462 year = date_str_items[1]; 330 year = date_str_items[1];
464 start_date = paste(year, "01", "01", sep="-"); 332 start_date = paste(year, "01", "01", sep="-");
465 # Properly set the end_date to be Dec 31 of the year. 333 # Properly set the end_date to be Dec 31 of the year.
466 end_date = paste(year, "12", "31", sep="-"); 334 end_date = paste(year, "12", "31", sep="-");
467 # Save the first DOY to later check if start_date is Jan 1. 335 # Save the first DOY to later check if start_date is Jan 1.
468 start_doy_ytd = as.integer(temperature_data_frame$DOY[1]); 336 start_doy_ytd = as.integer(temperature_data_frame$DOY[1]);
469 end_doy_ytd = as.integer(temperature_data_frame$DOY[num_rows]); 337 end_doy_ytd = as.integer(temperature_data_frame$DOY[num_ytd_rows]);
470 } 338 }
471 } else { 339 } else {
472 # We're processing only the 30 year normals data, so create an empty 340 # We're processing only the 30 year normals data, so create an empty
473 # data frame for containing temperature data after it is converted 341 # data frame for containing temperature data after it is converted
474 # from the 30 year normals format to the year-to-date format. 342 # from the 30 year normals format to the year-to-date format.
475 temperature_data_frame = data.frame(matrix(ncol=6, nrow=0)); 343 temperature_data_frame = get_new_temperature_data_frame();
476 colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
477 if (date_interval) { 344 if (date_interval) {
478 # We're plotting a date interval. 345 # We're plotting a date interval.
479 # Extract the year, month and day from the start date. 346 # Extract the year, month and day from the start date.
480 start_date_str = format(start_date); 347 start_date_str = format(start_date);
481 start_date_str_items = strsplit(start_date_str, "-")[[1]]; 348 start_date_str_items = strsplit(start_date_str, "-")[[1]];
489 end_date_month = end_date_str_items[2]; 356 end_date_month = end_date_str_items[2];
490 end_date_day = end_date_str_items[3]; 357 end_date_day = end_date_str_items[3];
491 end_date = paste(year, end_date_month, end_date_day, sep="-"); 358 end_date = paste(year, end_date_month, end_date_day, sep="-");
492 } else { 359 } else {
493 # We're plotting an entire year. 360 # We're plotting an entire year.
494 # Base all dates on the current date since 30 year
495 # normals data does not include any dates.
496 year = format(Sys.Date(), "%Y");
497 start_date = paste(year, "01", "01", sep="-"); 361 start_date = paste(year, "01", "01", sep="-");
498 end_date = paste(year, "12", "31", sep="-"); 362 end_date = paste(year, "12", "31", sep="-");
499 } 363 }
500 }
501 # See if we're in a leap year.
502 is_leap_year = is_leap_year(start_date);
503 # All normals data includes Feb 29 which is row 60 in
504 # the data, so delete that row if we're not in a leap year.
505 if (!is_leap_year) {
506 norm_data_frame = norm_data_frame[-c(60),];
507 } 364 }
508 # Set the location to be the station name if the user elected not to enter it. 365 # Set the location to be the station name if the user elected not to enter it.
509 if (is.null(location) | length(location) == 0) { 366 if (is.null(location) | length(location) == 0) {
510 location = norm_data_frame$NAME[1]; 367 location = norm_data_frame$NAME[1];
511 } 368 }
526 first_ytd_doy = temperature_data_frame$DOY[1]; 383 first_ytd_doy = temperature_data_frame$DOY[1];
527 # End DOY of input_norm data prepended to input_ytd. 384 # End DOY of input_norm data prepended to input_ytd.
528 prepend_end_doy_norm = first_ytd_doy - 1; 385 prepend_end_doy_norm = first_ytd_doy - 1;
529 # Get the number of rows for the restricted date interval 386 # Get the number of rows for the restricted date interval
530 # that are contained in temperature_data_frame. 387 # that are contained in temperature_data_frame.
531 temperature_data_frame_rows = end_date_ytd_row; 388 num_temperature_data_frame_rows = end_date_ytd_row;
532 # Get the last row needed from the 30 year normals data. 389 # Get the last row needed from the 30 year normals data.
533 last_norm_row = which(norm_data_frame$DOY==prepend_end_doy_norm); 390 last_norm_row = which(norm_data_frame$DOY==prepend_end_doy_norm);
534 # Get the number of rows for the restricted date interval 391 # Get the number of rows for the restricted date interval
535 # that are contained in norm_data_frame. 392 # that are contained in norm_data_frame.
536 norm_data_frame_rows = last_norm_row - first_norm_row; 393 num_norm_data_frame_rows = last_norm_row - first_norm_row;
537 # Create a temporary data frame to contain the 30 year normals 394 # Create a temporary data frame to contain the 30 year normals
538 # data from the start date to the date immediately prior to the 395 # data from the start date to the date immediately prior to the
539 # first row of the input_ytd data. 396 # first row of the input_ytd data.
540 tmp_norm_data_frame = data.frame(matrix(ncol=6, nrow=temperature_data_frame_rows+norm_data_frame_rows)); 397 tmp_norm_data_frame = get_new_temperature_data_frame(nrow=num_temperature_data_frame_rows+num_norm_data_frame_rows);
398 j = 1;
541 for (i in first_norm_row:last_norm_row) { 399 for (i in first_norm_row:last_norm_row) {
542 # Populate the temp_data_frame row with 400 # Populate the temp_data_frame row with
543 # values from norm_data_frame. 401 # values from norm_data_frame.
544 tmp_norm_data_frame[i,] = get_next_normals_row(norm_data_frame, year, i); 402 tmp_norm_data_frame[j,] = get_next_normals_row(norm_data_frame, year, i);
403 j = j + 1;
545 } 404 }
546 # Create a second temporary data frame containing the 405 # Create a second temporary data frame containing the
547 # appropriate rows from temperature_data_frame. 406 # appropriate rows from temperature_data_frame.
548 tmp_temperature_data_frame = temperature_data_frame[1:first_norm_row-1,]; 407 tmp_temperature_data_frame = temperature_data_frame[1:num_temperature_data_frame_rows,];
549 # Merge the 2 temporary data frames. 408 # Merge the 2 temporary data frames.
550 temperature_data_frame = rbind(tmp_norm_data_frame, tmp_temperature_data_frame); 409 temperature_data_frame = rbind(tmp_norm_data_frame, tmp_temperature_data_frame);
551 } else if (start_date_ytd_row > 0 & end_date_ytd_row == 0) { 410 } else if (start_date_ytd_row > 0 & end_date_ytd_row == 0) {
552 # The date interval starts in input_ytd and ends in input_norm, 411 # The date interval starts in input_ytd and ends in input_norm,
553 # so append appropriate rows from input_norm to appropriate rows 412 # so append appropriate rows from input_norm to appropriate rows
554 # from input_ytd. 413 # from input_ytd. First, get the number of rows for the restricted
555 num_rows = dim(temperature_data_frame)[1]; 414 # date interval that are contained in temperature_data_frame.
556 # Get the number of rows for the restricted date interval 415 num_temperature_data_frame_rows = num_ytd_rows - start_date_ytd_row + 1;
557 # that are contained in temperature_data_frame.
558 temperature_data_frame_rows = num_rows - start_date_ytd_row
559 # Get the DOY of the last row in the input_ytd data. 416 # Get the DOY of the last row in the input_ytd data.
560 last_ytd_doy = temperature_data_frame$DOY[num_rows]; 417 last_ytd_doy = temperature_data_frame$DOY[num_ytd_rows];
561 # Get the DOYs for the first and last rows from norm_data_frame 418 # Get the DOYs for the first and last rows from norm_data_frame
562 # that will be appended to temperature_data_frame. 419 # that will be appended to temperature_data_frame.
563 append_start_doy_norm = last_ytd_doy + 1; 420 append_start_doy_norm = last_ytd_doy + 1;
564 # Get the row from norm_data_frame containing first_norm_doy. 421 # Get the row from norm_data_frame containing first_norm_doy.
565 first_norm_row = which(norm_data_frame$DOY == append_start_doy_norm); 422 first_norm_row = which(norm_data_frame$DOY == append_start_doy_norm);
566 # Get the row from norm_data_frame containing end_date_doy. 423 # Get the row from norm_data_frame containing end_date_doy.
567 last_norm_row = which(norm_data_frame$DOY == end_date_doy); 424 last_norm_row = which(norm_data_frame$DOY == end_date_doy);
568 # Get the number of rows for the restricted date interval 425 # Get the number of rows for the restricted date interval
569 # that are contained in norm_data_frame. 426 # that are contained in norm_data_frame.
570 norm_data_frame_rows = last_norm_row - first_norm_row; 427 num_norm_data_frame_rows = last_norm_row - first_norm_row;
571 # Create a temporary data frame to contain the data 428 # Create a temporary data frame to contain the data
572 # taken from both temperateu_data_frame and norm_data_frame 429 # taken from both temperature_data_frame and norm_data_frame
573 # for the date interval. 430 # for the date interval.
574 tmp_data_frame = data.frame(matrix(ncol=6, nrow=temperature_data_frame_rows+norm_data_frame_rows)); 431 tmp_data_frame = get_new_temperature_data_frame(nrow=num_temperature_data_frame_rows+num_norm_data_frame_rows);
575 # Populate tmp_data_frame with the appropriate rows from temperature_data_frame. 432 # Populate tmp_data_frame with the appropriate rows from temperature_data_frame.
576 tmp_data_frame[temperature_data_frame_rows,] = temperature_data_frame[start_date_ytd_row:temperature_data_frame_rows,]; 433 j = start_date_ytd_row;
434 for (i in 1:num_temperature_data_frame_rows) {
435 tmp_data_frame[i,] = temperature_data_frame[j,];
436 j = j + 1;
437 }
577 # Apppend the appropriate rows from norm_data_frame to tmp_data_frame. 438 # Apppend the appropriate rows from norm_data_frame to tmp_data_frame.
578 for (i in first_norm_row:last_norm_row) { 439 current_iteration = num_temperature_data_frame_rows + 1;
579 tmp_data_frame[i,] = get_next_normals_row(norm_data_frame, year, i); 440 num_iterations = current_iteration + num_norm_data_frame_rows;
441 j = first_norm_row;
442 for (i in current_iteration:num_iterations) {
443 tmp_data_frame[i,] = get_next_normals_row(norm_data_frame, year, j);
444 j = j + 1;
580 } 445 }
581 temperature_data_frame = tmp_data_frame[,]; 446 temperature_data_frame = tmp_data_frame[,];
582 } else if (start_date_ytd_row == 0 & end_date_ytd_row == 0) { 447 } else if (start_date_ytd_row == 0 & end_date_ytd_row == 0) {
583 # The date interval is contained witin input_norm. 448 # The date interval is contained witin input_norm.
584 temperature_data_frame = from_30_year_normals(temperature_data_frame, norm_data_frame, start_date_doy, end_date_doy); 449 temperature_data_frame = from_30_year_normals(norm_data_frame, start_date_doy, end_date_doy, year);
585 } 450 }
586 } else { 451 } else {
587 # We're plotting an entire year. 452 # We're plotting an entire year.
588 if (start_doy_ytd > 1) { 453 if (start_doy_ytd > 1) {
589 # The input_ytd data starts after Jan 1, so prepend 454 # The input_ytd data starts after Jan 1, so prepend
614 } 479 }
615 } else { 480 } else {
616 # We're processing only the 30 year normals data. 481 # We're processing only the 30 year normals data.
617 if (date_interval) { 482 if (date_interval) {
618 # Populate temperature_data_frame from norm_data_frame. 483 # Populate temperature_data_frame from norm_data_frame.
619 temperature_data_frame = from_30_year_normals(temperature_data_frame, norm_data_frame, start_date_doy, end_date_doy); 484 temperature_data_frame = from_30_year_normals(norm_data_frame, start_date_doy, end_date_doy, year);
620 } else { 485 } else {
621 total_days = get_total_days(is_leap_year); 486 total_days = get_total_days(is_leap_year);
622 for (i in 1:total_days) { 487 for (i in 1:total_days) {
623 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, i); 488 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, i);
624 } 489 }
627 # Add a column containing the daylight length for each day. 492 # Add a column containing the daylight length for each day.
628 temperature_data_frame = add_daylight_length(temperature_data_frame); 493 temperature_data_frame = add_daylight_length(temperature_data_frame);
629 return(list(temperature_data_frame, start_date, end_date, prepend_end_doy_norm, append_start_doy_norm, is_leap_year, location)); 494 return(list(temperature_data_frame, start_date, end_date, prepend_end_doy_norm, append_start_doy_norm, is_leap_year, location));
630 } 495 }
631 496
632 render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval, 497 # Import the shared utility functions.
633 replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL, 498 utils_path <- paste(opt$script_dir, "utils.R", sep="/");
634 life_stages_adult=NULL, life_stages_nymph=NULL) { 499 source(utils_path);
635 if (chart_type=="pop_size_by_life_stage") {
636 if (life_stage=="Total") {
637 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
638 legend_text = c("Egg", "Nymph", "Adult");
639 columns = c(4, 2, 1);
640 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
641 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
642 lines(days, group2, lwd=2, lty=1, col=2);
643 lines(days, group3, lwd=2, lty=1, col=4);
644 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
645 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
646 if (plot_std_error=="yes") {
647 # Standard error for group.
648 lines(days, group+group_std_error, lty=2);
649 lines(days, group-group_std_error, lty=2);
650 # Standard error for group2.
651 lines(days, group2+group2_std_error, col=2, lty=2);
652 lines(days, group2-group2_std_error, col=2, lty=2);
653 # Standard error for group3.
654 lines(days, group3+group3_std_error, col=4, lty=2);
655 lines(days, group3-group3_std_error, col=4, lty=2);
656 }
657 } else {
658 if (life_stage=="Egg") {
659 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
660 legend_text = c(life_stage);
661 columns = c(4);
662 } else if (life_stage=="Nymph") {
663 stage = paste(life_stages_nymph, "Nymph Pop :", sep=" ");
664 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
665 legend_text = c(paste(life_stages_nymph, life_stage, sep=" "));
666 columns = c(2);
667 } else if (life_stage=="Adult") {
668 stage = paste(life_stages_adult, "Adult Pop", sep=" ");
669 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
670 legend_text = c(paste(life_stages_adult, life_stage, sep=" "));
671 columns = c(1);
672 }
673 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
674 legend("topleft", legend_text, lty=c(1), col="black", cex=3);
675 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
676 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
677 if (plot_std_error=="yes") {
678 # Standard error for group.
679 lines(days, group+group_std_error, lty=2);
680 lines(days, group-group_std_error, lty=2);
681 }
682 }
683 } else if (chart_type=="pop_size_by_generation") {
684 if (life_stage=="Total") {
685 title_str = ": Total Pop by Gen :";
686 } else if (life_stage=="Egg") {
687 title_str = ": Egg Pop by Gen :";
688 } else if (life_stage=="Nymph") {
689 title_str = paste(":", life_stages_nymph, "Nymph Pop by Gen", ":", sep=" ");
690 } else if (life_stage=="Adult") {
691 title_str = paste(":", life_stages_adult, "Adult Pop by Gen", ":", sep=" ");
692 }
693 title = paste(insect, ": Reps", replications, title_str, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
694 legend_text = c("P", "F1", "F2");
695 columns = c(1, 2, 4);
696 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3);
697 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
698 lines(days, group2, lwd=2, lty=1, col=2);
699 lines(days, group3, lwd=2, lty=1, col=4);
700 axis(side=1, at=ticks, labels=date_labels, las=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
701 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
702 if (plot_std_error=="yes") {
703 # Standard error for group.
704 lines(days, group+group_std_error, lty=2);
705 lines(days, group-group_std_error, lty=2);
706 # Standard error for group2.
707 lines(days, group2+group2_std_error, col=2, lty=2);
708 lines(days, group2-group2_std_error, col=2, lty=2);
709 # Standard error for group3.
710 lines(days, group3+group3_std_error, col=4, lty=2);
711 lines(days, group3-group3_std_error, col=4, lty=2);
712 }
713 }
714 }
715
716 stop_err = function(msg) {
717 cat(msg, file=stderr());
718 quit(save="no", status=1);
719 }
720
721 validate_date = function(date_str) {
722 valid_date = as.Date(date_str, format="%Y-%m-%d");
723 if( class(valid_date)=="try-error" || is.na(valid_date)) {
724 msg = paste("Invalid date: ", date_str, ", valid date format is yyyy-mm-dd.", sep="");
725 stop_err(msg);
726 }
727 return(valid_date);
728 }
729 500
730 if (is.null(opt$input_ytd)) { 501 if (is.null(opt$input_ytd)) {
731 processing_year_to_date_data = FALSE; 502 processing_year_to_date_data = FALSE;
732 } else { 503 } else {
733 processing_year_to_date_data = TRUE; 504 processing_year_to_date_data = TRUE;
748 prepend_end_doy_norm = data_list[[4]]; 519 prepend_end_doy_norm = data_list[[4]];
749 append_start_doy_norm = data_list[[5]]; 520 append_start_doy_norm = data_list[[5]];
750 is_leap_year = data_list[[6]]; 521 is_leap_year = data_list[[6]];
751 location = data_list[[7]]; 522 location = data_list[[7]];
752 523
753 if (is.null(opt$start_date) && is.null(opt$end_date)) { 524 # We're plotting an entire year.
754 # We're plotting an entire year. 525 # Display the total number of days in the Galaxy history item blurb.
755 date_interval = FALSE; 526 if (processing_year_to_date_data) {
756 # Display the total number of days in the Galaxy history item blurb. 527 cat("Number of days year-to-date: ", opt$num_days_ytd, "\n");
757 if (processing_year_to_date_data) { 528 } else {
758 cat("Number of days year-to-date: ", opt$num_days_ytd, "\n"); 529 if (is_leap_year) {
530 num_days = 366;
759 } else { 531 } else {
760 if (is_leap_year) { 532 num_days = 365;
761 num_days = 366; 533 }
762 } else { 534 cat("Number of days in year: ", num_days, "\n");
763 num_days = 365; 535 }
764 } 536
765 cat("Number of days in year: ", num_days, "\n");
766 }
767 } else {
768 # FIXME: currently custom date fields are free text, but
769 # Galaxy should soon include support for a date selector
770 # at which point this tool should be enhanced to use it.
771 # Validate start_date.
772 date_interval = TRUE;
773 # Calaculate the number of days in the date interval rather
774 # than using the number of rows in the input temperature data.
775 start_date = validate_date(opt$start_date);
776 # Validate end_date.
777 end_date = validate_date(opt$end_date);
778 if (start_date >= end_date) {
779 stop_err("The start date must be between 1 and 50 days before the end date when setting date intervals for plots.");
780 }
781 # Calculate the number of days in the date interval.
782 num_days = difftime(start_date, end_date, units=c("days"));
783 if (num_days > 50) {
784 # We need to restrict date intervals since
785 # plots render tick marks for each day.
786 stop_err("Date intervals for plotting cannot exceed 50 days.");
787 }
788 # Display the total number of days in the Galaxy history item blurb.
789 cat("Number of days in date interval: ", num_days, "\n");
790 }
791 # Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately. 537 # Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately.
792 if (plot_generations_separately) { 538 if (plot_generations_separately) {
793 temperature_data_frame_P = data.frame(temperature_data_frame); 539 temperature_data_frame_P = data.frame(temperature_data_frame);
794 temperature_data_frame_F1 = data.frame(temperature_data_frame); 540 temperature_data_frame_F1 = data.frame(temperature_data_frame);
795 temperature_data_frame_F2 = data.frame(temperature_data_frame); 541 temperature_data_frame_F2 = data.frame(temperature_data_frame);
796 } 542 }
797 543
798 # Get the ticks date labels for plots. 544 # Get the ticks date labels for plots.
799 ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, prepend_end_doy_norm, append_start_doy_norm, date_interval); 545 ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, prepend_end_doy_norm, append_start_doy_norm);
800 ticks = c(unlist(ticks_and_labels[1])); 546 ticks = c(unlist(ticks_and_labels[1]));
801 date_labels = c(unlist(ticks_and_labels[2])); 547 date_labels = c(unlist(ticks_and_labels[2]));
802 # All latitude values are the same, so get the value for plots from the first row. 548 # All latitude values are the same, so get the value for plots from the first row.
803 latitude = temperature_data_frame$LATITUDE[1]; 549 latitude = temperature_data_frame$LATITUDE[1];
804 550
1531 # Calculate adult populations for selected life stage. 1277 # Calculate adult populations for selected life stage.
1532 for (life_stage_adult in life_stages_adult) { 1278 for (life_stage_adult in life_stages_adult) {
1533 if (life_stage_adult == "Pre-vittelogenic") { 1279 if (life_stage_adult == "Pre-vittelogenic") {
1534 # Mean value for previttelogenic adults. 1280 # Mean value for previttelogenic adults.
1535 previttelogenic_adults = apply(Previttelogenic.replications, 1, mean); 1281 previttelogenic_adults = apply(Previttelogenic.replications, 1, mean);
1536 temperature_data_frame = append_vector(temperature_data_frame, previttelogenic_adults, "PRE-VITADULT"); 1282 temperature_data_frame = append_vector(temperature_data_frame, previttelogenic_adults, "PRE.VITADULT");
1537 # Standard error for previttelogenic adults. 1283 # Standard error for previttelogenic adults.
1538 previttelogenic_adults.std_error = apply(Previttelogenic.replications, 1, sd) / sqrt(opt$replications); 1284 previttelogenic_adults.std_error = apply(Previttelogenic.replications, 1, sd) / sqrt(opt$replications);
1539 temperature_data_frame = append_vector(temperature_data_frame, previttelogenic_adults.std_error, "PRE-VITADULTSE"); 1285 temperature_data_frame = append_vector(temperature_data_frame, previttelogenic_adults.std_error, "PRE.VITADULTSE");
1540 } else if (life_stage_adult == "Vittelogenic") { 1286 } else if (life_stage_adult == "Vittelogenic") {
1541 # Mean value for vittelogenic adults. 1287 # Mean value for vittelogenic adults.
1542 vittelogenic_adults = apply(Vittelogenic.replications, 1, mean); 1288 vittelogenic_adults = apply(Vittelogenic.replications, 1, mean);
1543 temperature_data_frame = append_vector(temperature_data_frame, vittelogenic_adults, "VITADULT"); 1289 temperature_data_frame = append_vector(temperature_data_frame, vittelogenic_adults, "VITADULT");
1544 # Standard error for vittelogenic adults. 1290 # Standard error for vittelogenic adults.
1572 F2.std_error = m_se[[6]]; 1318 F2.std_error = m_se[[6]];
1573 if (process_eggs) { 1319 if (process_eggs) {
1574 m_se = get_mean_and_std_error(P_eggs.replications, F1_eggs.replications, F2_eggs.replications); 1320 m_se = get_mean_and_std_error(P_eggs.replications, F1_eggs.replications, F2_eggs.replications);
1575 P_eggs = m_se[[1]]; 1321 P_eggs = m_se[[1]];
1576 P_eggs.std_error = m_se[[2]]; 1322 P_eggs.std_error = m_se[[2]];
1577 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_eggs, "EGG-P"); 1323 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_eggs, "EGG.P");
1578 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_eggs.std_error, "EGG-P-SE"); 1324 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_eggs.std_error, "EGG.P.SE");
1579 F1_eggs = m_se[[3]]; 1325 F1_eggs = m_se[[3]];
1580 F1_eggs.std_error = m_se[[4]]; 1326 F1_eggs.std_error = m_se[[4]];
1581 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_eggs, "EGG-F1"); 1327 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_eggs, "EGG.F1");
1582 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_eggs.std_error, "EGG-F1-SE"); 1328 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_eggs.std_error, "EGG.F1.SE");
1583 F2_eggs = m_se[[5]]; 1329 F2_eggs = m_se[[5]];
1584 F2_eggs.std_error = m_se[[6]]; 1330 F2_eggs.std_error = m_se[[6]];
1585 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_eggs, "EGG-F2"); 1331 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_eggs, "EGG.F2");
1586 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_eggs.std_error, "EGG-F2-SE"); 1332 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_eggs.std_error, "EGG.F2.SE");
1587 } 1333 }
1588 if (process_young_nymphs) { 1334 if (process_young_nymphs) {
1589 m_se = get_mean_and_std_error(P_young_nymphs.replications, F1_young_nymphs.replications, F2_young_nymphs.replications); 1335 m_se = get_mean_and_std_error(P_young_nymphs.replications, F1_young_nymphs.replications, F2_young_nymphs.replications);
1590 P_young_nymphs = m_se[[1]]; 1336 P_young_nymphs = m_se[[1]];
1591 P_young_nymphs.std_error = m_se[[2]]; 1337 P_young_nymphs.std_error = m_se[[2]];
1592 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_young_nymphs, "YOUNGNYMPH-P"); 1338 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_young_nymphs, "YOUNGNYMPH.P");
1593 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_young_nymphs.std_error, "YOUNGNYMPH-P-SE"); 1339 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_young_nymphs.std_error, "YOUNGNYMPH.P.SE");
1594 F1_young_nymphs = m_se[[3]]; 1340 F1_young_nymphs = m_se[[3]];
1595 F1_young_nymphs.std_error = m_se[[4]]; 1341 F1_young_nymphs.std_error = m_se[[4]];
1596 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_young_nymphs, "YOUNGNYMPH-F1"); 1342 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_young_nymphs, "YOUNGNYMPH.F1");
1597 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_young_nymphs.std_error, "YOUNGNYMPH-F1-SE"); 1343 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_young_nymphs.std_error, "YOUNGNYMPH.F1.SE");
1598 F2_young_nymphs = m_se[[5]]; 1344 F2_young_nymphs = m_se[[5]];
1599 F2_young_nymphs.std_error = m_se[[6]]; 1345 F2_young_nymphs.std_error = m_se[[6]];
1600 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_young_nymphs, "YOUNGNYMPH-F2"); 1346 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_young_nymphs, "YOUNGNYMPH.F2");
1601 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_young_nymphs.std_error, "YOUNGNYMPH-F2-SE"); 1347 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_young_nymphs.std_error, "YOUNGNYMPH.F2.SE");
1602 } 1348 }
1603 if (process_old_nymphs) { 1349 if (process_old_nymphs) {
1604 m_se = get_mean_and_std_error(P_old_nymphs.replications, F1_old_nymphs.replications, F2_old_nymphs.replications); 1350 m_se = get_mean_and_std_error(P_old_nymphs.replications, F1_old_nymphs.replications, F2_old_nymphs.replications);
1605 P_old_nymphs = m_se[[1]]; 1351 P_old_nymphs = m_se[[1]];
1606 P_old_nymphs.std_error = m_se[[2]]; 1352 P_old_nymphs.std_error = m_se[[2]];
1607 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_old_nymphs, "OLDNYMPH-P"); 1353 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_old_nymphs, "OLDNYMPH.P");
1608 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_old_nymphs.std_error, "OLDNYMPH-P-SE"); 1354 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_old_nymphs.std_error, "OLDNYMPH.P.SE");
1609 F1_old_nymphs = m_se[[3]]; 1355 F1_old_nymphs = m_se[[3]];
1610 F1_old_nymphs.std_error = m_se[[4]]; 1356 F1_old_nymphs.std_error = m_se[[4]];
1611 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_old_nymphs, "OLDNYMPH-F1"); 1357 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_old_nymphs, "OLDNYMPH.F1");
1612 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_old_nymphs.std_error, "OLDNYMPH-F1-SE"); 1358 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_old_nymphs.std_error, "OLDNYMPH.F1.SE");
1613 F2_old_nymphs = m_se[[5]]; 1359 F2_old_nymphs = m_se[[5]];
1614 F2_old_nymphs.std_error = m_se[[6]]; 1360 F2_old_nymphs.std_error = m_se[[6]];
1615 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_old_nymphs, "OLDNYMPH-F2"); 1361 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_old_nymphs, "OLDNYMPH.F2");
1616 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_old_nymphs.std_error, "OLDNYMPH-F2-SE"); 1362 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_old_nymphs.std_error, "OLDNYMPH.F2.SE");
1617 } 1363 }
1618 if (process_total_nymphs) { 1364 if (process_total_nymphs) {
1619 m_se = get_mean_and_std_error(P_total_nymphs.replications, F1_total_nymphs.replications, F2_total_nymphs.replications); 1365 m_se = get_mean_and_std_error(P_total_nymphs.replications, F1_total_nymphs.replications, F2_total_nymphs.replications);
1620 P_total_nymphs = m_se[[1]]; 1366 P_total_nymphs = m_se[[1]];
1621 P_total_nymphs.std_error = m_se[[2]]; 1367 P_total_nymphs.std_error = m_se[[2]];
1622 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_nymphs, "TOTALNYMPH-P"); 1368 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_nymphs, "TOTALNYMPH.P");
1623 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_nymphs.std_error, "TOTALNYMPH-P-SE"); 1369 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_nymphs.std_error, "TOTALNYMPH.P.SE");
1624 F1_total_nymphs = m_se[[3]]; 1370 F1_total_nymphs = m_se[[3]];
1625 F1_total_nymphs.std_error = m_se[[4]]; 1371 F1_total_nymphs.std_error = m_se[[4]];
1626 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_nymphs, "TOTALNYMPH-F1"); 1372 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_nymphs, "TOTALNYMPH.F1");
1627 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_nymphs.std_error, "TOTALNYMPH-F1-SE"); 1373 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_nymphs.std_error, "TOTALNYMPH.F1.SE");
1628 F2_total_nymphs = m_se[[5]]; 1374 F2_total_nymphs = m_se[[5]];
1629 F2_total_nymphs.std_error = m_se[[6]]; 1375 F2_total_nymphs.std_error = m_se[[6]];
1630 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_nymphs, "TOTALNYMPH-F2"); 1376 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_nymphs, "TOTALNYMPH.F2");
1631 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_nymphs.std_error, "TOTALNYMPH-F2-SE"); 1377 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_nymphs.std_error, "TOTALNYMPH.F2.SE");
1632 } 1378 }
1633 if (process_previttelogenic_adults) { 1379 if (process_previttelogenic_adults) {
1634 m_se = get_mean_and_std_error(P_previttelogenic_adults.replications, F1_previttelogenic_adults.replications, F2_previttelogenic_adults.replications); 1380 m_se = get_mean_and_std_error(P_previttelogenic_adults.replications, F1_previttelogenic_adults.replications, F2_previttelogenic_adults.replications);
1635 P_previttelogenic_adults = m_se[[1]]; 1381 P_previttelogenic_adults = m_se[[1]];
1636 P_previttelogenic_adults.std_error = m_se[[2]]; 1382 P_previttelogenic_adults.std_error = m_se[[2]];
1637 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_previttelogenic_adults, "PRE-VITADULT-P"); 1383 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_previttelogenic_adults, "PRE.VITADULT.P");
1638 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_previttelogenic_adults.std_error, "PRE-VITADULT-P-SE"); 1384 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_previttelogenic_adults.std_error, "PRE.VITADULT.P.SE");
1639 F1_previttelogenic_adults = m_se[[3]]; 1385 F1_previttelogenic_adults = m_se[[3]];
1640 F1_previttelogenic_adults.std_error = m_se[[4]]; 1386 F1_previttelogenic_adults.std_error = m_se[[4]];
1641 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_previttelogenic_adults, "PRE-VITADULT-F1"); 1387 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_previttelogenic_adults, "PRE.VITADULT.F1");
1642 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_previttelogenic_adults.std_error, "PRE-VITADULT-F1-SE"); 1388 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_previttelogenic_adults.std_error, "PRE.VITADULT.F1.SE");
1643 F2_previttelogenic_adults = m_se[[5]]; 1389 F2_previttelogenic_adults = m_se[[5]];
1644 F2_previttelogenic_adults.std_error = m_se[[6]]; 1390 F2_previttelogenic_adults.std_error = m_se[[6]];
1645 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_previttelogenic_adults, "PRE-VITADULT-F2"); 1391 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_previttelogenic_adults, "PRE.VITADULT.F2");
1646 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_previttelogenic_adults.std_error, "PRE-VITADULT-F2-SE"); 1392 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_previttelogenic_adults.std_error, "PRE.VITADULT.F2.SE");
1647 } 1393 }
1648 if (process_vittelogenic_adults) { 1394 if (process_vittelogenic_adults) {
1649 m_se = get_mean_and_std_error(P_vittelogenic_adults.replications, F1_vittelogenic_adults.replications, F2_vittelogenic_adults.replications); 1395 m_se = get_mean_and_std_error(P_vittelogenic_adults.replications, F1_vittelogenic_adults.replications, F2_vittelogenic_adults.replications);
1650 P_vittelogenic_adults = m_se[[1]]; 1396 P_vittelogenic_adults = m_se[[1]];
1651 P_vittelogenic_adults.std_error = m_se[[2]]; 1397 P_vittelogenic_adults.std_error = m_se[[2]];
1652 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_vittelogenic_adults, "VITADULT-P"); 1398 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_vittelogenic_adults, "VITADULT.P");
1653 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_vittelogenic_adults.std_error, "VITADULT-P-SE"); 1399 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_vittelogenic_adults.std_error, "VITADULT.P.SE");
1654 F1_vittelogenic_adults = m_se[[3]]; 1400 F1_vittelogenic_adults = m_se[[3]];
1655 F1_vittelogenic_adults.std_error = m_se[[4]]; 1401 F1_vittelogenic_adults.std_error = m_se[[4]];
1656 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_vittelogenic_adults, "VITADULT-F1"); 1402 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_vittelogenic_adults, "VITADULT.F1");
1657 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_vittelogenic_adults.std_error, "VITADULT-F1-SE"); 1403 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_vittelogenic_adults.std_error, "VITADULT.F1.SE");
1658 F2_vittelogenic_adults = m_se[[5]]; 1404 F2_vittelogenic_adults = m_se[[5]];
1659 F2_vittelogenic_adults.std_error = m_se[[6]]; 1405 F2_vittelogenic_adults.std_error = m_se[[6]];
1660 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_vittelogenic_adults, "VITADULT-F2"); 1406 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_vittelogenic_adults, "VITADULT.F2");
1661 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_vittelogenic_adults.std_error, "VITADULT-F2-SE"); 1407 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_vittelogenic_adults.std_error, "VITADULT.F2.SE");
1662 } 1408 }
1663 if (process_diapausing_adults) { 1409 if (process_diapausing_adults) {
1664 m_se = get_mean_and_std_error(P_diapausing_adults.replications, F1_diapausing_adults.replications, F2_diapausing_adults.replications); 1410 m_se = get_mean_and_std_error(P_diapausing_adults.replications, F1_diapausing_adults.replications, F2_diapausing_adults.replications);
1665 P_diapausing_adults = m_se[[1]]; 1411 P_diapausing_adults = m_se[[1]];
1666 P_diapausing_adults.std_error = m_se[[2]]; 1412 P_diapausing_adults.std_error = m_se[[2]];
1667 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_diapausing_adults, "DIAPAUSINGADULT-P"); 1413 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_diapausing_adults, "DIAPAUSINGADULT.P");
1668 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_diapausing_adults.std_error, "DIAPAUSINGADULT-P-SE"); 1414 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_diapausing_adults.std_error, "DIAPAUSINGADULT.P.SE");
1669 F1_diapausing_adults = m_se[[3]]; 1415 F1_diapausing_adults = m_se[[3]];
1670 F1_diapausing_adults.std_error = m_se[[4]]; 1416 F1_diapausing_adults.std_error = m_se[[4]];
1671 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_diapausing_adults, "DIAPAUSINGADULT-F1"); 1417 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_diapausing_adults, "DIAPAUSINGADULT.F1");
1672 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_diapausing_adults.std_error, "DIAPAUSINGADULT-F1-SE"); 1418 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_diapausing_adults.std_error, "DIAPAUSINGADULT.F1.SE");
1673 F2_diapausing_adults = m_se[[5]]; 1419 F2_diapausing_adults = m_se[[5]];
1674 F2_diapausing_adults.std_error = m_se[[6]]; 1420 F2_diapausing_adults.std_error = m_se[[6]];
1675 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_diapausing_adults, "DIAPAUSINGADULT-F2"); 1421 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_diapausing_adults, "DIAPAUSINGADULT.F2");
1676 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_diapausing_adults.std_error, "DIAPAUSINGADULT-F2-SE"); 1422 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_diapausing_adults.std_error, "DIAPAUSINGADULT.F2.SE");
1677 } 1423 }
1678 if (process_total_adults) { 1424 if (process_total_adults) {
1679 m_se = get_mean_and_std_error(P_total_adults.replications, F1_total_adults.replications, F2_total_adults.replications); 1425 m_se = get_mean_and_std_error(P_total_adults.replications, F1_total_adults.replications, F2_total_adults.replications);
1680 P_total_adults = m_se[[1]]; 1426 P_total_adults = m_se[[1]];
1681 P_total_adults.std_error = m_se[[2]]; 1427 P_total_adults.std_error = m_se[[2]];
1682 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_adults, "TOTALADULT-P"); 1428 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_adults, "TOTALADULT.P");
1683 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_adults.std_error, "TOTALADULT-P-SE"); 1429 temperature_data_frame_P = append_vector(temperature_data_frame_P, P_total_adults.std_error, "TOTALADULT.P.SE");
1684 F1_total_adults = m_se[[3]]; 1430 F1_total_adults = m_se[[3]];
1685 F1_total_adults.std_error = m_se[[4]]; 1431 F1_total_adults.std_error = m_se[[4]];
1686 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_adults, "TOTALADULT-F1"); 1432 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_adults, "TOTALADULT.F1");
1687 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_adults.std_error, "TOTALADULT-F1-SE"); 1433 temperature_data_frame_F1 = append_vector(temperature_data_frame_F1, F1_total_adults.std_error, "TOTALADULT.F1.SE");
1688 F2_total_adults = m_se[[5]]; 1434 F2_total_adults = m_se[[5]];
1689 F2_total_adults.std_error = m_se[[6]]; 1435 F2_total_adults.std_error = m_se[[6]];
1690 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_adults, "TOTALADULT-F2"); 1436 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_adults, "TOTALADULT.F2");
1691 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_adults.std_error, "TOTALADULT-F2-SE"); 1437 temperature_data_frame_F2 = append_vector(temperature_data_frame_F2, F2_total_adults.std_error, "TOTALADULT.F2.SE");
1692 } 1438 }
1693 } 1439 }
1694 1440
1695 # Save the analyzed data for combined generations. 1441 # Save the analyzed data for combined generations.
1696 file_path = paste("output_data_dir", "04_combined_generations.csv", sep="/"); 1442 file_path = paste("output_data_dir", "04_combined_generations.csv", sep="/");
1725 dev.off(); 1471 dev.off();
1726 } else if (life_stage == "Nymph") { 1472 } else if (life_stage == "Nymph") {
1727 for (life_stage_nymph in life_stages_nymph) { 1473 for (life_stage_nymph in life_stages_nymph) {
1728 # Start PDF device driver. 1474 # Start PDF device driver.
1729 dev.new(width=20, height=30); 1475 dev.new(width=20, height=30);
1730 file_path = get_file_path(life_stage, "nymph_pop_by_generation.pdf", life_stage_nymph=life_stage_nymph) 1476 file_path = get_file_path(life_stage, "nymph_pop_by_generation.pdf", sub_life_stage=life_stage_nymph)
1731 pdf(file=file_path, width=20, height=30, bg="white"); 1477 pdf(file=file_path, width=20, height=30, bg="white");
1732 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1478 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1733 if (life_stage_nymph=="Young") { 1479 if (life_stage_nymph=="Young") {
1734 # Young nymph population size by generation. 1480 # Young nymph population size by generation.
1735 maxval = max(P_young_nymphs+F1_young_nymphs+F2_young_nymphs) + 100; 1481 maxval = max(P_young_nymphs+F1_young_nymphs+F2_young_nymphs) + 100;
1758 group3 = F2_total_nymphs; 1504 group3 = F2_total_nymphs;
1759 group3_std_error = F2_total_nymphs.std_error; 1505 group3_std_error = F2_total_nymphs.std_error;
1760 } 1506 }
1761 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude, 1507 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
1762 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error, 1508 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1763 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_nymph=life_stage_nymph); 1509 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, sub_life_stage=life_stage_nymph);
1764 # Turn off device driver to flush output. 1510 # Turn off device driver to flush output.
1765 dev.off(); 1511 dev.off();
1766 } 1512 }
1767 } else if (life_stage == "Adult") { 1513 } else if (life_stage == "Adult") {
1768 for (life_stage_adult in life_stages_adult) { 1514 for (life_stage_adult in life_stages_adult) {
1769 # Start PDF device driver. 1515 # Start PDF device driver.
1770 dev.new(width=20, height=30); 1516 dev.new(width=20, height=30);
1771 file_path = get_file_path(life_stage, "adult_pop_by_generation.pdf", life_stage_adult=life_stage_adult) 1517 file_path = get_file_path(life_stage, "adult_pop_by_generation.pdf", sub_life_stage=life_stage_adult)
1772 pdf(file=file_path, width=20, height=30, bg="white"); 1518 pdf(file=file_path, width=20, height=30, bg="white");
1773 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1519 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1774 if (life_stage_adult=="Pre-vittelogenic") { 1520 if (life_stage_adult=="Pre-vittelogenic") {
1775 # Pre-vittelogenic adult population size by generation. 1521 # Pre-vittelogenic adult population size by generation.
1776 maxval = max(P_previttelogenic_adults+F1_previttelogenic_adults+F2_previttelogenic_adults) + 100; 1522 maxval = max(P_previttelogenic_adults+F1_previttelogenic_adults+F2_previttelogenic_adults) + 100;
1808 group3 = F2_total_adults; 1554 group3 = F2_total_adults;
1809 group3_std_error = F2_total_adults.std_error; 1555 group3_std_error = F2_total_adults.std_error;
1810 } 1556 }
1811 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude, 1557 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
1812 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error, 1558 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1813 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_adult=life_stage_adult); 1559 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, sub_life_stage=life_stage_adult);
1814 # Turn off device driver to flush output. 1560 # Turn off device driver to flush output.
1815 dev.off(); 1561 dev.off();
1816 } 1562 }
1817 } else if (life_stage == "Total") { 1563 } else if (life_stage == "Total") {
1818 # Start PDF device driver. 1564 # Start PDF device driver.
1847 dev.off(); 1593 dev.off();
1848 } else if (life_stage == "Nymph") { 1594 } else if (life_stage == "Nymph") {
1849 for (life_stage_nymph in life_stages_nymph) { 1595 for (life_stage_nymph in life_stages_nymph) {
1850 # Start PDF device driver. 1596 # Start PDF device driver.
1851 dev.new(width=20, height=30); 1597 dev.new(width=20, height=30);
1852 file_path = get_file_path(life_stage, "nymph_pop.pdf", life_stage_nymph=life_stage_nymph) 1598 file_path = get_file_path(life_stage, "nymph_pop.pdf", sub_life_stage=life_stage_nymph)
1853 pdf(file=file_path, width=20, height=30, bg="white"); 1599 pdf(file=file_path, width=20, height=30, bg="white");
1854 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1600 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1855 if (life_stage_nymph=="Total") { 1601 if (life_stage_nymph=="Total") {
1856 # Total nymph population size. 1602 # Total nymph population size.
1857 group = total_nymphs; 1603 group = total_nymphs;
1866 group_std_error = old_nymphs.std_error; 1612 group_std_error = old_nymphs.std_error;
1867 } 1613 }
1868 maxval = max(group+group_std_error) + 100; 1614 maxval = max(group+group_std_error) + 100;
1869 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude, 1615 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
1870 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error, 1616 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1871 life_stages_nymph=life_stage_nymph); 1617 sub_life_stage=life_stage_nymph);
1872 # Turn off device driver to flush output. 1618 # Turn off device driver to flush output.
1873 dev.off(); 1619 dev.off();
1874 } 1620 }
1875 } else if (life_stage == "Adult") { 1621 } else if (life_stage == "Adult") {
1876 for (life_stage_adult in life_stages_adult) { 1622 for (life_stage_adult in life_stages_adult) {
1877 # Start PDF device driver. 1623 # Start PDF device driver.
1878 dev.new(width=20, height=30); 1624 dev.new(width=20, height=30);
1879 file_path = get_file_path(life_stage, "adult_pop.pdf", life_stage_adult=life_stage_adult) 1625 file_path = get_file_path(life_stage, "adult_pop.pdf", sub_life_stage=life_stage_adult)
1880 pdf(file=file_path, width=20, height=30, bg="white"); 1626 pdf(file=file_path, width=20, height=30, bg="white");
1881 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1627 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1882 if (life_stage_adult=="Total") { 1628 if (life_stage_adult=="Total") {
1883 # Total adult population size. 1629 # Total adult population size.
1884 group = total_adults; 1630 group = total_adults;
1897 group_std_error = diapausing_adults.std_error 1643 group_std_error = diapausing_adults.std_error
1898 } 1644 }
1899 maxval = max(group+group_std_error) + 100; 1645 maxval = max(group+group_std_error) + 100;
1900 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude, 1646 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
1901 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error, 1647 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1902 life_stages_adult=life_stage_adult); 1648 sub_life_stage=life_stage_adult);
1903 # Turn off device driver to flush output. 1649 # Turn off device driver to flush output.
1904 dev.off(); 1650 dev.off();
1905 } 1651 }
1906 } else if (life_stage == "Total") { 1652 } else if (life_stage == "Total") {
1907 # Start PDF device driver. 1653 # Start PDF device driver.