Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 19:3c6e94e477cb draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Tue, 06 Mar 2018 14:26:45 -0500
parents f5ecff4800f2
children 214217142600
comparison
equal deleted inserted replaced
18:f5ecff4800f2 19:3c6e94e477cb
50 # Append daylight_length_vector as a new column to temperature_data_frame. 50 # Append daylight_length_vector as a new column to temperature_data_frame.
51 temperature_data_frame[, num_columns+1] = daylight_length_vector; 51 temperature_data_frame[, num_columns+1] = daylight_length_vector;
52 return(temperature_data_frame); 52 return(temperature_data_frame);
53 } 53 }
54 54
55 dev.egg = function(temperature) {
56 # This function is not used, but was
57 # in the original, so keep it for now.
58 dev.rate = -0.9843 * temperature + 33.438;
59 return(dev.rate);
60 }
61
62 dev.emerg = function(temperature) {
63 # This function is not used, but was
64 # in the original, so keep it for now.
65 emerg.rate = -0.5332 * temperature + 24.147;
66 return(emerg.rate);
67 }
68
69 dev.old = function(temperature) {
70 # This function is not used, but was
71 # in the original, so keep it for now.
72 n34 = -0.6119 * temperature + 17.602;
73 n45 = -0.4408 * temperature + 19.036;
74 dev.rate = mean(n34 + n45);
75 return(dev.rate);
76 }
77
78 dev.young = function(temperature) {
79 # This function is not used, but was
80 # in the original, so keep it for now.
81 n12 = -0.3728 * temperature + 14.68;
82 n23 = -0.6119 * temperature + 25.249;
83 dev.rate = mean(n12 + n23);
84 return(dev.rate);
85 }
86
87 get_date_labels = function(temperature_data_frame, num_rows) { 55 get_date_labels = function(temperature_data_frame, num_rows) {
88 # Keep track of the years to see if spanning years. 56 # Keep track of the years to see if spanning years.
89 month_labels = list(); 57 month_labels = list();
90 current_month_label = NULL; 58 current_month_label = NULL;
91 for (i in 1:num_rows) { 59 for (i in 1:num_rows) {
99 month_labels[length(month_labels)+1] = month_label; 67 month_labels[length(month_labels)+1] = month_label;
100 current_month_label = month_label; 68 current_month_label = month_label;
101 } 69 }
102 } 70 }
103 return(c(unlist(month_labels))); 71 return(c(unlist(month_labels)));
72 }
73
74
75 get_file_path = function(life_stage, base_name, life_stage_nymph=NULL, life_stage_adult=NULL) {
76 if (!is.null(life_stage_nymph)) {
77 lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph);
78 file_name = paste(lsi, tolower(life_stage_nymph), base_name, sep="_");
79 } else if (!is.null(life_stage_adult)) {
80 lsi = get_life_stage_index(life_stage, life_stage_adult=life_stage_adult);
81 file_name = paste(lsi, tolower(life_stage_adult), base_name, sep="_");
82 } else {
83 lsi = get_life_stage_index(life_stage);
84 file_name = paste(lsi, base_name, sep="_");
85 }
86 file_path = paste("output_dir", file_name, sep="/");
87 return(file_path);
104 } 88 }
105 89
106 get_life_stage_index = function(life_stage, life_stage_nymph=NULL, life_stage_adult=NULL) { 90 get_life_stage_index = function(life_stage, life_stage_nymph=NULL, life_stage_adult=NULL) {
107 # Name collection elements so that they 91 # Name collection elements so that they
108 # are displayed in logical order. 92 # are displayed in logical order.
344 } else { 328 } else {
345 plot_generations_separately = FALSE; 329 plot_generations_separately = FALSE;
346 } 330 }
347 # Read the temperature data into a data frame. 331 # Read the temperature data into a data frame.
348 temperature_data_frame = parse_input_data(opt$input, opt$num_days); 332 temperature_data_frame = parse_input_data(opt$input, opt$num_days);
349 output_dir = "output_dir";
350 # Get the date labels for plots. 333 # Get the date labels for plots.
351 date_labels = get_date_labels(temperature_data_frame, opt$num_days); 334 date_labels = get_date_labels(temperature_data_frame, opt$num_days);
352 # All latitude values are the same, so get the value for plots from the first row. 335 # All latitude values are the same, so get the value for plots from the first row.
353 latitude = temperature_data_frame$LATITUDE[1]; 336 latitude = temperature_data_frame$LATITUDE[1];
354 # Get the number of days for plots. 337 # Get the number of days for plots.
966 if (plot_generations_separately) { 949 if (plot_generations_separately) {
967 for (life_stage in life_stages) { 950 for (life_stage in life_stages) {
968 if (life_stage == "Egg") { 951 if (life_stage == "Egg") {
969 # Start PDF device driver. 952 # Start PDF device driver.
970 dev.new(width=20, height=30); 953 dev.new(width=20, height=30);
971 lsi = get_life_stage_index(life_stage); 954 file_path = get_file_path(life_stage, "egg_pop_by_generation.pdf")
972 file_name = paste(lsi, "egg_pop_by_generation.pdf", sep="_");
973 file_path = paste(output_dir, file_name, sep="/");
974 pdf(file=file_path, width=20, height=30, bg="white"); 955 pdf(file=file_path, width=20, height=30, bg="white");
975 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 956 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
976 # Egg population size by generation. 957 # Egg population size by generation.
977 maxval = max(P_eggs+F1_eggs+F2_eggs) + 100; 958 maxval = max(P_eggs+F1_eggs+F2_eggs) + 100;
978 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 959 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
982 dev.off(); 963 dev.off();
983 } else if (life_stage == "Nymph") { 964 } else if (life_stage == "Nymph") {
984 for (life_stage_nymph in life_stages_nymph) { 965 for (life_stage_nymph in life_stages_nymph) {
985 # Start PDF device driver. 966 # Start PDF device driver.
986 dev.new(width=20, height=30); 967 dev.new(width=20, height=30);
987 lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph); 968 file_path = get_file_path(life_stage, "nymph_pop_by_generation.pdf", life_stage_nymph=life_stage_nymph)
988 file_name = paste(lsi, tolower(life_stage_nymph), "nymph_pop_by_generation.pdf", sep="_");
989 file_path = paste(output_dir, file_name, sep="/");
990 pdf(file=file_path, width=20, height=30, bg="white"); 969 pdf(file=file_path, width=20, height=30, bg="white");
991 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 970 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
992 # Nymph population size by generation. 971 # Nymph population size by generation.
993 maxval = max(P_nymphs+F1_nymphs+F2_nymphs) + 100; 972 maxval = max(P_nymphs+F1_nymphs+F2_nymphs) + 100;
994 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 973 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
995 opt$replications, life_stage, group=P_nymphs, group_std_error=P_nymphs.std_error, group2=F1_nymphs, group2_std_error=F1_nymphs.std_error, 974 opt$replications, life_stage, group=P_nymphs, group_std_error=P_nymphs.std_error, group2=F1_nymphs, group2_std_error=F1_nymphs.std_error,
996 group3=F2_nymphs, group3_std_error=F2_nymphs.std_error, life_stages_nymph=life_stage_nymph); 975 group3=F2_nymphs, group3_std_error=F2_nymphs.std_error, life_stages_nymph=life_stage_nymph);
997 # Turn off device driver to flush output. 976 # Turn off device driver to flush output.
998 dev.off(); 977 dev.off();
999 } 978 }
1000 } else if (life_stage == "Adult") { 979 } else if (life_stage == "Adult") {
1001 for (life_stage_adult in life_stages_adult) { 980 for (life_stage_adult in life_stages_adult) {
1002 # Start PDF device driver. 981 # Start PDF device driver.
1003 dev.new(width=20, height=30); 982 dev.new(width=20, height=30);
1004 lsi = get_life_stage_index(life_stage, life_stage_adult=life_stage_adult); 983 file_path = get_file_path(life_stage, "adult_pop_by_generation.pdf", life_stage_adult=life_stage_adult)
1005 file_name = paste(lsi, tolower(life_stage_adult), "adult_pop_by_generation.pdf", sep="_");
1006 file_path = paste(output_dir, file_name, sep="/");
1007 pdf(file=file_path, width=20, height=30, bg="white"); 984 pdf(file=file_path, width=20, height=30, bg="white");
1008 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 985 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1009 # Adult population size by generation. 986 # Adult population size by generation.
1010 maxval = max(P_adults+F1_adults+F2_adults) + 100; 987 maxval = max(P_adults+F1_adults+F2_adults) + 100;
1011 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 988 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
1017 } else if (life_stage == "Total") { 994 } else if (life_stage == "Total") {
1018 # Start PDF device driver. 995 # Start PDF device driver.
1019 # Name collection elements so that they 996 # Name collection elements so that they
1020 # are displayed in logical order. 997 # are displayed in logical order.
1021 dev.new(width=20, height=30); 998 dev.new(width=20, height=30);
1022 lsi = get_life_stage_index(life_stage); 999 file_path = get_file_path(life_stage, "total_pop_by_generation.pdf")
1023 file_name = paste(lsi, "total_pop_by_generation.pdf", sep="_");
1024 file_path = paste(output_dir, file_name, sep="/");
1025 pdf(file=file_path, width=20, height=30, bg="white"); 1000 pdf(file=file_path, width=20, height=30, bg="white");
1026 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1001 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1027 # Total population size by generation. 1002 # Total population size by generation.
1028 maxval = max(P+F1+F2) + 100; 1003 maxval = max(P+F1+F2) + 100;
1029 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1004 render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
1035 } else { 1010 } else {
1036 for (life_stage in life_stages) { 1011 for (life_stage in life_stages) {
1037 if (life_stage == "Egg") { 1012 if (life_stage == "Egg") {
1038 # Start PDF device driver. 1013 # Start PDF device driver.
1039 dev.new(width=20, height=30); 1014 dev.new(width=20, height=30);
1040 lsi = get_life_stage_index(life_stage); 1015 file_path = get_file_path(life_stage, "egg_pop.pdf")
1041 file_name = paste(lsi, "egg_pop.pdf", sep="_");
1042 file_path = paste(output_dir, file_name, sep="/");
1043 pdf(file=file_path, width=20, height=30, bg="white"); 1016 pdf(file=file_path, width=20, height=30, bg="white");
1044 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1017 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1045 # Egg population size. 1018 # Egg population size.
1046 maxval = max(eggs+eggs.std_error) + 100; 1019 maxval = max(eggs+eggs.std_error) + 100;
1047 render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1020 render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,
1050 dev.off(); 1023 dev.off();
1051 } else if (life_stage == "Nymph") { 1024 } else if (life_stage == "Nymph") {
1052 for (life_stage_nymph in life_stages_nymph) { 1025 for (life_stage_nymph in life_stages_nymph) {
1053 # Start PDF device driver. 1026 # Start PDF device driver.
1054 dev.new(width=20, height=30); 1027 dev.new(width=20, height=30);
1055 lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph); 1028 file_path = get_file_path(life_stage, "nymph_pop.pdf", life_stage_nymph=life_stage_nymph)
1056 file_name = paste(lsi, tolower(life_stage_nymph), "nymph_pop.pdf", sep="_");
1057 file_path = paste(output_dir, file_name, sep="/");
1058 pdf(file=file_path, width=20, height=30, bg="white"); 1029 pdf(file=file_path, width=20, height=30, bg="white");
1059 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1030 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1060 if (life_stage_nymph=="Total") { 1031 if (life_stage_nymph=="Total") {
1061 # Total nymph population size. 1032 # Total nymph population size.
1062 group = total_nymphs; 1033 group = total_nymphs;
1078 } 1049 }
1079 } else if (life_stage == "Adult") { 1050 } else if (life_stage == "Adult") {
1080 for (life_stage_adult in life_stages_adult) { 1051 for (life_stage_adult in life_stages_adult) {
1081 # Start PDF device driver. 1052 # Start PDF device driver.
1082 dev.new(width=20, height=30); 1053 dev.new(width=20, height=30);
1083 lsi = get_life_stage_index(life_stage, life_stage_adult=life_stage_adult); 1054 file_path = get_file_path(life_stage, "adult_pop.pdf", life_stage_adult=life_stage_adult)
1084 file_name = paste(lsi, tolower(life_stage_adult), "adult_pop.pdf", sep="_");
1085 file_path = paste(output_dir, file_name, sep="/");
1086 pdf(file=file_path, width=20, height=30, bg="white"); 1055 pdf(file=file_path, width=20, height=30, bg="white");
1087 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1056 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1088 if (life_stage_adult=="Total") { 1057 if (life_stage_adult=="Total") {
1089 # Total adult population size. 1058 # Total adult population size.
1090 group = total_adults; 1059 group = total_adults;
1109 dev.off(); 1078 dev.off();
1110 } 1079 }
1111 } else if (life_stage == "Total") { 1080 } else if (life_stage == "Total") {
1112 # Start PDF device driver. 1081 # Start PDF device driver.
1113 dev.new(width=20, height=30); 1082 dev.new(width=20, height=30);
1114 lsi = get_life_stage_index(life_stage); 1083 file_path = get_file_path(life_stage, "total_pop.pdf")
1115 file_name = paste(lsi, "total_pop.pdf", sep="_");
1116 file_path = paste(output_dir, file_name, sep="/");
1117 pdf(file=file_path, width=20, height=30, bg="white"); 1084 pdf(file=file_path, width=20, height=30, bg="white");
1118 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1085 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1119 # Total population size. 1086 # Total population size.
1120 maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100; 1087 maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100;
1121 render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, 1088 render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval,