Mercurial > repos > greg > insect_phenology_model
comparison insect_phenology_model.R @ 35:29ec818b1c29 draft
Uploaded
author | greg |
---|---|
date | Tue, 20 Mar 2018 09:45:26 -0400 |
parents | 7aa848b0e55c |
children | 5097cfeedc4f |
comparison
equal
deleted
inserted
replaced
34:7aa848b0e55c | 35:29ec818b1c29 |
---|---|
58 # Append vector vec as a new column to data_frame. | 58 # Append vector vec as a new column to data_frame. |
59 data_frame[,num_columns+1] = vec; | 59 data_frame[,num_columns+1] = vec; |
60 # Reset the column names with the additional column for later access. | 60 # Reset the column names with the additional column for later access. |
61 colnames(data_frame) = append(current_column_names, new_column_name); | 61 colnames(data_frame) = append(current_column_names, new_column_name); |
62 return(data_frame); | 62 return(data_frame); |
63 } | |
64 | |
65 get_x_axis_ticks_and_labels = function(temperature_data_frame, num_rows) { | |
66 # Keep track of the years to see if spanning years. | |
67 month_labels = list(); | |
68 ticks = list(); | |
69 current_month_label = NULL; | |
70 for (i in 1:num_rows) { | |
71 # Get the year and month from the date which | |
72 # has the format YYYY-MM-DD. | |
73 date = format(temperature_data_frame$DATE[i]); | |
74 # Get the month label. | |
75 items = strsplit(date, "-")[[1]]; | |
76 month = items[2]; | |
77 month_label = month.abb[as.integer(month)]; | |
78 if (!identical(current_month_label, month_label)) { | |
79 # Add an x-axis tick for the month. | |
80 ticks[length(ticks)+1] = i; | |
81 month_labels[length(month_labels)+1] = month_label; | |
82 current_month_label = month_label; | |
83 } | |
84 # Get the day. | |
85 day = weekdays(as.Date(date)); | |
86 if (day=="Sunday") { | |
87 # Add an x-axis tick if we're on a Sunday. | |
88 ticks[length(ticks)+1] = i; | |
89 # Add a blank month label so it is not displayed. | |
90 month_labels[length(month_labels)+1] = ""; | |
91 } | |
92 } | |
93 return(list(ticks, month_labels)); | |
94 } | 63 } |
95 | 64 |
96 get_file_path = function(life_stage, base_name, life_stage_nymph=NULL, life_stage_adult=NULL) { | 65 get_file_path = function(life_stage, base_name, life_stage_nymph=NULL, life_stage_adult=NULL) { |
97 if (!is.null(life_stage_nymph)) { | 66 if (!is.null(life_stage_nymph)) { |
98 lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph); | 67 lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph); |
222 averages = sum(dh) / 24; | 191 averages = sum(dh) / 24; |
223 } | 192 } |
224 return(c(curr_mean_temp, averages)) | 193 return(c(curr_mean_temp, averages)) |
225 } | 194 } |
226 | 195 |
196 get_tick_index = function(index, last_tick, ticks, month_labels) { | |
197 # The R code tries hard not to draw overlapping tick labels, and so | |
198 # will omit labels where they would abut or overlap previously drawn | |
199 # labels. This can result in, for example, every other tick being | |
200 # labelled. We'll keep track of the last tick to make sure all of | |
201 # the month labels are displayed, and missing ticks are restricted | |
202 # to Sundays which have no labels anyway. | |
203 if (last_tick==0) { | |
204 return(length(ticks)+1); | |
205 } | |
206 last_saved_tick = ticks[[length(ticks)]]; | |
207 if (index-last_saved_tick<6) { | |
208 last_saved_month = month_labels[[length(month_labels)]]; | |
209 if (last_saved_month=="") { | |
210 # We're safe overwriting a tick | |
211 # with no label (i.e., a Sunday tick). | |
212 return(length(ticks)); | |
213 } else { | |
214 # Don't eliminate a Month label. | |
215 return(NULL); | |
216 } | |
217 } | |
218 return(length(ticks)+1); | |
219 } | |
220 | |
221 get_x_axis_ticks_and_labels = function(temperature_data_frame, num_rows) { | |
222 # Keep track of the years to see if spanning years. | |
223 month_labels = list(); | |
224 ticks = list(); | |
225 current_month_label = NULL; | |
226 last_tick = 0; | |
227 for (i in 1:num_rows) { | |
228 # Get the year and month from the date which | |
229 # has the format YYYY-MM-DD. | |
230 date = format(temperature_data_frame$DATE[i]); | |
231 # Get the month label. | |
232 items = strsplit(date, "-")[[1]]; | |
233 month = items[2]; | |
234 month_label = month.abb[as.integer(month)]; | |
235 tick_index = get_tick_index(i, last_tick, ticks, month_labels) | |
236 if (!is.null(tick_index)) { | |
237 if (!identical(current_month_label, month_label)) { | |
238 # Add an x-axis tick for the month. | |
239 ticks[tick_index] = i; | |
240 month_labels[tick_index] = month_label; | |
241 current_month_label = month_label; | |
242 last_tick = i; | |
243 } | |
244 # Get the day. | |
245 day = weekdays(as.Date(date)); | |
246 if (day=="Sunday") { | |
247 # Add an x-axis tick if we're on a Sunday. | |
248 ticks[tick_index] = i; | |
249 # Add a blank month label so it is not displayed. | |
250 month_labels[tick_index] = ""; | |
251 last_tick = i; | |
252 } | |
253 } | |
254 } | |
255 return(list(ticks, month_labels)); | |
256 } | |
257 | |
227 mortality.adult = function(temperature) { | 258 mortality.adult = function(temperature) { |
228 if (temperature < 12.7) { | 259 if (temperature < 12.7) { |
229 mortality.probability = 0.002; | 260 mortality.probability = 0.002; |
230 } | 261 } |
231 else { | 262 else { |
279 if (chart_type=="pop_size_by_life_stage") { | 310 if (chart_type=="pop_size_by_life_stage") { |
280 if (life_stage=="Total") { | 311 if (life_stage=="Total") { |
281 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); | 312 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); |
282 legend_text = c("Egg", "Nymph", "Adult"); | 313 legend_text = c("Egg", "Nymph", "Adult"); |
283 columns = c(4, 2, 1); | 314 columns = c(4, 2, 1); |
284 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3); | 315 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); |
285 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3); | 316 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3); |
286 lines(days, group2, lwd=2, lty=1, col=2); | 317 lines(days, group2, lwd=2, lty=1, col=2); |
287 lines(days, group3, lwd=2, lty=1, col=4); | 318 lines(days, group3, lwd=2, lty=1, col=4); |
288 axis(side=1, at=ticks, labels=date_labels); | 319 axis(side=1, at=ticks, labels=date_labels, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); |
289 axis(side=2); | 320 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); |
290 if (plot_std_error=="yes") { | 321 if (plot_std_error=="yes") { |
291 # Standard error for group. | 322 # Standard error for group. |
292 lines(days, group+group_std_error, lty=2); | 323 lines(days, group+group_std_error, lty=2); |
293 lines(days, group-group_std_error, lty=2); | 324 lines(days, group-group_std_error, lty=2); |
294 # Standard error for group2. | 325 # Standard error for group2. |
312 stage = paste(life_stages_adult, "Adult Pop", sep=" "); | 343 stage = paste(life_stages_adult, "Adult Pop", sep=" "); |
313 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); | 344 title = paste(insect, ": Reps", replications, ":", stage, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); |
314 legend_text = c(paste(life_stages_adult, life_stage, sep=" ")); | 345 legend_text = c(paste(life_stages_adult, life_stage, sep=" ")); |
315 columns = c(1); | 346 columns = c(1); |
316 } | 347 } |
317 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3); | 348 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); |
318 legend("topleft", legend_text, lty=c(1), col="black", cex=3); | 349 legend("topleft", legend_text, lty=c(1), col="black", cex=3); |
319 axis(side=1, at=ticks, labels=date_labels); | 350 axis(side=1, at=ticks, labels=date_labels, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); |
320 axis(side=2); | 351 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); |
321 if (plot_std_error=="yes") { | 352 if (plot_std_error=="yes") { |
322 # Standard error for group. | 353 # Standard error for group. |
323 lines(days, group+group_std_error, lty=2); | 354 lines(days, group+group_std_error, lty=2); |
324 lines(days, group-group_std_error, lty=2); | 355 lines(days, group-group_std_error, lty=2); |
325 } | 356 } |
335 title_str = paste(":", life_stages_adult, "Adult Pop by Gen", ":", sep=" "); | 366 title_str = paste(":", life_stages_adult, "Adult Pop by Gen", ":", sep=" "); |
336 } | 367 } |
337 title = paste(insect, ": Reps", replications, title_str, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); | 368 title = paste(insect, ": Reps", replications, title_str, location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" "); |
338 legend_text = c("P", "F1", "F2"); | 369 legend_text = c("P", "F1", "F2"); |
339 columns = c(1, 2, 4); | 370 columns = c(1, 2, 4); |
340 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3); | 371 plot(days, group, main=title, type="l", ylim=c(0, maxval), axes=FALSE, lwd=2, xlab="", ylab=""); |
341 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3); | 372 legend("topleft", legend_text, lty=c(1, 1, 1), col=columns, cex=3); |
342 lines(days, group2, lwd=2, lty=1, col=2); | 373 lines(days, group2, lwd=2, lty=1, col=2); |
343 lines(days, group3, lwd=2, lty=1, col=4); | 374 lines(days, group3, lwd=2, lty=1, col=4); |
344 axis(side=1, at=ticks, labels=date_labels); | 375 axis(side=1, at=ticks, labels=date_labels, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); |
345 axis(side=2); | 376 axis(side=2, font.axis=3, xpd=TRUE, cex=3, cex.lab=3, cex.axis=3, cex.main=3); |
346 if (plot_std_error=="yes") { | 377 if (plot_std_error=="yes") { |
347 # Standard error for group. | 378 # Standard error for group. |
348 lines(days, group+group_std_error, lty=2); | 379 lines(days, group+group_std_error, lty=2); |
349 lines(days, group-group_std_error, lty=2); | 380 lines(days, group-group_std_error, lty=2); |
350 # Standard error for group2. | 381 # Standard error for group2. |