Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 44:5260db6479db draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Sun, 12 Nov 2017 10:19:10 -0500
parents 758c6e5f29a9
children d7e406de8094
comparison
equal deleted inserted replaced
43:1b3020d54cbe 44:5260db6479db
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 convert_csv_to_rdata=function(start_doy, end_doy, temperature_data, data_matrix)
29 { 29 {
30 # Integer day of the year. 30 # Make sure the first column starts with the start date in
31 data_matrix[,1] <- c(1:opt$num_days) 31 # the raw csv data converted to the integer day of the year
32 # and continues through the end date in the raw csv data
33 # converted to the integer day of the year.
34 data_matrix[,1] <- c(start_doy:end_doy)
32 # Minimum 35 # Minimum
33 data_matrix[,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 data_matrix[,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(data_matrix, 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, start_doy, end_doy)
42 { 45 {
43 # From Forsythe 1995. 46 # From Forsythe 1995.
44 p=0.8333 47 p=0.8333
45 dl <- NULL 48 dl <- NULL
46 for (i in 1:num_days) { 49 for (i in start_doy:end_doy) {
47 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186))) 50 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186)))
48 phi <- asin(0.39795 * cos(theta)) 51 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))) 52 dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)))
50 } 53 }
51 # Return a vector of daylength for the number of 54 # Return a vector of daylength for the number of
68 dd <- 0 71 dd <- 0
69 } 72 }
70 else { 73 else {
71 # Extract daylength data for the number of 74 # Extract daylength data for the number of
72 # days specified in the input temperature data. 75 # days specified in the input temperature data.
73 dlprofile <- daylength(latitude, num_days) 76 dlprofile <- daylength(latitude, start_doy, end_doy)
74 # Initialize hourly temperature. 77 # Initialize hourly temperature.
75 T <- NULL 78 T <- NULL
76 # Initialize degree hour vector. 79 # Initialize degree hour vector.
77 dh <- NULL 80 dh <- NULL
78 # Calculate daylength in given date. 81 # Calculate daylength in given date.
197 return = mort.prob 200 return = mort.prob
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 # The input data currently must have 6 columns:
206 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
207 temperature_data <- read.csv(file=opt$input, header=T, strip.white=TRUE, sep=",")
203 start_date <- temperature_data[c(1:1), 3] 208 start_date <- temperature_data[c(1:1), 3]
209 start_doy <- temperature_data[c(1:1), 4]
204 end_date <- temperature_data[c(opt$num_days:opt$num_days), 3] 210 end_date <- temperature_data[c(opt$num_days:opt$num_days), 3]
211 end_doy <- temperature_data[c(opt$num_days:opt$num_days), 4]
205 raw_data_matrix <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days) 212 raw_data_matrix <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days)
206 temperature_file_path <- convert_csv_to_rdata(temperature_data, raw_data_matrix) 213 temperature_file_path <- convert_csv_to_rdata(start_doy, end_doy, temperature_data, raw_data_matrix)
207 latitude <- temperature_data[1, 1] 214 latitude <- temperature_data[1, 1]
208 215
216 cat("Start date: ", start_date, "\n")
217 cat("End date: ", end_date, "\n")
209 cat("Number of days: ", opt$num_days, "\n") 218 cat("Number of days: ", opt$num_days, "\n")
210 219
211 # Initialize matrix for results from all replications. 220 # 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) 221 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications)
213 newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol=opt$replications) 222 newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol=opt$replications)
221 vec.ini <- c(0, 3, 0, 0, 0) 230 vec.ini <- c(0, 3, 0, 0, 0)
222 # Overwintering, previttelogenic, DD=0, T=0, no-diapause. 231 # Overwintering, previttelogenic, DD=0, T=0, no-diapause.
223 vec.mat <- rep(vec.ini, n) 232 vec.mat <- rep(vec.ini, n)
224 # Complete matrix for the population. 233 # Complete matrix for the population.
225 vec.mat <- base::t(matrix(vec.mat, nrow=5)) 234 vec.mat <- base::t(matrix(vec.mat, nrow=5))
226 # Complete photoperiod profile in a year, requires daylength function. 235 # Complete photoperiod profile for the total
227 ph.p <- daylength(latitude, opt$num_days) 236 # number of days in the input data.
237 ph.p <- daylength(latitude, start_doy, end_doy)
228 238
229 # Time series of population size. 239 # Time series of population size.
230 tot.pop <- NULL 240 tot.pop <- NULL
231 gen0.pop <- rep(0, opt$num_days) 241 gen0.pop <- rep(0, opt$num_days)
232 gen1.pop <- rep(0, opt$num_days) 242 gen1.pop <- rep(0, opt$num_days)
235 g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days) 245 g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days)
236 N.newborn <- N.death <- N.adult <- rep(0, opt$num_days) 246 N.newborn <- N.death <- N.adult <- rep(0, opt$num_days)
237 dd.day <- rep(0, opt$num_days) 247 dd.day <- rep(0, opt$num_days)
238 248
239 # All the days included in the input temperature dataset. 249 # All the days included in the input temperature dataset.
240 for (day in 1:opt$num_days) { 250 for (day in start_doy:end_doy) {
241 # Photoperiod in the day. 251 # Photoperiod in the day.
242 photoperiod <- ph.p[day] 252 photoperiod <- ph.p[day]
243 temp.profile <- hourtemp(latitude, day, temperature_file_path, opt$num_days) 253 temp.profile <- hourtemp(latitude, day, temperature_file_path, opt$num_days)
244 mean.temp <- temp.profile[1] 254 mean.temp <- temp.profile[1]
245 dd.temp <- temp.profile[2] 255 dd.temp <- temp.profile[2]
532 g0a.rep[,N.rep] <- g0.adult 542 g0a.rep[,N.rep] <- g0.adult
533 g1a.rep[,N.rep] <- g1.adult 543 g1a.rep[,N.rep] <- g1.adult
534 g2a.rep[,N.rep] <- g2.adult 544 g2a.rep[,N.rep] <- g2.adult
535 } 545 }
536 546
537 # Data analysis and visualization 547 # Data analysis and visualization can currently
538 # default: plot 1 year of result 548 # plot only within a single calendar year. TODO:
539 # but can be expanded to accommodate multiple years 549 # enhance this to accomodate multiple calendar years.
540 n.yr <- 1 550 #n.yr <- 1
541 day.all <- c(1:opt$num_days * n.yr) 551 #day.all <- c(1:opt$num_days * n.yr)
552 day.all <- c(start_doy:end_doy)
542 553
543 # mean value for adults 554 # mean value for adults
544 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean) 555 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean)
545 # mean value for nymphs 556 # mean value for nymphs
546 sn <- apply((S1.rep + S2.rep), 1,mean) 557 sn <- apply((S1.rep + S2.rep), 1,mean)