Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 64:e2f2780f5c14 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Tue, 14 Nov 2017 11:35:11 -0500
parents 46f360025893
children 8f3dba5dddb5
comparison
equal deleted inserted replaced
63:3e94de63a382 64:e2f2780f5c14
9 make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"), 9 make_option(c("-e", "--location"), action="store", dest="location", help="Selected location"),
10 make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), 10 make_option(c("-f", "--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"),
11 make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), 11 make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"),
12 make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"), 12 make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"),
13 make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"), 13 make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"),
14 make_option(c("-m", "--num_columns"), action="store", dest="num_columns", type="integer", help="Total number of columns in the temperature dataset"),
14 make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"), 15 make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"),
15 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"), 16 make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"),
16 make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), 17 make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"),
17 make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), 18 make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"),
18 make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"), 19 make_option(c("-s", "--replications"), action="store", dest="replications", type="integer", help="Number of replications"),
23 24
24 parser <- OptionParser(usage="%prog [options] file", option_list=option_list) 25 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
25 args <- parse_args(parser, positional_arguments=TRUE) 26 args <- parse_args(parser, positional_arguments=TRUE)
26 opt <- args$options 27 opt <- args$options
27 28
28 get_daylight_length = function(latitude, temperature_data, num_days) 29 add_daylight_length = function(latitude, temperature_data_matrix, num_days) {
29 {
30 # Return a vector of daylight length (photoperido profile) for 30 # Return a vector of daylight length (photoperido profile) for
31 # the number of days specified in the input temperature data 31 # the number of days specified in the input temperature data
32 # (from Forsythe 1995). 32 # (from Forsythe 1995).
33 p = 0.8333 33 p = 0.8333
34 daylight_length_vector <- NULL 34 daylight_length_vector <- NULL
35 for (i in 1:num_days) { 35 for (i in 1:num_days) {
36 # Get the day of the year from the current row 36 # Get the day of the year from the current row
37 # of the temperature data for computation. 37 # of the temperature data for computation.
38 doy <- temperature_data[i, 4] 38 doy <- temperature_data_matrix[i, 4]
39 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186))) 39 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)))
40 phi <- asin(0.39795 * cos(theta)) 40 phi <- asin(0.39795 * cos(theta))
41 # Compute the length of daylight for the day of the year. 41 # Compute the length of daylight for the day of the year.
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)))) 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))))
43 } 43 }
44 daylight_length_vector 44 # Append daylight_length_vector as a new column to temperature_data_matrix.
45 } 45 temperature_data_matrix <- cbind(temperature_data_matrix, daylight_length_vector)
46 46 }
47 get_temperature_at_hour = function(latitude, temperature_data, daylight_length_vector, row, num_days) 47
48 { 48 get_temperature_at_hour = function(latitude, temperature_data_matrix, row, num_days) {
49 # Base development threshold for Brown Marmolated Stink Bug 49 # Base development threshold for Brown Marmolated Stink Bug
50 # insect phenology model. 50 # insect phenology model.
51 # TODO: Pass insect on the command line to accomodate more 51 # TODO: Pass insect on the command line to accomodate more
52 # the just the Brown Marmolated Stink Bub. 52 # the just the Brown Marmolated Stink Bub.
53 threshold <- 14.17 53 threshold <- 14.17
54 54
55 # Input temperature currently has the following columns. 55 # Input temperature currently has the following columns.
56 # # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX 56 # # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
57 # Minimum temperature for current row. 57 # Minimum temperature for current row.
58 dnp <- temperature_data[row, 5] 58 dnp <- temperature_data_matrix[row, 5]
59 # Maximum temperature for current row. 59 # Maximum temperature for current row.
60 dxp <- temperature_data[row, 6] 60 dxp <- temperature_data_matrix[row, 6]
61 # Mean temperature for current row. 61 # Mean temperature for current row.
62 dmean <- 0.5 * (dnp + dxp) 62 dmean <- 0.5 * (dnp + dxp)
63 # Initialize degree day accumulation 63 # Initialize degree day accumulation
64 dd <- 0 64 dd <- 0
65 if (dxp < threshold) { 65 if (dxp < threshold) {
69 # Initialize hourly temperature. 69 # Initialize hourly temperature.
70 T <- NULL 70 T <- NULL
71 # Initialize degree hour vector. 71 # Initialize degree hour vector.
72 dh <- NULL 72 dh <- NULL
73 # Daylight length for current row. 73 # Daylight length for current row.
74 y <- daylight_length_vector[row] 74 y <- temperature_data_matrix[row, 7]
75 # Darkness length. 75 # Darkness length.
76 z <- 24 - y 76 z <- 24 - y
77 # Lag coefficient. 77 # Lag coefficient.
78 a <- 1.86 78 a <- 1.86
79 # Darkness coefficient. 79 # Darkness coefficient.
120 } 120 }
121 return=c(dmean, dd) 121 return=c(dmean, dd)
122 return 122 return
123 } 123 }
124 124
125 dev.egg = function(temperature) 125 dev.egg = function(temperature) {
126 {
127 dev.rate= -0.9843 * temperature + 33.438 126 dev.rate= -0.9843 * temperature + 33.438
128 return = dev.rate 127 return = dev.rate
129 return 128 return
130 } 129 }
131 130
132 dev.young = function(temperature) 131 dev.young = function(temperature) {
133 {
134 n12 <- -0.3728 * temperature + 14.68 132 n12 <- -0.3728 * temperature + 14.68
135 n23 <- -0.6119 * temperature + 25.249 133 n23 <- -0.6119 * temperature + 25.249
136 dev.rate = mean(n12 + n23) 134 dev.rate = mean(n12 + n23)
137 return = dev.rate 135 return = dev.rate
138 return 136 return
139 } 137 }
140 138
141 dev.old = function(temperature) 139 dev.old = function(temperature) {
142 {
143 n34 <- -0.6119 * temperature + 17.602 140 n34 <- -0.6119 * temperature + 17.602
144 n45 <- -0.4408 * temperature + 19.036 141 n45 <- -0.4408 * temperature + 19.036
145 dev.rate = mean(n34 + n45) 142 dev.rate = mean(n34 + n45)
146 return = dev.rate 143 return = dev.rate
147 return 144 return
148 } 145 }
149 146
150 dev.emerg = function(temperature) 147 dev.emerg = function(temperature) {
151 {
152 emerg.rate <- -0.5332 * temperature + 24.147 148 emerg.rate <- -0.5332 * temperature + 24.147
153 return = emerg.rate 149 return = emerg.rate
154 return 150 return
155 } 151 }
156 152
157 mortality.egg = function(temperature) 153 mortality.egg = function(temperature) {
158 {
159 if (temperature < 12.7) { 154 if (temperature < 12.7) {
160 mort.prob = 0.8 155 mort.prob = 0.8
161 } 156 }
162 else { 157 else {
163 mort.prob = 0.8 - temperature / 40.0 158 mort.prob = 0.8 - temperature / 40.0
167 } 162 }
168 return = mort.prob 163 return = mort.prob
169 return 164 return
170 } 165 }
171 166
172 mortality.nymph = function(temperature) 167 mortality.nymph = function(temperature) {
173 {
174 if (temperature < 12.7) { 168 if (temperature < 12.7) {
175 mort.prob = 0.03 169 mort.prob = 0.03
176 } 170 }
177 else { 171 else {
178 mort.prob = temperature * 0.0008 + 0.03 172 mort.prob = temperature * 0.0008 + 0.03
179 } 173 }
180 return = mort.prob 174 return = mort.prob
181 return 175 return
182 } 176 }
183 177
184 mortality.adult = function(temperature) 178 mortality.adult = function(temperature) {
185 {
186 if (temperature < 12.7) { 179 if (temperature < 12.7) {
187 mort.prob = 0.002 180 mort.prob = 0.002
188 } 181 }
189 else { 182 else {
190 mort.prob = temperature * 0.0005 + 0.02 183 mort.prob = temperature * 0.0005 + 0.02
194 } 187 }
195 188
196 # Read in the input temperature datafile into a Data Frame object. 189 # Read in the input temperature datafile into a Data Frame object.
197 # The input data currently must have 6 columns: 190 # The input data currently must have 6 columns:
198 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX 191 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
199 temperature_data <- read.csv(file=opt$input, header=T, strip.white=TRUE, sep=",") 192 temperature_data_matrix <- read.csv(file=opt$input, header=T, strip.white=TRUE, sep=",")
200 start_date <- temperature_data[1, 3] 193 start_date <- temperature_data_matrix[1, 3]
201 end_date <- temperature_data[opt$num_days, 3] 194 end_date <- temperature_data_matrix[opt$num_days, 3]
202 latitude <- temperature_data[1, 1] 195 latitude <- temperature_data_matrix[1, 1]
203 daylight_length_vector <- get_daylight_length(latitude, temperature_data, opt$num_days) 196 add_daylight_length(latitude, temperature_data_matrix, opt$num_days)
204 197
205 cat("Number of days: ", opt$num_days, "\n") 198 cat("Number of days: ", opt$num_days, "\n")
206 199
207 # Initialize matrix for results from all replications. 200 # Initialize matrix for results from all 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) 201 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications)
230 dd.day <- rep(0, opt$num_days) 223 dd.day <- rep(0, opt$num_days)
231 224
232 # All the days included in the input temperature dataset. 225 # All the days included in the input temperature dataset.
233 for (row in 1:opt$num_days) { 226 for (row in 1:opt$num_days) {
234 # Get the integer day of the year for the current row. 227 # Get the integer day of the year for the current row.
235 doy <- temperature_data[row, 4] 228 doy <- temperature_data_matrix[row, 4]
236 # Photoperiod in the day. 229 # Photoperiod in the day.
237 photoperiod <- daylight_length_vector[row] 230 photoperiod <- temperature_data_matrix[row, 7]
238 temp.profile <- get_temperature_at_hour(latitude, temperature_data, daylight_length_vector, row, opt$num_days) 231 temp.profile <- get_temperature_at_hour(latitude, temperature_data_matrix, row, opt$num_days)
239 mean.temp <- temp.profile[1] 232 mean.temp <- temp.profile[1]
240 dd.temp <- temp.profile[2] 233 dd.temp <- temp.profile[2]
241 dd.day[row] <- dd.temp 234 dd.day[row] <- dd.temp
242 # Trash bin for death. 235 # Trash bin for death.
243 death.vec <- NULL 236 death.vec <- NULL