Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 50:927321ed0322 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:59:06 -0400
parents 1b6864c5b50a
children 3f827fe31756
comparison
equal deleted inserted replaced
49:1b6864c5b50a 50:927321ed0322
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) {
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(norm_data_frame, start_date_doy, end_date_doy, year) { 68 from_30_year_normals = function(norm_data_frame, start_date_doy, end_date_doy, year) {
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 # Add 1 to the number of rows to ensure that the end date is included. 72 # Add 1 to the number of rows to ensure that the end date is included.
82 for (i in first_norm_row:last_norm_row) { 76 for (i in first_norm_row:last_norm_row) {
83 j = j + 1; 77 j = j + 1;
84 tmp_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);
85 } 79 }
86 return (tmp_data_frame); 80 return (tmp_data_frame);
87 }
88
89 get_file_path = function(life_stage, base_name, life_stage_nymph=NULL, life_stage_adult=NULL) {
90 if (!is.null(life_stage_nymph)) {
91 lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph);
92 file_name = paste(lsi, tolower(life_stage_nymph), base_name, sep="_");
93 } else if (!is.null(life_stage_adult)) {
94 lsi = get_life_stage_index(life_stage, life_stage_adult=life_stage_adult);
95 file_name = paste(lsi, tolower(life_stage_adult), base_name, sep="_");
96 } else {
97 lsi = get_life_stage_index(life_stage);
98 file_name = paste(lsi, base_name, sep="_");
99 }
100 file_path = paste("output_plots_dir", file_name, sep="/");
101 return(file_path);
102 }
103
104 get_life_stage_index = function(life_stage, life_stage_nymph=NULL, life_stage_adult=NULL) {
105 # Name collection elements so that they
106 # are displayed in logical order.
107 if (life_stage=="Egg") {
108 lsi = "01";
109 } else if (life_stage=="Nymph") {
110 if (life_stage_nymph=="Young") {
111 lsi = "02";
112 } else if (life_stage_nymph=="Old") {
113 lsi = "03";
114 } else if (life_stage_nymph=="Total") {
115 lsi="04";
116 }
117 } else if (life_stage=="Adult") {
118 if (life_stage_adult=="Pre-vittelogenic") {
119 lsi = "05";
120 } else if (life_stage_adult=="Vittelogenic") {
121 lsi = "06";
122 } else if (life_stage_adult=="Diapausing") {
123 lsi = "07";
124 } else if (life_stage_adult=="Total") {
125 lsi = "08";
126 }
127 } else if (life_stage=="Total") {
128 lsi = "09";
129 }
130 return(lsi);
131 }
132
133 get_mean_and_std_error = function(p_replications, f1_replications, f2_replications) {
134 # P mean.
135 p_m = apply(p_replications, 1, mean);
136 # P standard error.
137 p_se = apply(p_replications, 1, sd) / sqrt(opt$replications);
138 # F1 mean.
139 f1_m = apply(f1_replications, 1, mean);
140 # F1 standard error.
141 f1_se = apply(f1_replications, 1, sd) / sqrt(opt$replications);
142 # F2 mean.
143 f2_m = apply(f2_replications, 1, mean);
144 # F2 standard error.
145 f2_se = apply(f2_replications, 1, sd) / sqrt(opt$replications);
146 return(list(p_m, p_se, f1_m, f1_se, f2_m, f2_se))
147 } 81 }
148 82
149 get_new_norm_data_frame = function(is_leap_year, input_norm=NULL, nrow=0) { 83 get_new_norm_data_frame = function(is_leap_year, input_norm=NULL, nrow=0) {
150 # The input_norm data has the following 10 columns: 84 # The input_norm data has the following 10 columns:
151 # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX 85 # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX
269 } 203 }
270 } 204 }
271 averages = sum(dh) / 24; 205 averages = sum(dh) / 24;
272 } 206 }
273 return(c(curr_mean_temp, averages)) 207 return(c(curr_mean_temp, averages))
274 }
275
276 get_tick_index = function(index, last_tick, ticks, tick_labels, tick_sep) {
277 # The R code tries hard not to draw overlapping tick labels, and so
278 # will omit labels where they would abut or overlap previously drawn
279 # labels. This can result in, for example, every other tick being
280 # labelled. We'll keep track of the last tick to make sure all of
281 # the month labels are displayed, and missing ticks are restricted
282 # to Sundays which have no labels anyway.
283 if (last_tick==0) {
284 return(length(ticks)+1);
285 }
286 last_saved_tick = ticks[[length(ticks)]];
287 if (index-last_saved_tick<tick_sep) {
288 last_saved_month = tick_labels[[length(tick_labels)]];
289 if (last_saved_month=="") {
290 # We're safe overwriting a tick
291 # with no label (i.e., a Sunday tick).
292 return(length(ticks));
293 } else {
294 # Don't eliminate a Month label.
295 return(NULL);
296 }
297 }
298 return(length(ticks)+1);
299 }
300
301 get_total_days = function(is_leap_year) {
302 # Get the total number of days in the current year.
303 if (is_leap_year) {
304 return(366);
305 } else {
306 return(365);
307 }
308 }
309
310 get_x_axis_ticks_and_labels = function(temperature_data_frame, prepend_end_doy_norm, append_start_doy_norm, date_interval) {
311 # Generate a list of ticks and labels for plotting the x axis.
312 if (prepend_end_doy_norm > 0) {
313 prepend_end_norm_row = which(temperature_data_frame$DOY==prepend_end_doy_norm);
314 } else {
315 prepend_end_norm_row = 0;
316 }
317 if (append_start_doy_norm > 0) {
318 append_start_norm_row = which(temperature_data_frame$DOY==append_start_doy_norm);
319 } else {
320 append_start_norm_row = 0;
321 }
322 num_rows = dim(temperature_data_frame)[1];
323 tick_labels = list();
324 ticks = list();
325 current_month_label = NULL;
326 last_tick = 0;
327 if (date_interval) {
328 tick_sep = 0;
329 } else {
330 tick_sep = 3;
331 }
332 for (i in 1:num_rows) {
333 # Get the year and month from the date which
334 # has the format YYYY-MM-DD.
335 date = format(temperature_data_frame$DATE[i]);
336 # Get the month label.
337 items = strsplit(date, "-")[[1]];
338 month = items[2];
339 month_label = month.abb[as.integer(month)];
340 day = as.integer(items[3]);
341 doy = as.integer(temperature_data_frame$DOY[i]);
342 # We're plotting the entire year, so ticks will
343 # occur on Sundays and the first of each month.
344 if (i == prepend_end_norm_row) {
345 # Add a tick for the end of the 30 year normnals data
346 # that was prepended to the year-to-date data.
347 label_str = "End prepended 30 year normals";
348 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
349 ticks[tick_index] = i;
350 if (date_interval) {
351 # Append the day to label_str
352 tick_labels[tick_index] = paste(label_str, day, sep=" ");
353 } else {
354 tick_labels[tick_index] = label_str;
355 }
356 last_tick = i;
357 } else if (doy == append_start_doy_norm) {
358 # Add a tick for the start of the 30 year normnals data
359 # that was appended to the year-to-date data.
360 label_str = "Start appended 30 year normals";
361 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
362 ticks[tick_index] = i;
363 if (!identical(current_month_label, month_label)) {
364 # Append the month to label_str.
365 label_str = paste(label_str, month_label, spe=" ");
366 current_month_label = month_label;
367 }
368 if (date_interval) {
369 # Append the day to label_str
370 label_str = paste(label_str, day, sep=" ");
371 }
372 tick_labels[tick_index] = label_str;
373 last_tick = i;
374 } else if (i==num_rows) {
375 # Add a tick for the last day of the year.
376 label_str = "";
377 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
378 ticks[tick_index] = i;
379 if (!identical(current_month_label, month_label)) {
380 # Append the month to label_str.
381 label_str = month_label;
382 current_month_label = month_label;
383 }
384 if (date_interval) {
385 # Append the day to label_str
386 label_str = paste(label_str, day, sep=" ");
387 }
388 tick_labels[tick_index] = label_str;
389 } else {
390 if (!identical(current_month_label, month_label)) {
391 # Add a tick for the month.
392 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
393 ticks[tick_index] = i;
394 if (date_interval) {
395 # Append the day to the month.
396 tick_labels[tick_index] = paste(month_label, day, sep=" ");
397 } else {
398 tick_labels[tick_index] = month_label;
399 }
400 current_month_label = month_label;
401 last_tick = i;
402 }
403 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
404 if (!is.null(tick_index)) {
405 if (date_interval) {
406 # Add a tick for every day. The first tick is the
407 # month label, so add a tick only if i is not 1
408 if (i>1 & day>1) {
409 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep)
410 ticks[tick_index] = i;
411 # Add the day as the label.
412 tick_labels[tick_index] = day;
413 last_tick = i;
414 }
415 } else {
416 # Get the day.
417 day = weekdays(as.Date(date));
418 if (day=="Sunday") {
419 # Add a tick if we're on a Sunday.
420 ticks[tick_index] = i;
421 # Add a blank month label so it is not displayed.
422 tick_labels[tick_index] = "";
423 last_tick = i;
424 }
425 }
426 }
427 }
428 }
429 return(list(ticks, tick_labels));
430 }
431
432 get_year_from_date = function(date_str) {
433 date_str_items = strsplit(date_str, "-")[[1]];
434 return (date_str_items[1]);
435 } 208 }
436 209
437 is_leap_year = function(date_str) { 210 is_leap_year = function(date_str) {
438 # Extract the year from the date_str. 211 # Extract the year from the date_str.
439 date = format(date_str); 212 date = format(date_str);
719 # Add a column containing the daylight length for each day. 492 # Add a column containing the daylight length for each day.
720 temperature_data_frame = add_daylight_length(temperature_data_frame); 493 temperature_data_frame = add_daylight_length(temperature_data_frame);
721 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));
722 } 495 }
723 496
724 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.
725 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="/");
726 life_stages_adult=NULL, life_stages_nymph=NULL) { 499 source(utils_path);
727 if (chart_type=="pop_size_by_life_stage") {
728 if (life_stage=="Total") {
729 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
730 legend_text = c("Egg", "Nymph", "Adult");
731 columns = c(4, 2, 1);
732 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);
733 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
734 lines(days, group2, lwd=2, lty=1, col=2);
735 lines(days, group3, lwd=2, lty=1, col=4);
736 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);
737 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
738 if (plot_std_error=="yes") {
739 # Standard error for group.
740 lines(days, group+group_std_error, lty=2);
741 lines(days, group-group_std_error, lty=2);
742 # Standard error for group2.
743 lines(days, group2+group2_std_error, col=2, lty=2);
744 lines(days, group2-group2_std_error, col=2, lty=2);
745 # Standard error for group3.
746 lines(days, group3+group3_std_error, col=4, lty=2);
747 lines(days, group3-group3_std_error, col=4, lty=2);
748 }
749 } else {
750 if (life_stage=="Egg") {
751 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
752 legend_text = c(life_stage);
753 columns = c(4);
754 } else if (life_stage=="Nymph") {
755 stage = paste(life_stages_nymph, "Nymph Pop :", sep=" ");
756 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
757 legend_text = c(paste(life_stages_nymph, life_stage, sep=" "));
758 columns = c(2);
759 } else if (life_stage=="Adult") {
760 stage = paste(life_stages_adult, "Adult Pop", sep=" ");
761 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
762 legend_text = c(paste(life_stages_adult, life_stage, sep=" "));
763 columns = c(1);
764 }
765 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);
766 legend("topleft", legend_text, lty=c(1), col="black", cex=3);
767 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);
768 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
769 if (plot_std_error=="yes") {
770 # Standard error for group.
771 lines(days, group+group_std_error, lty=2);
772 lines(days, group-group_std_error, lty=2);
773 }
774 }
775 } else if (chart_type=="pop_size_by_generation") {
776 if (life_stage=="Total") {
777 title_str = ": Total Pop by Gen :";
778 } else if (life_stage=="Egg") {
779 title_str = ": Egg Pop by Gen :";
780 } else if (life_stage=="Nymph") {
781 title_str = paste(":", life_stages_nymph, "Nymph Pop by Gen", ":", sep=" ");
782 } else if (life_stage=="Adult") {
783 title_str = paste(":", life_stages_adult, "Adult Pop by Gen", ":", sep=" ");
784 }
785 title = paste(insect, ": Reps", replications, title_str, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
786 legend_text = c("P", "F1", "F2");
787 columns = c(1, 2, 4);
788 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);
789 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3);
790 lines(days, group2, lwd=2, lty=1, col=2);
791 lines(days, group3, lwd=2, lty=1, col=4);
792 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);
793 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3);
794 if (plot_std_error=="yes") {
795 # Standard error for group.
796 lines(days, group+group_std_error, lty=2);
797 lines(days, group-group_std_error, lty=2);
798 # Standard error for group2.
799 lines(days, group2+group2_std_error, col=2, lty=2);
800 lines(days, group2-group2_std_error, col=2, lty=2);
801 # Standard error for group3.
802 lines(days, group3+group3_std_error, col=4, lty=2);
803 lines(days, group3-group3_std_error, col=4, lty=2);
804 }
805 }
806 }
807
808 stop_err = function(msg) {
809 cat(msg, file=stderr());
810 quit(save="no", status=1);
811 }
812
813 validate_date = function(date_str) {
814 valid_date = as.Date(date_str, format="%Y-%m-%d");
815 if( class(valid_date)=="try-error" || is.na(valid_date)) {
816 msg = paste("Invalid date: ", date_str, ", valid date format is yyyy-mm-dd.", sep="");
817 stop_err(msg);
818 }
819 return(valid_date);
820 }
821 500
822 if (is.null(opt$input_ytd)) { 501 if (is.null(opt$input_ytd)) {
823 processing_year_to_date_data = FALSE; 502 processing_year_to_date_data = FALSE;
824 } else { 503 } else {
825 processing_year_to_date_data = TRUE; 504 processing_year_to_date_data = TRUE;
840 prepend_end_doy_norm = data_list[[4]]; 519 prepend_end_doy_norm = data_list[[4]];
841 append_start_doy_norm = data_list[[5]]; 520 append_start_doy_norm = data_list[[5]];
842 is_leap_year = data_list[[6]]; 521 is_leap_year = data_list[[6]];
843 location = data_list[[7]]; 522 location = data_list[[7]];
844 523
845 if (is.null(opt$start_date) && is.null(opt$end_date)) { 524 # We're plotting an entire year.
846 # We're plotting an entire year. 525 # Display the total number of days in the Galaxy history item blurb.
847 date_interval = FALSE; 526 if (processing_year_to_date_data) {
848 # 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");
849 if (processing_year_to_date_data) { 528 } else {
850 cat("Number of days year-to-date: ", opt$num_days_ytd, "\n"); 529 if (is_leap_year) {
530 num_days = 366;
851 } else { 531 } else {
852 if (is_leap_year) { 532 num_days = 365;
853 num_days = 366; 533 }
854 } else { 534 cat("Number of days in year: ", num_days, "\n");
855 num_days = 365; 535 }
856 } 536
857 cat("Number of days in year: ", num_days, "\n");
858 }
859 } else {
860 # FIXME: currently custom date fields are free text, but
861 # Galaxy should soon include support for a date selector
862 # at which point this tool should be enhanced to use it.
863 # Validate start_date.
864 date_interval = TRUE;
865 # Calaculate the number of days in the date interval rather
866 # than using the number of rows in the input temperature data.
867 start_date = validate_date(opt$start_date);
868 # Validate end_date.
869 end_date = validate_date(opt$end_date);
870 if (start_date >= end_date) {
871 stop_err("The start date must be between 1 and 50 days before the end date when setting date intervals for plots.");
872 }
873 # Calculate the number of days in the date interval.
874 num_days = difftime(end_date, start_date, units=c("days"));
875 # Add 1 to the number of days to make the dates inclusive. For
876 # example, if the user enters a date range of 2018-01-01 to
877 # 2018-01-31, they likely expect the end date to be included.
878 num_days = num_days + 1;
879 if (num_days > 50) {
880 # We need to restrict date intervals since
881 # plots render tick marks for each day.
882 stop_err("Date intervals for plotting cannot exceed 50 days.");
883 }
884 # Display the total number of days in the Galaxy history item blurb.
885 cat("Number of days in date interval: ", num_days, "\n");
886 }
887 # 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.
888 if (plot_generations_separately) { 538 if (plot_generations_separately) {
889 temperature_data_frame_P = data.frame(temperature_data_frame); 539 temperature_data_frame_P = data.frame(temperature_data_frame);
890 temperature_data_frame_F1 = data.frame(temperature_data_frame); 540 temperature_data_frame_F1 = data.frame(temperature_data_frame);
891 temperature_data_frame_F2 = data.frame(temperature_data_frame); 541 temperature_data_frame_F2 = data.frame(temperature_data_frame);
892 } 542 }
893 543
894 # Get the ticks date labels for plots. 544 # Get the ticks date labels for plots.
895 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);
896 ticks = c(unlist(ticks_and_labels[1])); 546 ticks = c(unlist(ticks_and_labels[1]));
897 date_labels = c(unlist(ticks_and_labels[2])); 547 date_labels = c(unlist(ticks_and_labels[2]));
898 # 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.
899 latitude = temperature_data_frame$LATITUDE[1]; 549 latitude = temperature_data_frame$LATITUDE[1];
900 550
1627 # Calculate adult populations for selected life stage. 1277 # Calculate adult populations for selected life stage.
1628 for (life_stage_adult in life_stages_adult) { 1278 for (life_stage_adult in life_stages_adult) {
1629 if (life_stage_adult == "Pre-vittelogenic") { 1279 if (life_stage_adult == "Pre-vittelogenic") {
1630 # Mean value for previttelogenic adults. 1280 # Mean value for previttelogenic adults.
1631 previttelogenic_adults = apply(Previttelogenic.replications, 1, mean); 1281 previttelogenic_adults = apply(Previttelogenic.replications, 1, mean);
1632 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");
1633 # Standard error for previttelogenic adults. 1283 # Standard error for previttelogenic adults.
1634 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);
1635 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");
1636 } else if (life_stage_adult == "Vittelogenic") { 1286 } else if (life_stage_adult == "Vittelogenic") {
1637 # Mean value for vittelogenic adults. 1287 # Mean value for vittelogenic adults.
1638 vittelogenic_adults = apply(Vittelogenic.replications, 1, mean); 1288 vittelogenic_adults = apply(Vittelogenic.replications, 1, mean);
1639 temperature_data_frame = append_vector(temperature_data_frame, vittelogenic_adults, "VITADULT"); 1289 temperature_data_frame = append_vector(temperature_data_frame, vittelogenic_adults, "VITADULT");
1640 # Standard error for vittelogenic adults. 1290 # Standard error for vittelogenic adults.
1668 F2.std_error = m_se[[6]]; 1318 F2.std_error = m_se[[6]];
1669 if (process_eggs) { 1319 if (process_eggs) {
1670 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);
1671 P_eggs = m_se[[1]]; 1321 P_eggs = m_se[[1]];
1672 P_eggs.std_error = m_se[[2]]; 1322 P_eggs.std_error = m_se[[2]];
1673 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");
1674 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");
1675 F1_eggs = m_se[[3]]; 1325 F1_eggs = m_se[[3]];
1676 F1_eggs.std_error = m_se[[4]]; 1326 F1_eggs.std_error = m_se[[4]];
1677 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");
1678 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");
1679 F2_eggs = m_se[[5]]; 1329 F2_eggs = m_se[[5]];
1680 F2_eggs.std_error = m_se[[6]]; 1330 F2_eggs.std_error = m_se[[6]];
1681 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");
1682 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");
1683 } 1333 }
1684 if (process_young_nymphs) { 1334 if (process_young_nymphs) {
1685 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);
1686 P_young_nymphs = m_se[[1]]; 1336 P_young_nymphs = m_se[[1]];
1687 P_young_nymphs.std_error = m_se[[2]]; 1337 P_young_nymphs.std_error = m_se[[2]];
1688 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");
1689 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");
1690 F1_young_nymphs = m_se[[3]]; 1340 F1_young_nymphs = m_se[[3]];
1691 F1_young_nymphs.std_error = m_se[[4]]; 1341 F1_young_nymphs.std_error = m_se[[4]];
1692 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");
1693 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");
1694 F2_young_nymphs = m_se[[5]]; 1344 F2_young_nymphs = m_se[[5]];
1695 F2_young_nymphs.std_error = m_se[[6]]; 1345 F2_young_nymphs.std_error = m_se[[6]];
1696 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");
1697 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");
1698 } 1348 }
1699 if (process_old_nymphs) { 1349 if (process_old_nymphs) {
1700 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);
1701 P_old_nymphs = m_se[[1]]; 1351 P_old_nymphs = m_se[[1]];
1702 P_old_nymphs.std_error = m_se[[2]]; 1352 P_old_nymphs.std_error = m_se[[2]];
1703 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");
1704 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");
1705 F1_old_nymphs = m_se[[3]]; 1355 F1_old_nymphs = m_se[[3]];
1706 F1_old_nymphs.std_error = m_se[[4]]; 1356 F1_old_nymphs.std_error = m_se[[4]];
1707 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");
1708 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");
1709 F2_old_nymphs = m_se[[5]]; 1359 F2_old_nymphs = m_se[[5]];
1710 F2_old_nymphs.std_error = m_se[[6]]; 1360 F2_old_nymphs.std_error = m_se[[6]];
1711 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");
1712 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");
1713 } 1363 }
1714 if (process_total_nymphs) { 1364 if (process_total_nymphs) {
1715 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);
1716 P_total_nymphs = m_se[[1]]; 1366 P_total_nymphs = m_se[[1]];
1717 P_total_nymphs.std_error = m_se[[2]]; 1367 P_total_nymphs.std_error = m_se[[2]];
1718 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");
1719 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");
1720 F1_total_nymphs = m_se[[3]]; 1370 F1_total_nymphs = m_se[[3]];
1721 F1_total_nymphs.std_error = m_se[[4]]; 1371 F1_total_nymphs.std_error = m_se[[4]];
1722 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");
1723 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");
1724 F2_total_nymphs = m_se[[5]]; 1374 F2_total_nymphs = m_se[[5]];
1725 F2_total_nymphs.std_error = m_se[[6]]; 1375 F2_total_nymphs.std_error = m_se[[6]];
1726 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");
1727 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");
1728 } 1378 }
1729 if (process_previttelogenic_adults) { 1379 if (process_previttelogenic_adults) {
1730 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);
1731 P_previttelogenic_adults = m_se[[1]]; 1381 P_previttelogenic_adults = m_se[[1]];
1732 P_previttelogenic_adults.std_error = m_se[[2]]; 1382 P_previttelogenic_adults.std_error = m_se[[2]];
1733 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");
1734 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");
1735 F1_previttelogenic_adults = m_se[[3]]; 1385 F1_previttelogenic_adults = m_se[[3]];
1736 F1_previttelogenic_adults.std_error = m_se[[4]]; 1386 F1_previttelogenic_adults.std_error = m_se[[4]];
1737 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");
1738 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");
1739 F2_previttelogenic_adults = m_se[[5]]; 1389 F2_previttelogenic_adults = m_se[[5]];
1740 F2_previttelogenic_adults.std_error = m_se[[6]]; 1390 F2_previttelogenic_adults.std_error = m_se[[6]];
1741 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");
1742 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");
1743 } 1393 }
1744 if (process_vittelogenic_adults) { 1394 if (process_vittelogenic_adults) {
1745 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);
1746 P_vittelogenic_adults = m_se[[1]]; 1396 P_vittelogenic_adults = m_se[[1]];
1747 P_vittelogenic_adults.std_error = m_se[[2]]; 1397 P_vittelogenic_adults.std_error = m_se[[2]];
1748 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");
1749 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");
1750 F1_vittelogenic_adults = m_se[[3]]; 1400 F1_vittelogenic_adults = m_se[[3]];
1751 F1_vittelogenic_adults.std_error = m_se[[4]]; 1401 F1_vittelogenic_adults.std_error = m_se[[4]];
1752 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");
1753 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");
1754 F2_vittelogenic_adults = m_se[[5]]; 1404 F2_vittelogenic_adults = m_se[[5]];
1755 F2_vittelogenic_adults.std_error = m_se[[6]]; 1405 F2_vittelogenic_adults.std_error = m_se[[6]];
1756 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");
1757 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");
1758 } 1408 }
1759 if (process_diapausing_adults) { 1409 if (process_diapausing_adults) {
1760 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);
1761 P_diapausing_adults = m_se[[1]]; 1411 P_diapausing_adults = m_se[[1]];
1762 P_diapausing_adults.std_error = m_se[[2]]; 1412 P_diapausing_adults.std_error = m_se[[2]];
1763 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");
1764 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");
1765 F1_diapausing_adults = m_se[[3]]; 1415 F1_diapausing_adults = m_se[[3]];
1766 F1_diapausing_adults.std_error = m_se[[4]]; 1416 F1_diapausing_adults.std_error = m_se[[4]];
1767 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");
1768 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");
1769 F2_diapausing_adults = m_se[[5]]; 1419 F2_diapausing_adults = m_se[[5]];
1770 F2_diapausing_adults.std_error = m_se[[6]]; 1420 F2_diapausing_adults.std_error = m_se[[6]];
1771 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");
1772 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");
1773 } 1423 }
1774 if (process_total_adults) { 1424 if (process_total_adults) {
1775 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);
1776 P_total_adults = m_se[[1]]; 1426 P_total_adults = m_se[[1]];
1777 P_total_adults.std_error = m_se[[2]]; 1427 P_total_adults.std_error = m_se[[2]];
1778 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");
1779 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");
1780 F1_total_adults = m_se[[3]]; 1430 F1_total_adults = m_se[[3]];
1781 F1_total_adults.std_error = m_se[[4]]; 1431 F1_total_adults.std_error = m_se[[4]];
1782 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");
1783 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");
1784 F2_total_adults = m_se[[5]]; 1434 F2_total_adults = m_se[[5]];
1785 F2_total_adults.std_error = m_se[[6]]; 1435 F2_total_adults.std_error = m_se[[6]];
1786 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");
1787 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");
1788 } 1438 }
1789 } 1439 }
1790 1440
1791 # Save the analyzed data for combined generations. 1441 # Save the analyzed data for combined generations.
1792 file_path = paste("output_data_dir", "04_combined_generations.csv", sep="/"); 1442 file_path = paste("output_data_dir", "04_combined_generations.csv", sep="/");
1821 dev.off(); 1471 dev.off();
1822 } else if (life_stage == "Nymph") { 1472 } else if (life_stage == "Nymph") {
1823 for (life_stage_nymph in life_stages_nymph) { 1473 for (life_stage_nymph in life_stages_nymph) {
1824 # Start PDF device driver. 1474 # Start PDF device driver.
1825 dev.new(width=20, height=30); 1475 dev.new(width=20, height=30);
1826 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)
1827 pdf(file=file_path, width=20, height=30, bg="white"); 1477 pdf(file=file_path, width=20, height=30, bg="white");
1828 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1478 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1829 if (life_stage_nymph=="Young") { 1479 if (life_stage_nymph=="Young") {
1830 # Young nymph population size by generation. 1480 # Young nymph population size by generation.
1831 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;
1854 group3 = F2_total_nymphs; 1504 group3 = F2_total_nymphs;
1855 group3_std_error = F2_total_nymphs.std_error; 1505 group3_std_error = F2_total_nymphs.std_error;
1856 } 1506 }
1857 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,
1858 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,
1859 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);
1860 # Turn off device driver to flush output. 1510 # Turn off device driver to flush output.
1861 dev.off(); 1511 dev.off();
1862 } 1512 }
1863 } else if (life_stage == "Adult") { 1513 } else if (life_stage == "Adult") {
1864 for (life_stage_adult in life_stages_adult) { 1514 for (life_stage_adult in life_stages_adult) {
1865 # Start PDF device driver. 1515 # Start PDF device driver.
1866 dev.new(width=20, height=30); 1516 dev.new(width=20, height=30);
1867 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)
1868 pdf(file=file_path, width=20, height=30, bg="white"); 1518 pdf(file=file_path, width=20, height=30, bg="white");
1869 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1519 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1870 if (life_stage_adult=="Pre-vittelogenic") { 1520 if (life_stage_adult=="Pre-vittelogenic") {
1871 # Pre-vittelogenic adult population size by generation. 1521 # Pre-vittelogenic adult population size by generation.
1872 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;
1904 group3 = F2_total_adults; 1554 group3 = F2_total_adults;
1905 group3_std_error = F2_total_adults.std_error; 1555 group3_std_error = F2_total_adults.std_error;
1906 } 1556 }
1907 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,
1908 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,
1909 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);
1910 # Turn off device driver to flush output. 1560 # Turn off device driver to flush output.
1911 dev.off(); 1561 dev.off();
1912 } 1562 }
1913 } else if (life_stage == "Total") { 1563 } else if (life_stage == "Total") {
1914 # Start PDF device driver. 1564 # Start PDF device driver.
1943 dev.off(); 1593 dev.off();
1944 } else if (life_stage == "Nymph") { 1594 } else if (life_stage == "Nymph") {
1945 for (life_stage_nymph in life_stages_nymph) { 1595 for (life_stage_nymph in life_stages_nymph) {
1946 # Start PDF device driver. 1596 # Start PDF device driver.
1947 dev.new(width=20, height=30); 1597 dev.new(width=20, height=30);
1948 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)
1949 pdf(file=file_path, width=20, height=30, bg="white"); 1599 pdf(file=file_path, width=20, height=30, bg="white");
1950 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1600 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1951 if (life_stage_nymph=="Total") { 1601 if (life_stage_nymph=="Total") {
1952 # Total nymph population size. 1602 # Total nymph population size.
1953 group = total_nymphs; 1603 group = total_nymphs;
1962 group_std_error = old_nymphs.std_error; 1612 group_std_error = old_nymphs.std_error;
1963 } 1613 }
1964 maxval = max(group+group_std_error) + 100; 1614 maxval = max(group+group_std_error) + 100;
1965 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,
1966 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,
1967 life_stages_nymph=life_stage_nymph); 1617 sub_life_stage=life_stage_nymph);
1968 # Turn off device driver to flush output. 1618 # Turn off device driver to flush output.
1969 dev.off(); 1619 dev.off();
1970 } 1620 }
1971 } else if (life_stage == "Adult") { 1621 } else if (life_stage == "Adult") {
1972 for (life_stage_adult in life_stages_adult) { 1622 for (life_stage_adult in life_stages_adult) {
1973 # Start PDF device driver. 1623 # Start PDF device driver.
1974 dev.new(width=20, height=30); 1624 dev.new(width=20, height=30);
1975 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)
1976 pdf(file=file_path, width=20, height=30, bg="white"); 1626 pdf(file=file_path, width=20, height=30, bg="white");
1977 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1627 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1978 if (life_stage_adult=="Total") { 1628 if (life_stage_adult=="Total") {
1979 # Total adult population size. 1629 # Total adult population size.
1980 group = total_adults; 1630 group = total_adults;
1993 group_std_error = diapausing_adults.std_error 1643 group_std_error = diapausing_adults.std_error
1994 } 1644 }
1995 maxval = max(group+group_std_error) + 100; 1645 maxval = max(group+group_std_error) + 100;
1996 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,
1997 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,
1998 life_stages_adult=life_stage_adult); 1648 sub_life_stage=life_stage_adult);
1999 # Turn off device driver to flush output. 1649 # Turn off device driver to flush output.
2000 dev.off(); 1650 dev.off();
2001 } 1651 }
2002 } else if (life_stage == "Total") { 1652 } else if (life_stage == "Total") {
2003 # Start PDF device driver. 1653 # Start PDF device driver.