insect_phenology_model: insect_phenology

comparison insect_phenology_model.R @ 73:9bc7ab4da9e7 draft

Uploaded

author	greg
date	Tue, 14 Nov 2017 13:54:27 -0500
parents	3c5a149d4255
children	9843894544a8

comparison

equal deleted inserted replaced

-:3c5a149d4255
+:9bc7ab4da9e7
 parser <- OptionParser(usage="%prog [options] file", option_list=option_list)
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
-add_daylight_length = function(latitude, temperature_data_matrix, num_days) {
+parse_input_data = function(input_file, num_rows, num_columns) {
+# Read in the input temperature datafile into a data frame.
+if (num_rows == 6) {
+# The input data has the following 6 columns:
+# LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
+temperature_data_frame <- read.csv(file=input_file, header=T, strip.white=TRUE, sep=",")
+# Add a column containing the daylight length for each day.
+temperature_data_frame <- add_daylight_length(latitude, temperature_data_frame, opt$num_days)
+temperature_data_frame
+}
+}
+add_daylight_length = function(latitude, temperature_data_frame, num_days) {
 # Return a vector of daylight length (photoperido profile) for
 # the number of days specified in the input temperature data
 # (from Forsythe 1995).
 p = 0.8333
 daylight_length_vector <- NULL
 for (i in 1:num_days) {
 # Get the day of the year from the current row
 # of the temperature data for computation.
-doy <- temperature_data_matrix[i, 4]
+doy <- temperature_data_frame[i, 4]
 theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)))
 phi <- asin(0.39795 * cos(theta))
 # Compute the length of daylight for the day of the year.
 daylight_length_vector[i] <- 24 - (24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi))))
 }
-# Append daylight_length_vector as a new column to temperature_data_matrix.
+# Append daylight_length_vector as a new column to temperature_data_frame.
-temperature_data_matrix[, 7] <- daylight_length_vector
+temperature_data_frame[, 7] <- daylight_length_vector
-# Return the altered temperature_data_matrix.
+# Return the altered temperature_data_frame.
-temperature_data_matrix
+temperature_data_frame
 }
-get_temperature_at_hour = function(latitude, temperature_data_matrix, row, num_days) {
+get_temperature_at_hour = function(latitude, temperature_data_frame, row, num_days) {
 # Base development threshold for Brown Marmolated Stink Bug
 # insect phenology model.
 # TODO: Pass insect on the command line to accomodate more
 # the just the Brown Marmolated Stink Bub.
 threshold <- 14.17
 # Input temperature currently has the following columns.
 # # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
 # Minimum temperature for current row.
-dnp <- temperature_data_matrix[row, 5]
+dnp <- temperature_data_frame[row, 5]
 # Maximum temperature for current row.
-dxp <- temperature_data_matrix[row, 6]
+dxp <- temperature_data_frame[row, 6]
 # Mean temperature for current row.
 dmean <- 0.5 * (dnp + dxp)
 # Initialize degree day accumulation
 dd <- 0
 if (dxp < threshold) {
 # Initialize hourly temperature.
 T <- NULL
 # Initialize degree hour vector.
 dh <- NULL
 # Daylight length for current row.
-y <- temperature_data_matrix[row, 7]
+y <- temperature_data_frame[row, 7]
 # Darkness length.
 z <- 24 - y
 # Lag coefficient.
 a <- 1.86
 # Darkness coefficient.
 }
 return = mort.prob
 return
 }
-# Read in the input temperature datafile into a Data Frame object.
+temperature_data_frame <- parse_input_data(opt$input, opt$num_rows, opt$num_columns)
-# The input data currently must have 6 columns:
+latitude <- temperature_data_frame[1, 1]
-# LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX
-temperature_data_matrix <- read.csv(file=opt$input, header=T, strip.white=TRUE, sep=",")
-start_date <- temperature_data_matrix[1, 3]
-end_date <- temperature_data_matrix[opt$num_days, 3]
-latitude <- temperature_data_matrix[1, 1]
-temperature_data_matrix <- add_daylight_length(latitude, temperature_data_matrix, opt$num_days)
 cat("Number of days: ", opt$num_days, "\n")
-print(temperature_data_matrix)
 # Initialize matrix for results from all replications.
 S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications)
 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)
 dd.day <- rep(0, opt$num_days)
 # All the days included in the input temperature dataset.
 for (row in 1:opt$num_days) {
 # Get the integer day of the year for the current row.
-doy <- temperature_data_matrix[row, 4]
+doy <- temperature_data_frame[row, 4]
 # Photoperiod in the day.
-photoperiod <- temperature_data_matrix[row, 7]
+photoperiod <- temperature_data_frame[row, 7]
-temp.profile <- get_temperature_at_hour(latitude, temperature_data_matrix, row, opt$num_days)
+temp.profile <- get_temperature_at_hour(latitude, temperature_data_frame, row, opt$num_days)
 mean.temp <- temp.profile[1]
 dd.temp <- temp.profile[2]
 dd.day[row] <- dd.temp
 # Trash bin for death.
 death.vec <- NULL
 }
 # Data analysis and visualization can currently
 # plot only within a single calendar year.
 # TODO: enhance this to accomodate multiple calendar years.
+start_date <- temperature_data_frame[1, 3]
+end_date <- temperature_data_frame[opt$num_days, 3]
 n.yr <- 1
 day.all <- c(1:opt$num_days * n.yr)
 # mean value for adults
 sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean)

Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 73:9bc7ab4da9e7 draft