Mercurial > repos > greg > extract_ipm_date_interval
comparison utils.R @ 0:1868b7913590 draft
Uploaded
author | greg |
---|---|
date | Tue, 07 Aug 2018 13:05:10 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:1868b7913590 |
---|---|
1 #!/usr/bin/env Rscript | |
2 | |
3 get_file_path = function(life_stage, base_name, sub_life_stage=NULL) { | |
4 if (is.null(sub_life_stage)) { | |
5 lsi = get_life_stage_index(life_stage); | |
6 file_name = paste(lsi, base_name, sep="_"); | |
7 } else { | |
8 lsi = get_life_stage_index(life_stage, sub_life_stage=sub_life_stage); | |
9 file_name = paste(lsi, tolower(sub_life_stage), base_name, sep="_"); | |
10 } | |
11 file_path = paste("output_plots_dir", file_name, sep="/"); | |
12 return(file_path); | |
13 } | |
14 | |
15 get_year_from_date = function(date_str) { | |
16 date_str_items = strsplit(date_str, "-")[[1]]; | |
17 return (date_str_items[1]); | |
18 } | |
19 | |
20 get_life_stage_index = function(life_stage, sub_life_stage=NULL) { | |
21 # Name collection elements so that they | |
22 # are displayed in logical order. | |
23 if (life_stage=="Egg") { | |
24 lsi = "01"; | |
25 } else if (life_stage=="Nymph") { | |
26 if (sub_life_stage=="Young") { | |
27 lsi = "02"; | |
28 } else if (sub_life_stage=="Old") { | |
29 lsi = "03"; | |
30 } else if (sub_life_stage=="Total") { | |
31 lsi="04"; | |
32 } | |
33 } else if (life_stage=="Adult") { | |
34 if (sub_life_stage=="Pre-vittelogenic") { | |
35 lsi = "05"; | |
36 } else if (sub_life_stage=="Vittelogenic") { | |
37 lsi = "06"; | |
38 } else if (sub_life_stage=="Diapausing") { | |
39 lsi = "07"; | |
40 } else if (sub_life_stage=="Total") { | |
41 lsi = "08"; | |
42 } | |
43 } else if (life_stage=="Total") { | |
44 lsi = "09"; | |
45 } | |
46 return(lsi); | |
47 } | |
48 | |
49 get_mean_and_std_error = function(p_replications, f1_replications, f2_replications) { | |
50 # P mean. | |
51 p_m = apply(p_replications, 1, mean); | |
52 # P standard error. | |
53 p_se = apply(p_replications, 1, sd) / sqrt(opt$replications); | |
54 # F1 mean. | |
55 f1_m = apply(f1_replications, 1, mean); | |
56 # F1 standard error. | |
57 f1_se = apply(f1_replications, 1, sd) / sqrt(opt$replications); | |
58 # F2 mean. | |
59 f2_m = apply(f2_replications, 1, mean); | |
60 # F2 standard error. | |
61 f2_se = apply(f2_replications, 1, sd) / sqrt(opt$replications); | |
62 return(list(p_m, p_se, f1_m, f1_se, f2_m, f2_se)) | |
63 } | |
64 | |
65 get_tick_index = function(index, last_tick, ticks, tick_labels, tick_sep) { | |
66 # The R code tries hard not to draw overlapping tick labels, and so | |
67 # will omit labels where they would abut or overlap previously drawn | |
68 # labels. This can result in, for example, every other tick being | |
69 # labelled. We'll keep track of the last tick to make sure all of | |
70 # the month labels are displayed, and missing ticks are restricted | |
71 # to Sundays which have no labels anyway. | |
72 if (last_tick==0) { | |
73 return(length(ticks)+1); | |
74 } | |
75 last_saved_tick = ticks[[length(ticks)]]; | |
76 if (index-last_saved_tick<tick_sep) { | |
77 last_saved_month = tick_labels[[length(tick_labels)]]; | |
78 if (last_saved_month=="") { | |
79 # We're safe overwriting a tick | |
80 # with no label (i.e., a Sunday tick). | |
81 return(length(ticks)); | |
82 } else { | |
83 # Don't eliminate a Month label. | |
84 return(NULL); | |
85 } | |
86 } | |
87 return(length(ticks)+1); | |
88 } | |
89 | |
90 get_total_days = function(is_leap_year) { | |
91 # Get the total number of days in the current year. | |
92 if (is_leap_year) { | |
93 return(366); | |
94 } else { | |
95 return(365); | |
96 } | |
97 } | |
98 | |
99 get_x_axis_ticks_and_labels = function(temperature_data_frame, prepend_end_doy_norm=0, append_start_doy_norm=0, date_interval=FALSE) { | |
100 # Generate a list of ticks and labels for plotting the x axis. | |
101 if (prepend_end_doy_norm > 0) { | |
102 prepend_end_norm_row = which(temperature_data_frame$DOY==prepend_end_doy_norm); | |
103 } else { | |
104 prepend_end_norm_row = 0; | |
105 } | |
106 if (append_start_doy_norm > 0) { | |
107 append_start_norm_row = which(temperature_data_frame$DOY==append_start_doy_norm); | |
108 } else { | |
109 append_start_norm_row = 0; | |
110 } | |
111 num_rows = dim(temperature_data_frame)[1]; | |
112 tick_labels = list(); | |
113 ticks = list(); | |
114 current_month_label = NULL; | |
115 last_tick = 0; | |
116 if (date_interval) { | |
117 tick_sep = 0; | |
118 } else { | |
119 tick_sep = 3; | |
120 } | |
121 for (i in 1:num_rows) { | |
122 # Get the year and month from the date which | |
123 # has the format YYYY-MM-DD. | |
124 date = format(temperature_data_frame$DATE[i]); | |
125 # Get the month label. | |
126 items = strsplit(date, "-")[[1]]; | |
127 month = items[2]; | |
128 month_label = month.abb[as.integer(month)]; | |
129 day = as.integer(items[3]); | |
130 doy = as.integer(temperature_data_frame$DOY[i]); | |
131 # We're plotting the entire year, so ticks will | |
132 # occur on Sundays and the first of each month. | |
133 if (i == prepend_end_norm_row) { | |
134 # Add a tick for the end of the 30 year normnals data | |
135 # that was prepended to the year-to-date data. | |
136 label_str = "End prepended 30 year normals"; | |
137 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep) | |
138 ticks[tick_index] = i; | |
139 if (date_interval) { | |
140 # Append the day to label_str | |
141 tick_labels[tick_index] = paste(label_str, day, sep=" "); | |
142 } else { | |
143 tick_labels[tick_index] = label_str; | |
144 } | |
145 last_tick = i; | |
146 } else if (doy == append_start_doy_norm) { | |
147 # Add a tick for the start of the 30 year normnals data | |
148 # that was appended to the year-to-date data. | |
149 label_str = "Start appended 30 year normals"; | |
150 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep) | |
151 ticks[tick_index] = i; | |
152 if (!identical(current_month_label, month_label)) { | |
153 # Append the month to label_str. | |
154 label_str = paste(label_str, month_label, spe=" "); | |
155 current_month_label = month_label; | |
156 } | |
157 if (date_interval) { | |
158 # Append the day to label_str | |
159 label_str = paste(label_str, day, sep=" "); | |
160 } | |
161 tick_labels[tick_index] = label_str; | |
162 last_tick = i; | |
163 } else if (i==num_rows) { | |
164 # Add a tick for the last day of the year. | |
165 label_str = ""; | |
166 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep) | |
167 ticks[tick_index] = i; | |
168 if (!identical(current_month_label, month_label)) { | |
169 # Append the month to label_str. | |
170 label_str = month_label; | |
171 current_month_label = month_label; | |
172 } | |
173 if (date_interval) { | |
174 # Append the day to label_str | |
175 label_str = paste(label_str, day, sep=" "); | |
176 } | |
177 tick_labels[tick_index] = label_str; | |
178 } else { | |
179 if (!identical(current_month_label, month_label)) { | |
180 # Add a tick for the month. | |
181 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep) | |
182 ticks[tick_index] = i; | |
183 if (date_interval) { | |
184 # Append the day to the month. | |
185 tick_labels[tick_index] = paste(month_label, day, sep=" "); | |
186 } else { | |
187 tick_labels[tick_index] = month_label; | |
188 } | |
189 current_month_label = month_label; | |
190 last_tick = i; | |
191 } | |
192 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep) | |
193 if (!is.null(tick_index)) { | |
194 if (date_interval) { | |
195 # Add a tick for every day. The first tick is the | |
196 # month label, so add a tick only if i is not 1 | |
197 if (i>1 & day>1) { | |
198 tick_index = get_tick_index(i, last_tick, ticks, tick_labels, tick_sep) | |
199 ticks[tick_index] = i; | |
200 # Add the day as the label. | |
201 tick_labels[tick_index] = day; | |
202 last_tick = i; | |
203 } | |
204 } else { | |
205 # Get the day. | |
206 day = weekdays(as.Date(date)); | |
207 if (day=="Sunday") { | |
208 # Add a tick if we're on a Sunday. | |
209 ticks[tick_index] = i; | |
210 # Add a blank month label so it is not displayed. | |
211 tick_labels[tick_index] = ""; | |
212 last_tick = i; | |
213 } | |
214 } | |
215 } | |
216 } | |
217 } | |
218 return(list(ticks, tick_labels)); | |
219 } | |
220 | |
221 render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval, | |
222 replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL, | |
223 sub_life_stage=NULL) { | |
224 if (chart_type=="pop_size_by_life_stage") { | |
225 if (life_stage=="Total") { | |
226 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); | |
227 legend_text = c("Egg", "Nymph", "Adult"); | |
228 columns = c(4, 2, 1); | |
229 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); | |
230 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3); | |
231 lines(days, group2, lwd=2, lty=1, col=2); | |
232 lines(days, group3, lwd=2, lty=1, col=4); | |
233 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); | |
234 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); | |
235 if (plot_std_error=="yes") { | |
236 # Standard error for group. | |
237 lines(days, group+group_std_error, lty=2); | |
238 lines(days, group-group_std_error, lty=2); | |
239 # Standard error for group2. | |
240 lines(days, group2+group2_std_error, col=2, lty=2); | |
241 lines(days, group2-group2_std_error, col=2, lty=2); | |
242 # Standard error for group3. | |
243 lines(days, group3+group3_std_error, col=4, lty=2); | |
244 lines(days, group3-group3_std_error, col=4, lty=2); | |
245 } | |
246 } else { | |
247 if (life_stage=="Egg") { | |
248 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); | |
249 legend_text = c(life_stage); | |
250 columns = c(4); | |
251 } else if (life_stage=="Nymph") { | |
252 stage = paste(sub_life_stage, "Nymph Pop :", sep=" "); | |
253 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); | |
254 legend_text = c(paste(sub_life_stage, life_stage, sep=" ")); | |
255 columns = c(2); | |
256 } else if (life_stage=="Adult") { | |
257 stage = paste(sub_life_stage, "Adult Pop", sep=" "); | |
258 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); | |
259 legend_text = c(paste(sub_life_stage, life_stage, sep=" ")); | |
260 columns = c(1); | |
261 } | |
262 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); | |
263 legend("topleft", legend_text, lty=c(1), col="black", cex=3); | |
264 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); | |
265 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); | |
266 if (plot_std_error=="yes") { | |
267 # Standard error for group. | |
268 lines(days, group+group_std_error, lty=2); | |
269 lines(days, group-group_std_error, lty=2); | |
270 } | |
271 } | |
272 } else if (chart_type=="pop_size_by_generation") { | |
273 if (life_stage=="Total") { | |
274 title_str = ": Total Pop by Gen :"; | |
275 } else if (life_stage=="Egg") { | |
276 title_str = ": Egg Pop by Gen :"; | |
277 } else if (life_stage=="Nymph") { | |
278 title_str = paste(":", sub_life_stage, "Nymph Pop by Gen", ":", sep=" "); | |
279 } else if (life_stage=="Adult") { | |
280 title_str = paste(":", sub_life_stage, "Adult Pop by Gen", ":", sep=" "); | |
281 } | |
282 title = paste(insect, ": Reps", replications, title_str, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); | |
283 legend_text = c("P", "F1", "F2"); | |
284 columns = c(1, 2, 4); | |
285 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); | |
286 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3); | |
287 lines(days, group2, lwd=2, lty=1, col=2); | |
288 lines(days, group3, lwd=2, lty=1, col=4); | |
289 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); | |
290 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); | |
291 if (plot_std_error=="yes") { | |
292 # Standard error for group. | |
293 lines(days, group+group_std_error, lty=2); | |
294 lines(days, group-group_std_error, lty=2); | |
295 # Standard error for group2. | |
296 lines(days, group2+group2_std_error, col=2, lty=2); | |
297 lines(days, group2-group2_std_error, col=2, lty=2); | |
298 # Standard error for group3. | |
299 lines(days, group3+group3_std_error, col=4, lty=2); | |
300 lines(days, group3-group3_std_error, col=4, lty=2); | |
301 } | |
302 } | |
303 } | |
304 | |
305 stop_err = function(msg) { | |
306 cat(msg, file=stderr()); | |
307 quit(save="no", status=1); | |
308 } | |
309 | |
310 validate_date = function(date_str) { | |
311 valid_date = as.Date(date_str, format="%Y-%m-%d"); | |
312 if( class(valid_date)=="try-error" || is.na(valid_date)) { | |
313 msg = paste("Invalid date: ", date_str, ", valid date format is yyyy-mm-dd.", sep=""); | |
314 stop_err(msg); | |
315 } | |
316 return(valid_date); | |
317 } |