Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 4:e7b1fc0133bb draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Mon, 13 Nov 2017 12:57:46 -0500
parents 24fa0d35a8bf
children 1878a03f9c9f
comparison
equal deleted inserted replaced
3:24fa0d35a8bf 4:e7b1fc0133bb
23 23
24 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 24 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
25 args <- parse_args(parser, positional_arguments=TRUE) 25 args <- parse_args(parser, positional_arguments=TRUE)
26 opt <- args$options 26 opt <- args$options
27 27
28 convert_csv_to_rdata=function(temperature_data, data_matrix) 28 get_daylight_length = function(latitude, temperature_data, num_days)
29 { 29 {
30 # Integer day of the year. 30 # Return a vector of daylight length (photoperido profile) for
31 data_matrix[,1] <- c(1:opt$num_days) 31 # the number of days specified in the input temperature data
32 # Minimum 32 # (from Forsythe 1995).
33 data_matrix[,2] <- temperature_data[c(1:opt$num_days), 5] 33 p = 0.8333
34 # Maximum 34 daylight_length_vector <- NULL
35 data_matrix[,3] <- temperature_data[c(1:opt$num_days), 6]
36 namedat <- "tempdata.Rdat"
37 save(data_matrix, file=namedat)
38 namedat
39 }
40
41 daylength=function(latitude, num_days)
42 {
43 # From Forsythe 1995.
44 p=0.8333
45 dl <- NULL
46 for (i in 1:num_days) { 35 for (i in 1:num_days) {
47 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186))) 36 # Get the day of the year from the current row
37 # of the temperature data for computation.
38 doy <- temperature_data[i, 4]
39 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)))
48 phi <- asin(0.39795 * cos(theta)) 40 phi <- asin(0.39795 * cos(theta))
49 dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi))) 41 # Compute the length of daylight for the day of the year.
50 } 42 daylight_length_vector[i] <- 24 - (24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi))))
51 # Return a vector of daylength for the number of 43 }
52 # days specified in the input temperature data. 44 daylight_length_vector
53 dl 45 }
54 } 46
55 47 get_temperature_at_hour = function(latitude, temperature_data, daylight_length_vector, row, num_days)
56 hourtemp=function(latitude, date, temperature_file_path, num_days) 48 {
57 {
58 load(temperature_file_path)
59 # Base development threshold for Brown Marmolated Stink Bug 49 # Base development threshold for Brown Marmolated Stink Bug
60 # insect phenology model. 50 # insect phenology model.
51 # TODO: Pass insect on the command line to accomodate more
52 # the just the Brown Marmolated Stink Bub.
61 threshold <- 14.17 53 threshold <- 14.17
62 dnp <- data_matrix[date, 2] # daily minimum 54
63 dxp <- data_matrix[date, 3] # daily maximum 55 # Input temperature currently has the following columns.
56 # # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
57 # Minimum temperature for current row.
58 dnp <- temperature_data[row, 5]
59 # Maximum temperature for current row.
60 dxp <- temperature_data[row, 6]
61 # Mean temperature for current row.
64 dmean <- 0.5 * (dnp + dxp) 62 dmean <- 0.5 * (dnp + dxp)
65 dd <- 0 # initialize degree day accumulation 63 # Initialize degree day accumulation
66 64 dd <- 0
67 if (dxp<threshold) { 65 if (dxp < threshold) {
68 dd <- 0 66 dd <- 0
69 } 67 }
70 else { 68 else {
71 # Extract daylength data for the number of
72 # days specified in the input temperature data.
73 dlprofile <- daylength(latitude, num_days)
74 # Initialize hourly temperature. 69 # Initialize hourly temperature.
75 T <- NULL 70 T <- NULL
76 # Initialize degree hour vector. 71 # Initialize degree hour vector.
77 dh <- NULL 72 dh <- NULL
78 # Calculate daylength in given date. 73 # Daylight length for current row.
79 y <- dlprofile[date] 74 y <- daylight_length_vector[row]
80 # Night length. 75 # Darkness length.
81 z <- 24 - y 76 z <- 24 - y
82 # Lag coefficient. 77 # Lag coefficient.
83 a <- 1.86 78 a <- 1.86
84 # Night coefficient. 79 # Darkness coefficient.
85 b <- 2.20 80 b <- 2.20
86 # Sunrise time. 81 # Sunrise time.
87 risetime <- 12 - y / 2 82 risetime <- 12 - y / 2
88 # Sunset time. 83 # Sunset time.
89 settime <- 12 + y / 2 84 settime <- 12 + y / 2
90 ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp 85 ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp
91 for (i in 1:24) { 86 for (i in 1:24) {
92 if (i > risetime && i<settime) { 87 if (i > risetime && i < settime) {
93 # Number of hours after Tmin until sunset. 88 # Number of hours after Tmin until sunset.
94 m <- i - 5 89 m <- i - 5
95 T[i]=(dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp 90 T[i] = (dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp
96 if (T[i]<8.4) { 91 if (T[i] < 8.4) {
97 dh[i] <- 0 92 dh[i] <- 0
98 } 93 }
99 else { 94 else {
100 dh[i] <- T[i] - 8.4 95 dh[i] <- T[i] - 8.4
101 } 96 }
102 } 97 }
103 else if (i > settime) { 98 else if (i > settime) {
104 n <- i - settime 99 n <- i - settime
105 T[i]=dnp + (ts - dnp) * exp( - b * n / z) 100 T[i] = dnp + (ts - dnp) * exp( - b * n / z)
106 if (T[i]<8.4) { 101 if (T[i] < 8.4) {
107 dh[i] <- 0 102 dh[i] <- 0
108 } 103 }
109 else { 104 else {
110 dh[i] <- T[i] - 8.4 105 dh[i] <- T[i] - 8.4
111 } 106 }
112 } 107 }
113 else { 108 else {
114 n <- i + 24 - settime 109 n <- i + 24 - settime
115 T[i]=dnp + (ts - dnp) * exp( - b * n / z) 110 T[i]=dnp + (ts - dnp) * exp( - b * n / z)
116 if (T[i]<8.4) { 111 if (T[i] < 8.4) {
117 dh[i] <- 0 112 dh[i] <- 0
118 } 113 }
119 else { 114 else {
120 dh[i] <- T[i] - 8.4 115 dh[i] <- T[i] - 8.4
121 } 116 }
197 return = mort.prob 192 return = mort.prob
198 return 193 return
199 } 194 }
200 195
201 # Read in the input temperature datafile into a Data Frame object. 196 # Read in the input temperature datafile into a Data Frame object.
202 temperature_data <- read.csv(file=opt$input, header=T, sep=",") 197 # The input data currently must have 6 columns:
203 start_date <- temperature_data[c(1:1), 3] 198 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
204 end_date <- temperature_data[c(opt$num_days:opt$num_days), 3] 199 temperature_data <- read.csv(file=opt$input, header=T, strip.white=TRUE, sep=",")
205 raw_data_matrix <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days) 200 start_date <- temperature_data[1, 3]
206 temperature_file_path <- convert_csv_to_rdata(temperature_data, raw_data_matrix) 201 end_date <- temperature_data[opt$num_days, 3]
207 latitude <- temperature_data[1, 1] 202 latitude <- temperature_data[1, 1]
203 daylight_length_vector <- get_daylight_length(latitude, temperature_data, opt$num_days)
208 204
209 cat("Number of days: ", opt$num_days, "\n") 205 cat("Number of days: ", opt$num_days, "\n")
210 206
211 # Initialize matrix for results from all replications. 207 # Initialize matrix for results from all replications.
212 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications) 208 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications)
221 vec.ini <- c(0, 3, 0, 0, 0) 217 vec.ini <- c(0, 3, 0, 0, 0)
222 # Overwintering, previttelogenic, DD=0, T=0, no-diapause. 218 # Overwintering, previttelogenic, DD=0, T=0, no-diapause.
223 vec.mat <- rep(vec.ini, n) 219 vec.mat <- rep(vec.ini, n)
224 # Complete matrix for the population. 220 # Complete matrix for the population.
225 vec.mat <- base::t(matrix(vec.mat, nrow=5)) 221 vec.mat <- base::t(matrix(vec.mat, nrow=5))
226 # Complete photoperiod profile in a year, requires daylength function.
227 ph.p <- daylength(latitude, opt$num_days)
228
229 # Time series of population size. 222 # Time series of population size.
230 tot.pop <- NULL 223 tot.pop <- NULL
231 gen0.pop <- rep(0, opt$num_days) 224 gen0.pop <- rep(0, opt$num_days)
232 gen1.pop <- rep(0, opt$num_days) 225 gen1.pop <- rep(0, opt$num_days)
233 gen2.pop <- rep(0, opt$num_days) 226 gen2.pop <- rep(0, opt$num_days)
235 g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days) 228 g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days)
236 N.newborn <- N.death <- N.adult <- rep(0, opt$num_days) 229 N.newborn <- N.death <- N.adult <- rep(0, opt$num_days)
237 dd.day <- rep(0, opt$num_days) 230 dd.day <- rep(0, opt$num_days)
238 231
239 # All the days included in the input temperature dataset. 232 # All the days included in the input temperature dataset.
240 for (day in 1:opt$num_days) { 233 for (row in 1:opt$num_days) {
234 # Get the integer day of the year for the current row.
235 doy <- temperature_data[row, 4]
241 # Photoperiod in the day. 236 # Photoperiod in the day.
242 photoperiod <- ph.p[day] 237 photoperiod <- daylight_length_vector[row]
243 temp.profile <- hourtemp(latitude, day, temperature_file_path, opt$num_days) 238 temp.profile <- get_temperature_at_hour(latitude, temperature_data, daylight_length_vector, row, opt$num_days)
244 mean.temp <- temp.profile[1] 239 mean.temp <- temp.profile[1]
245 dd.temp <- temp.profile[2] 240 dd.temp <- temp.profile[2]
246 dd.day[day] <- dd.temp 241 dd.day[row] <- dd.temp
247 # Trash bin for death. 242 # Trash bin for death.
248 death.vec <- NULL 243 death.vec <- NULL
249 # Newborn. 244 # Newborn.
250 birth.vec <- NULL 245 birth.vec <- NULL
251 246
270 else if (vec.ind[2] == 1 | vec.ind[2] == 2) { 265 else if (vec.ind[2] == 1 | vec.ind[2] == 2) {
271 death.prob = opt$nymph_mort * mortality.nymph(mean.temp) 266 death.prob = opt$nymph_mort * mortality.nymph(mean.temp)
272 } 267 }
273 else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) { 268 else if (vec.ind[2] == 3 | vec.ind[2] == 4 | vec.ind[2] == 5) {
274 # For adult. 269 # For adult.
275 if (day < day.kill) { 270 if (doy < day.kill) {
276 death.prob = opt$adult_mort * mortality.adult(mean.temp) 271 death.prob = opt$adult_mort * mortality.adult(mean.temp)
277 } 272 }
278 else { 273 else {
279 # Increase adult mortality after fall equinox. 274 # Increase adult mortality after fall equinox.
280 death.prob = opt$adult_mort * post.mort * mortality.adult(mean.temp) 275 death.prob = opt$adult_mort * post.mort * mortality.adult(mean.temp)
288 else { 283 else {
289 # Aggregrate index of dead bug. 284 # Aggregrate index of dead bug.
290 # Event 1 end of diapause. 285 # Event 1 end of diapause.
291 if (vec.ind[1] == 0 && vec.ind[2] == 3) { 286 if (vec.ind[1] == 0 && vec.ind[2] == 3) {
292 # Overwintering adult (previttelogenic). 287 # Overwintering adult (previttelogenic).
293 if (photoperiod > opt$photoperiod && vec.ind[3] > 68 && day < 180) { 288 if (photoperiod > opt$photoperiod && vec.ind[3] > 68 && doy < 180) {
294 # Add 68C to become fully reproductively matured. 289 # Add 68C to become fully reproductively matured.
295 # Transfer to vittelogenic. 290 # Transfer to vittelogenic.
296 vec.ind <- c(0, 4, 0, 0, 0) 291 vec.ind <- c(0, 4, 0, 0, 0)
297 vec.mat[i,] <- vec.ind 292 vec.mat[i,] <- vec.ind
298 } 293 }
354 birth.vec <- rbind(birth.vec, new.vec) 349 birth.vec <- rbind(birth.vec, new.vec)
355 } 350 }
356 } 351 }
357 352
358 # Event 2 oviposition -- for gen 1. 353 # Event 2 oviposition -- for gen 1.
359 if (vec.ind[2] == 4 && vec.ind[1] == 1 && mean.temp > 12.5 && day < 222) { 354 if (vec.ind[2] == 4 && vec.ind[1] == 1 && mean.temp > 12.5 && doy < 222) {
360 # Vittelogenic stage, 1st generation 355 # Vittelogenic stage, 1st generation
361 if (vec.ind[4] == 0) { 356 if (vec.ind[4] == 0) {
362 # Just turned in vittelogenic stage. 357 # Just turned in vittelogenic stage.
363 n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size)) 358 n.birth=round(runif(1, 2 + opt$min_clutch_size, 8 + opt$max_clutch_size))
364 } 359 }
415 if (vec.ind[3] >= (250 + opt$old_nymph_accum)) { 410 if (vec.ind[3] >= (250 + opt$old_nymph_accum)) {
416 # From young to old nymph, dd requirement met. 411 # From young to old nymph, dd requirement met.
417 current.gen <- vec.ind[1] 412 current.gen <- vec.ind[1]
418 # Transfer to old nym stage. 413 # Transfer to old nym stage.
419 vec.ind <- c(current.gen, 2, 0, 0, 0) 414 vec.ind <- c(current.gen, 2, 0, 0, 0)
420 if (photoperiod < opt$photoperiod && day > 180) { 415 if (photoperiod < opt$photoperiod && doy > 180) {
421 vec.ind[5] <- 1 416 vec.ind[5] <- 1
422 } # Prepare for diapausing. 417 } # Prepare for diapausing.
423 } 418 }
424 else { 419 else {
425 # Add 1 day in current stage. 420 # Add 1 day in current stage.
493 gen1 <- sum(vec.mat[,1] == 1) 488 gen1 <- sum(vec.mat[,1] == 1)
494 # Second generation. 489 # Second generation.
495 gen2 <- sum(vec.mat[,1] == 2) 490 gen2 <- sum(vec.mat[,1] == 2)
496 # Sum of all adults. 491 # Sum of all adults.
497 n.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5) 492 n.adult <- sum(vec.mat[,2] == 3) + sum(vec.mat[,2] == 4) + sum(vec.mat[,2] == 5)
498 # Gen eration 0 pop size. 493
499 gen0.pop[day] <- gen0 494 # Generation 0 pop size.
500 gen1.pop[day] <- gen1 495 gen0.pop[row] <- gen0
501 gen2.pop[day] <- gen2 496 gen1.pop[row] <- gen1
502 S0[day] <- s0 497 gen2.pop[row] <- gen2
503 S1[day] <- s1 498
504 S2[day] <- s2 499 S0[row] <- s0
505 S3[day] <- s3 500 S1[row] <- s1
506 S4[day] <- s4 501 S2[row] <- s2
507 S5[day] <- s5 502 S3[row] <- s3
508 g0.adult[day] <- sum(vec.mat[,1] == 0) 503 S4[row] <- s4
509 g1.adult[day] <- sum((vec.mat[,1] == 1 & vec.mat[,2] == 3) | (vec.mat[,1] == 1 & vec.mat[,2] == 4) | (vec.mat[,1] == 1 & vec.mat[,2] == 5)) 504 S5[row] <- s5
510 g2.adult[day] <- sum((vec.mat[,1]== 2 & vec.mat[,2] == 3) | (vec.mat[,1] == 2 & vec.mat[,2] == 4) | (vec.mat[,1] == 2 & vec.mat[,2] == 5)) 505
511 506 g0.adult[row] <- sum(vec.mat[,1] == 0)
512 N.newborn[day] <- n.newborn 507 g1.adult[row] <- sum((vec.mat[,1] == 1 & vec.mat[,2] == 3) | (vec.mat[,1] == 1 & vec.mat[,2] == 4) | (vec.mat[,1] == 1 & vec.mat[,2] == 5))
513 N.death[day] <- n.death 508 g2.adult[row] <- sum((vec.mat[,1]== 2 & vec.mat[,2] == 3) | (vec.mat[,1] == 2 & vec.mat[,2] == 4) | (vec.mat[,1] == 2 & vec.mat[,2] == 5))
514 N.adult[day] <- n.adult 509
510 N.newborn[row] <- n.newborn
511 N.death[row] <- n.death
512 N.adult[row] <- n.adult
515 } # end of days specified in the input temperature data 513 } # end of days specified in the input temperature data
516 514
517 dd.cum <- cumsum(dd.day) 515 dd.cum <- cumsum(dd.day)
516
518 # Collect all the outputs. 517 # Collect all the outputs.
519 S0.rep[,N.rep] <- S0 518 S0.rep[,N.rep] <- S0
520 S1.rep[,N.rep] <- S1 519 S1.rep[,N.rep] <- S1
521 S2.rep[,N.rep] <- S2 520 S2.rep[,N.rep] <- S2
522 S3.rep[,N.rep] <- S3 521 S3.rep[,N.rep] <- S3
532 g0a.rep[,N.rep] <- g0.adult 531 g0a.rep[,N.rep] <- g0.adult
533 g1a.rep[,N.rep] <- g1.adult 532 g1a.rep[,N.rep] <- g1.adult
534 g2a.rep[,N.rep] <- g2.adult 533 g2a.rep[,N.rep] <- g2.adult
535 } 534 }
536 535
537 # Data analysis and visualization 536 # Data analysis and visualization can currently
538 # default: plot 1 year of result 537 # plot only within a single calendar year.
539 # but can be expanded to accommodate multiple years 538 # TODO: enhance this to accomodate multiple calendar years.
540 n.yr <- 1 539 n.yr <- 1
541 day.all <- c(1:opt$num_days * n.yr) 540 day.all <- c(1:opt$num_days * n.yr)
542 541
543 # mean value for adults 542 # mean value for adults
544 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean) 543 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean)
576 # SE for F1 adult 575 # SE for F1 adult
577 g1a.se <- apply(g1a.rep, 1, sd) / sqrt(opt$replications) 576 g1a.se <- apply(g1a.rep, 1, sd) / sqrt(opt$replications)
578 # SE for F2 adult 577 # SE for F2 adult
579 g2a.se <- apply(g2a.rep, 1, sd) / sqrt(opt$replications) 578 g2a.se <- apply(g2a.rep, 1, sd) / sqrt(opt$replications)
580 579
581 dev.new(width=20, height=20) 580 dev.new(width=20, height=30)
582 581
583 # Start PDF device driver to save charts to output. 582 # Start PDF device driver to save charts to output.
584 pdf(file=opt$output, height=20, width=20, bg="white") 583 pdf(file=opt$output, width=20, height=30, bg="white")
585 584
586 par(mar = c(5, 6, 4, 4), mfrow=c(3, 1)) 585 par(mar = c(5, 6, 4, 4), mfrow=c(3, 1))
587 586
588 # Subfigure 2: population size by life stage 587 # Subfigure 1: population size by life stage
589 title <- paste("BSMB Total Population Size by Life Stage:", opt$location, ", Latitude:", latitude, ", Temperature Dates:", start_date, "to", end_date, sep=" ") 588 title <- paste("BSMB total population by life stage :", opt$location, ": Lat:", latitude, ":", start_date, "to", end_date, sep=" ")
590 plot(day.all, sa, main=title, type="l", ylim=c(0, max(se + se.se, sn + sn.se, sa + sa.se)), axes=F, lwd=2, xlab="", ylab="Number", cex=2, cex.lab=2, cex.axis=2, cex.main=2) 589 plot(day.all, sa, main=title, type="l", ylim=c(0, max(se + se.se, sn + sn.se, sa + sa.se)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
591 # Young and old nymphs. 590 # Young and old nymphs.
592 lines(day.all, sn, lwd=2, lty=1, col=2) 591 lines(day.all, sn, lwd=2, lty=1, col=2)
593 # Eggs 592 # Eggs
594 lines(day.all, se, lwd=2, lty=1, col=4) 593 lines(day.all, se, lwd=2, lty=1, col=4)
595 axis(1, at = c(1:12) * 30 - 15, cex.axis=2, labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) 594 axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels=c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
596 axis(2, cex.axis = 2) 595 axis(2, cex.axis=3)
597 leg.text <- c("Egg", "Nymph", "Adult") 596 leg.text <- c("Egg", "Nymph", "Adult")
598 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=2) 597 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=3)
599 if (opt$se_plot == 1) { 598 if (opt$se_plot == 1) {
600 # add SE lines to plot 599 # Add SE lines to plot
601 # SE for adults 600 # SE for adults
602 lines (day.all, sa + sa.se, lty=2) 601 lines (day.all, sa + sa.se, lty=2)
603 lines (day.all, sa - sa.se, lty=2) 602 lines (day.all, sa - sa.se, lty=2)
604 # SE for nymphs 603 # SE for nymphs
605 lines (day.all, sn + sn.se, col=2, lty=2) 604 lines (day.all, sn + sn.se, col=2, lty=2)
606 lines (day.all, sn - sn.se, col=2, lty=2) 605 lines (day.all, sn - sn.se, col=2, lty=2)
607 # SE for eggs 606 # SE for eggs
608 lines (day.all, se + se.se, col=4, lty=2) 607 lines (day.all, se + se.se, col=4, lty=2)
609 lines (day.all, se - se.se, col=4, lty=2) 608 lines (day.all, se - se.se, col=4, lty=2)
610 } 609 }
611 610
612 # Subfigure 3: population size by generation 611 # Subfigure 2: population size by generation
613 title <- paste("BSMB Total Population Size by Generation:", opt$location, ", Latitude:", latitude, ", Temperature Dates:", start_date, "to", end_date, sep=" ") 612 title <- paste("BSMB total population by generation :", opt$location, ": Lat:", latitude, ":", start_date, "to", end_date, sep=" ")
614 plot(day.all, g0, main=title, type="l", ylim=c(0, max(g2)), axes=F, lwd=2, xlab="", ylab="Number", cex=2, cex.lab=2, cex.axis=2, cex.main=2) 613 plot(day.all, g0, main=title, type="l", ylim=c(0, max(g2)), axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
615 lines(day.all, g1, lwd = 2, lty = 1, col = 2) 614 lines(day.all, g1, lwd = 2, lty = 1, col=2)
616 lines(day.all, g2, lwd = 2, lty = 1, col = 4) 615 lines(day.all, g2, lwd = 2, lty = 1, col=4)
617 axis(1, at = c(1:12) * 30 - 15, cex.axis = 2, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) 616 axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
618 axis(2, cex.axis = 2) 617 axis(2, cex.axis=3)
619 leg.text <- c("P", "F1", "F2") 618 leg.text <- c("P", "F1", "F2")
620 legend("topleft", leg.text, lty = c(1, 1, 1), col =c(1, 2, 4), cex = 2) 619 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
621 if (opt$se_plot == 1) { 620 if (opt$se_plot == 1) {
622 # add SE lines to plot 621 # Add SE lines to plot
623 # SE for adults 622 # SE for adults
624 lines (day.all, g0 + g0.se, lty = 2) 623 lines (day.all, g0+g0.se, lty=2)
625 lines (day.all, g0 - g0.se, lty = 2) 624 lines (day.all, g0-g0.se, lty=2)
626 # SE for nymphs 625 # SE for nymphs
627 lines (day.all, g1 + g1.se, col = 2, lty = 2) 626 lines (day.all, g1+g1.se, col=2, lty=2)
628 lines (day.all, g1 - g1.se, col = 2, lty = 2) 627 lines (day.all, g1-g1.se, col=2, lty=2)
629 # SE for eggs 628 # SE for eggs
630 lines (day.all, g2 + g2.se, col = 4, lty = 2) 629 lines (day.all, g2+g2.se, col=4, lty=2)
631 lines (day.all, g2 - g2.se, col = 4, lty = 2) 630 lines (day.all, g2-g2.se, col=4, lty=2)
632 } 631 }
633 632
634 # Subfigure 4: adult population size by generation 633 # Subfigure 3: adult population size by generation
635 title <- paste("BSMB Adult Population Size by Generation:", opt$location, ", Latitude:", latitude, ", Temperature Dates:", start_date, "to", end_date, sep=" ") 634 title <- paste("BSMB adult population by generation :", opt$location, ": Lat:", latitude, ":", start_date, "to", end_date, sep=" ")
636 plot(day.all, g0a, ylim=c(0, max(g2a) + 100), main=title, type="l", axes=F, lwd=2, xlab="Year", ylab="Number", cex=2, cex.lab=2, cex.axis=2, cex.main=2) 635 plot(day.all, g0a, ylim=c(0, max(g2a) + 100), main=title, type="l", axes=F, lwd=2, xlab="", ylab="", cex=3, cex.lab=3, cex.axis=3, cex.main=3)
637 lines(day.all, g1a, lwd = 2, lty = 1, col = 2) 636 lines(day.all, g1a, lwd = 2, lty = 1, col=2)
638 lines(day.all, g2a, lwd = 2, lty = 1, col = 4) 637 lines(day.all, g2a, lwd = 2, lty = 1, col=4)
639 axis(1, at = c(1:12) * 30 - 15, cex.axis = 2, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")) 638 axis(1, at=c(1:12) * 30 - 15, cex.axis=3, labels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))
640 axis(2, cex.axis = 2) 639 axis(2, cex.axis=3)
641 leg.text <- c("P", "F1", "F2") 640 leg.text <- c("P", "F1", "F2")
642 legend("topleft", leg.text, lty = c(1, 1, 1), col = c(1, 2, 4), cex = 2) 641 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(1, 2, 4), cex=3)
643 if (opt$se_plot == 1) { 642 if (opt$se_plot == 1) {
644 # add SE lines to plot 643 # Add SE lines to plot
645 # SE for adults 644 # SE for adults
646 lines (day.all, g0a + g0a.se, lty = 2) 645 lines (day.all, g0a+g0a.se, lty=2)
647 lines (day.all, g0a - g0a.se, lty = 2) 646 lines (day.all, g0a-g0a.se, lty=2)
648 # SE for nymphs 647 # SE for nymphs
649 lines (day.all, g1a + g1a.se, col = 2, lty = 2) 648 lines (day.all, g1a+g1a.se, col=2, lty=2)
650 lines (day.all, g1a - g1a.se, col = 2, lty = 2) 649 lines (day.all, g1a-g1a.se, col=2, lty=2)
651 # SE for eggs 650 # SE for eggs
652 lines (day.all, g2a + g2a.se, col = 4, lty = 2) 651 lines (day.all, g2a+g2a.se, col=4, lty=2)
653 lines (day.all, g2a - g2a.se, col = 4, lty = 2) 652 lines (day.all, g2a-g2a.se, col=4, lty=2)
654 } 653 }
655 654
656 # Turn off device driver to flush output. 655 # Turn off device driver to flush output.
657 dev.off() 656 dev.off()