Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 31:e56529d02c3e draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Thu, 09 Nov 2017 13:25:06 -0500
parents 9146801c5ba1
children 5517d380c337
comparison
equal deleted inserted replaced
30:9146801c5ba1 31:e56529d02c3e
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(loc, temperature_data) 28 convert_csv_to_rdata=function(data_matrix)
29 { 29 {
30 expdata <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days) 30 # Also stores the current year.
31 expdata[,1] <- c(1:opt$num_days) 31 date <- data_matrix[1, 3]
32 year <- substr(date, 1, 4)
33 # Integer day of the year.
34 data_matrix[,1] <- c(1:opt$num_days)
32 # Minimum 35 # Minimum
33 expdata[,2] <- temperature_data[c(1:opt$num_days), 5] 36 data_matrix[,2] <- temperature_data[c(1:opt$num_days), 5]
34 # Maximum 37 # Maximum
35 expdata[,3] <- temperature_data[c(1:opt$num_days), 6] 38 data_matrix[,3] <- temperature_data[c(1:opt$num_days), 6]
36 namedat <- "tempdata.Rdat" 39 namedat <- "tempdata.Rdat"
37 save(expdata, file=namedat) 40 save(data_matrix, file=namedat)
38 namedat 41 namedat
39 } 42 }
40 43
41 daylength=function(latitude, num_days) 44 daylength=function(latitude, num_days)
42 { 45 {
57 { 60 {
58 load(temperature_file_path) 61 load(temperature_file_path)
59 # Base development threshold for Brown Marmolated Stink Bug 62 # Base development threshold for Brown Marmolated Stink Bug
60 # insect phenology model. 63 # insect phenology model.
61 threshold <- 14.17 64 threshold <- 14.17
62 dnp <- expdata[date, 2] # daily minimum 65 dnp <- data_matrix[date, 2] # daily minimum
63 dxp <- expdata[date, 3] # daily maximum 66 dxp <- data_matrix[date, 3] # daily maximum
64 dmean <- 0.5 * (dnp + dxp) 67 dmean <- 0.5 * (dnp + dxp)
65 dd <- 0 # initialize degree day accumulation 68 dd <- 0 # initialize degree day accumulation
66 69
67 if (dxp<threshold) { 70 if (dxp<threshold) {
68 dd <- 0 71 dd <- 0
198 return 201 return
199 } 202 }
200 203
201 # Read in the input temperature datafile into a Data Frame object. 204 # Read in the input temperature datafile into a Data Frame object.
202 temperature_data <- read.csv(file=opt$input, header=T, sep=",") 205 temperature_data <- read.csv(file=opt$input, header=T, sep=",")
203 temperature_file_path <- convert_csv_to_rdata(opt$location, temperature_data) 206 raw_data_matrix <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days)
207 temperature_file_path <- convert_csv_to_rdata(raw_data_matrix)
204 latitude <- temperature_data[1, 1] 208 latitude <- temperature_data[1, 1]
205 209
206 cat("Number of days: ", opt$num_days, "\n") 210 cat("Number of days: ", opt$num_days, "\n")
207 cat("Latitude: ", latitude, "\n") 211 cat("Latitude: ", latitude, "\n")
208 cat("Replications: ", opt$replications, "\n") 212 cat("Replications: ", opt$replications, "\n")
576 g1a.se <- apply(g1a.rep, 1, sd) / sqrt(opt$replications) 580 g1a.se <- apply(g1a.rep, 1, sd) / sqrt(opt$replications)
577 # SE for F2 adult 581 # SE for F2 adult
578 g2a.se <- apply(g2a.rep, 1, sd) / sqrt(opt$replications) 582 g2a.se <- apply(g2a.rep, 1, sd) / sqrt(opt$replications)
579 583
580 dev.new(width=20, height=20) 584 dev.new(width=20, height=20)
581 sink("/dev/null")
582 585
583 # Start PDF device driver to save charts to output. 586 # Start PDF device driver to save charts to output.
584 pdf(file=opt$output, height=20, width=20, bg="white") 587 pdf(file=opt$output, height=20, width=20, bg="white")
585 588
586 par(mar = c(5, 6, 4, 4), mfrow=c(3, 1)) 589 par(mar = c(5, 6, 4, 4), mfrow=c(3, 1))
587 590
588 # Subfigure 2: population size by life stage 591 # Subfigure 2: population size by life stage
589 plot(day.all, sa, main = "BSMB Total Population Size by Life Stage", 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) 592 title <- paste(opt$location, year, "BSMB Total Population Size by Life Stage", sep=" ")
590 # Young and old nymphs 593 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)
591 lines(day.all, sn, lwd = 2, lty = 1, col = 2) 594 # Young and old nymphs.
595 lines(day.all, sn, lwd=2, lty=1, col=2)
592 # Eggs 596 # Eggs
593 lines(day.all, se, lwd = 2, lty = 1, col = 4) 597 lines(day.all, se, lwd=2, lty=1, col=4)
594 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")) 598 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"))
595 axis(2, cex.axis = 2) 599 axis(2, cex.axis = 2)
596 leg.text <- c("Egg", "Nymph", "Adult") 600 leg.text <- c("Egg", "Nymph", "Adult")
597 legend("topleft", leg.text, lty = c(1, 1, 1), col = c(4, 2, 1), cex = 2) 601 legend("topleft", leg.text, lty=c(1, 1, 1), col=c(4, 2, 1), cex=2)
598 if (opt$se_plot == 1) { 602 if (opt$se_plot == 1) {
599 # add SE lines to plot 603 # add SE lines to plot
600 # SE for adults 604 # SE for adults
601 lines (day.all, sa + sa.se, lty = 2) 605 lines (day.all, sa + sa.se, lty=2)
602 lines (day.all, sa - sa.se, lty = 2) 606 lines (day.all, sa - sa.se, lty=2)
603 # SE for nymphs 607 # SE for nymphs
604 lines (day.all, sn + sn.se, col = 2, lty = 2) 608 lines (day.all, sn + sn.se, col=2, lty=2)
605 lines (day.all, sn - sn.se, col = 2, lty = 2) 609 lines (day.all, sn - sn.se, col=2, lty=2)
606 # SE for eggs 610 # SE for eggs
607 lines (day.all, se + se.se, col = 4, lty = 2) 611 lines (day.all, se + se.se, col=4, lty=2)
608 lines (day.all, se - se.se, col = 4, lty = 2) 612 lines (day.all, se - se.se, col=4, lty=2)
609 } 613 }
610 614
611 # Subfigure 3: population size by generation 615 # Subfigure 3: population size by generation
612 plot(day.all, g0, main = "BSMB Total Population Size by Generation", 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) 616 title <- paste(opt$location, year, "BSMB Total Population Size by Generation", sep=" ")
617 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 lines(day.all, g1, lwd = 2, lty = 1, col = 2) 618 lines(day.all, g1, lwd = 2, lty = 1, col = 2)
614 lines(day.all, g2, lwd = 2, lty = 1, col = 4) 619 lines(day.all, g2, lwd = 2, lty = 1, col = 4)
615 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")) 620 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(2, cex.axis = 2) 621 axis(2, cex.axis = 2)
617 leg.text <- c("P", "F1", "F2") 622 leg.text <- c("P", "F1", "F2")
628 lines (day.all, g2 + g2.se, col = 4, lty = 2) 633 lines (day.all, g2 + g2.se, col = 4, lty = 2)
629 lines (day.all, g2 - g2.se, col = 4, lty = 2) 634 lines (day.all, g2 - g2.se, col = 4, lty = 2)
630 } 635 }
631 636
632 # Subfigure 4: adult population size by generation 637 # Subfigure 4: adult population size by generation
633 plot(day.all, g0a, ylim = c(0, max(g2a) + 100), main = "BSMB Adult Population Size by Generation", type = "l", axes = F, lwd = 2, xlab = "Year", ylab = "Number", cex = 2, cex.lab = 2, cex.axis = 2, cex.main = 2) 638 title <- paste(opt$location, year, "BSMB Adult Population Size by Generation", sep=" ")
639 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)
634 lines(day.all, g1a, lwd = 2, lty = 1, col = 2) 640 lines(day.all, g1a, lwd = 2, lty = 1, col = 2)
635 lines(day.all, g2a, lwd = 2, lty = 1, col = 4) 641 lines(day.all, g2a, lwd = 2, lty = 1, col = 4)
636 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")) 642 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"))
637 axis(2, cex.axis = 2) 643 axis(2, cex.axis = 2)
638 leg.text <- c("P", "F1", "F2") 644 leg.text <- c("P", "F1", "F2")
648 # SE for eggs 654 # SE for eggs
649 lines (day.all, g2a + g2a.se, col = 4, lty = 2) 655 lines (day.all, g2a + g2a.se, col = 4, lty = 2)
650 lines (day.all, g2a - g2a.se, col = 4, lty = 2) 656 lines (day.all, g2a - g2a.se, col = 4, lty = 2)
651 } 657 }
652 658
653 sink()
654 # Turn off device driver to flush output. 659 # Turn off device driver to flush output.
655 dev.off() 660 dev.off()