Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 46:d7e406de8094 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 10:42:31 -0500
parents 5260db6479db
children d6482112bf35
comparison
equal deleted inserted replaced
45:16a913a6f7d7 46:d7e406de8094
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(start_doy, end_doy, temperature_data, data_matrix) 28 get_daylight_length = function(latitude, temperature_data, num_days)
29 { 29 {
30 # Make sure the first column starts with the start date in 30 # Return a vector of daylight length (photoperido profile) for
31 # the raw csv data converted to the integer day of the year 31 # the number of days specified in the input temperature data
32 # and continues through the end date in the raw csv data 32 # (from Forsythe 1995).
33 # converted to the integer day of the year. 33 p = 0.8333
34 data_matrix[,1] <- c(start_doy:end_doy) 34 daylight_length_vector <- NULL
35 # Minimum 35 for (i in 1:num_days) {
36 data_matrix[,2] <- temperature_data[c(1:opt$num_days), 5] 36 # Get the day of the year from the current row
37 # Maximum 37 # of the temperature data for computation.
38 data_matrix[,3] <- temperature_data[c(1:opt$num_days), 6] 38 doy <- temperature_data[c(i:i), 4]
39 namedat <- "tempdata.Rdat" 39 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)))
40 save(data_matrix, file=namedat)
41 namedat
42 }
43
44 daylength=function(latitude, start_doy, end_doy)
45 {
46 # From Forsythe 1995.
47 p=0.8333
48 dl <- NULL
49 for (i in start_doy:end_doy) {
50 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186)))
51 phi <- asin(0.39795 * cos(theta)) 40 phi <- asin(0.39795 * cos(theta))
52 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.
53 } 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))))
54 # Return a vector of daylength for the number of 43 }
55 # days specified in the input temperature data. 44 daylight_length_vector
56 dl 45 }
57 } 46
58 47 get_temperature_at_hour = function(latitude, temperature_data, daylight_length_vector, row, num_days)
59 hourtemp=function(latitude, date, temperature_file_path, num_days) 48 {
60 {
61 load(temperature_file_path)
62 # Base development threshold for Brown Marmolated Stink Bug 49 # Base development threshold for Brown Marmolated Stink Bug
63 # 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.
64 threshold <- 14.17 53 threshold <- 14.17
65 dnp <- data_matrix[date, 2] # daily minimum 54
66 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.
67 dmean <- 0.5 * (dnp + dxp) 62 dmean <- 0.5 * (dnp + dxp)
68 dd <- 0 # initialize degree day accumulation 63 dd <- 0 # initialize degree day accumulation
69 64 if (dxp < threshold) {
70 if (dxp<threshold) {
71 dd <- 0 65 dd <- 0
72 } 66 }
73 else { 67 else {
74 # Extract daylength data for the number of
75 # days specified in the input temperature data.
76 dlprofile <- daylength(latitude, start_doy, end_doy)
77 # Initialize hourly temperature. 68 # Initialize hourly temperature.
78 T <- NULL 69 T <- NULL
79 # Initialize degree hour vector. 70 # Initialize degree hour vector.
80 dh <- NULL 71 dh <- NULL
81 # Calculate daylength in given date. 72 # Daylight length for current row.
82 y <- dlprofile[date] 73 y <- daylight_length_vector[row]
83 # Night length. 74 # Darkness length.
84 z <- 24 - y 75 z <- 24 - y
85 # Lag coefficient. 76 # Lag coefficient.
86 a <- 1.86 77 a <- 1.86
87 # Night coefficient. 78 # Darkness coefficient.
88 b <- 2.20 79 b <- 2.20
89 # Sunrise time. 80 # Sunrise time.
90 risetime <- 12 - y / 2 81 risetime <- 12 - y / 2
91 # Sunset time. 82 # Sunset time.
92 settime <- 12 + y / 2 83 settime <- 12 + y / 2
93 ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp 84 ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp
94 for (i in 1:24) { 85 for (i in 1:24) {
95 if (i > risetime && i<settime) { 86 if (i > risetime && i < settime) {
96 # Number of hours after Tmin until sunset. 87 # Number of hours after Tmin until sunset.
97 m <- i - 5 88 m <- i - 5
98 T[i]=(dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp 89 T[i] = (dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp
99 if (T[i]<8.4) { 90 if (T[i] < 8.4) {
100 dh[i] <- 0 91 dh[i] <- 0
101 } 92 }
102 else { 93 else {
103 dh[i] <- T[i] - 8.4 94 dh[i] <- T[i] - 8.4
104 } 95 }
105 } 96 }
106 else if (i > settime) { 97 else if (i > settime) {
107 n <- i - settime 98 n <- i - settime
108 T[i]=dnp + (ts - dnp) * exp( - b * n / z) 99 T[i] = dnp + (ts - dnp) * exp( - b * n / z)
109 if (T[i]<8.4) { 100 if (T[i] < 8.4) {
110 dh[i] <- 0 101 dh[i] <- 0
111 } 102 }
112 else { 103 else {
113 dh[i] <- T[i] - 8.4 104 dh[i] <- T[i] - 8.4
114 } 105 }
115 } 106 }
116 else { 107 else {
117 n <- i + 24 - settime 108 n <- i + 24 - settime
118 T[i]=dnp + (ts - dnp) * exp( - b * n / z) 109 T[i]=dnp + (ts - dnp) * exp( - b * n / z)
119 if (T[i]<8.4) { 110 if (T[i] < 8.4) {
120 dh[i] <- 0 111 dh[i] <- 0
121 } 112 }
122 else { 113 else {
123 dh[i] <- T[i] - 8.4 114 dh[i] <- T[i] - 8.4
124 } 115 }
204 # Read in the input temperature datafile into a Data Frame object. 195 # Read in the input temperature datafile into a Data Frame object.
205 # The input data currently must have 6 columns: 196 # The input data currently must have 6 columns:
206 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX 197 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
207 temperature_data <- read.csv(file=opt$input, header=T, strip.white=TRUE, sep=",") 198 temperature_data <- read.csv(file=opt$input, header=T, strip.white=TRUE, sep=",")
208 start_date <- temperature_data[c(1:1), 3] 199 start_date <- temperature_data[c(1:1), 3]
209 start_doy <- temperature_data[c(1:1), 4]
210 end_date <- temperature_data[c(opt$num_days:opt$num_days), 3] 200 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]
212 raw_data_matrix <- matrix(rep(0, opt$num_days * 6), nrow=opt$num_days)
213 temperature_file_path <- convert_csv_to_rdata(start_doy, end_doy, temperature_data, raw_data_matrix)
214 latitude <- temperature_data[1, 1] 201 latitude <- temperature_data[1, 1]
202 daylight_length_vector <- get_daylight_length(latitude, temperature_data, opt$num_days)
215 203
216 cat("Start date: ", start_date, "\n") 204 cat("Start date: ", start_date, "\n")
217 cat("End date: ", end_date, "\n") 205 cat("End date: ", end_date, "\n")
218 cat("Number of days: ", opt$num_days, "\n") 206 cat("Number of days: ", opt$num_days, "\n")
219 207
230 vec.ini <- c(0, 3, 0, 0, 0) 218 vec.ini <- c(0, 3, 0, 0, 0)
231 # Overwintering, previttelogenic, DD=0, T=0, no-diapause. 219 # Overwintering, previttelogenic, DD=0, T=0, no-diapause.
232 vec.mat <- rep(vec.ini, n) 220 vec.mat <- rep(vec.ini, n)
233 # Complete matrix for the population. 221 # Complete matrix for the population.
234 vec.mat <- base::t(matrix(vec.mat, nrow=5)) 222 vec.mat <- base::t(matrix(vec.mat, nrow=5))
235 # Complete photoperiod profile for the total
236 # number of days in the input data.
237 ph.p <- daylength(latitude, start_doy, end_doy)
238
239 # Time series of population size. 223 # Time series of population size.
240 tot.pop <- NULL 224 tot.pop <- NULL
241 gen0.pop <- rep(0, opt$num_days) 225 gen0.pop <- rep(0, opt$num_days)
242 gen1.pop <- rep(0, opt$num_days) 226 gen1.pop <- rep(0, opt$num_days)
243 gen2.pop <- rep(0, opt$num_days) 227 gen2.pop <- rep(0, opt$num_days)
245 g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days) 229 g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days)
246 N.newborn <- N.death <- N.adult <- rep(0, opt$num_days) 230 N.newborn <- N.death <- N.adult <- rep(0, opt$num_days)
247 dd.day <- rep(0, opt$num_days) 231 dd.day <- rep(0, opt$num_days)
248 232
249 # All the days included in the input temperature dataset. 233 # All the days included in the input temperature dataset.
250 for (day in start_doy:end_doy) { 234 for (row in 1:opt$num_days) {
235 # Get the integer day of the year for the current row.
236 doy <- temperature_data[c(row:row), 4]
251 # Photoperiod in the day. 237 # Photoperiod in the day.
252 photoperiod <- ph.p[day] 238 photoperiod <- daylight_length_vector[row]
253 temp.profile <- hourtemp(latitude, day, temperature_file_path, opt$num_days) 239 temp.profile <- get_temperature_at_hour(latitude, temperature_data, daylight_length_vector, row, opt$num_days)
254 mean.temp <- temp.profile[1] 240 mean.temp <- temp.profile[1]
255 dd.temp <- temp.profile[2] 241 dd.temp <- temp.profile[2]
256 dd.day[day] <- dd.temp 242 dd.day[day] <- dd.temp
257 # Trash bin for death. 243 # Trash bin for death.
258 death.vec <- NULL 244 death.vec <- NULL
543 g1a.rep[,N.rep] <- g1.adult 529 g1a.rep[,N.rep] <- g1.adult
544 g2a.rep[,N.rep] <- g2.adult 530 g2a.rep[,N.rep] <- g2.adult
545 } 531 }
546 532
547 # Data analysis and visualization can currently 533 # Data analysis and visualization can currently
548 # plot only within a single calendar year. TODO: 534 # plot only within a single calendar year.
549 # enhance this to accomodate multiple calendar years. 535 # TODO: enhance this to accomodate multiple calendar years.
550 #n.yr <- 1 536 n.yr <- 1
551 #day.all <- c(1:opt$num_days * n.yr) 537 day.all <- c(1:opt$num_days * n.yr)
552 day.all <- c(start_doy:end_doy)
553 538
554 # mean value for adults 539 # mean value for adults
555 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean) 540 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean)
556 # mean value for nymphs 541 # mean value for nymphs
557 sn <- apply((S1.rep + S2.rep), 1,mean) 542 sn <- apply((S1.rep + S2.rep), 1,mean)