Mercurial > repos > greg > insect_phenology_model
diff insect_phenology_model.R @ 142:b06b3881ecf0 draft
Uploaded
author | greg |
---|---|
date | Thu, 20 Dec 2018 09:08:50 -0500 |
parents | 14afb6e85581 |
children | 7628ba67c7ff |
line wrap: on
line diff
--- a/insect_phenology_model.R Wed Nov 21 11:33:40 2018 -0500 +++ b/insect_phenology_model.R Thu Dec 20 09:08:50 2018 -0500 @@ -250,50 +250,62 @@ return(mortality.probability) } -mortality.egg = function(temperature, adj=0) { - # If no input from adjustment, default - # value is 0 (data from Nielsen, 2008). - T.mortality = c(15, 17, 20, 25, 27, 30, 33, 35); - egg.mortality = c(50, 2, 1, 0, 0, 0, 5, 100); - # Calculates slopes and intercepts for lines. - slopes = NULL; - intercepts = NULL; - for (i in 1:length(T.mortality)) { - slopes[i] = (egg.mortality[i+1] - egg.mortality[i]) / (T.mortality[i+1] - T.mortality[i]); - intercepts[i] = -slopes[i] * T.mortality[i] + egg.mortality[i]; +#mortality.egg = function(temperature, adj=0) { +# # If no input from adjustment, default +# # value is 0 (data from Nielsen, 2008). +# T.mortality = c(15, 17, 20, 25, 27, 30, 33, 35); +# egg.mortality = c(50, 2, 1, 0, 0, 0, 5, 100); +# # Calculates slopes and intercepts for lines. +# slopes = NULL; +# intercepts = NULL; +# for (i in 1:length(T.mortality)) { +# slopes[i] = (egg.mortality[i+1] - egg.mortality[i]) / (T.mortality[i+1] - T.mortality[i]); +# intercepts[i] = -slopes[i] * T.mortality[i] + egg.mortality[i]; +# } +# # Calculates mortality based on temperature. +# mortality.probability = NULL; +# for (j in 1:length(temperature)) { +# mortality.probability[j] = if(temperature[j] <= T.mortality[2]) { +# temperature[j] * slopes[1] + intercepts[1]; +# } else if (temperature[j] > T.mortality[2] && temperature[j] <= T.mortality[3]) { +# temperature[j] * slopes[2] + intercepts[2]; +# } else if (temperature[j] > T.mortality[3] && temperature[j] <= T.mortality[4]) { +# temperature[j] * slopes[3] + intercepts[3]; +# } else if (temperature[j] > T.mortality[4] && temperature[j] <= T.mortality[5]) { +# temperature[j] * slopes[4] + intercepts[4]; +# } else if (temperature[j] > T.mortality[5] && temperature[j] <= T.mortality[6]) { +# temperature[j] * slopes[5] + intercepts[5]; +# } else if (temperature[j] > T.mortality[6] && temperature[j] <= T.mortality[7]) { +# temperature[j] * slopes[6] + intercepts[6]; +# } else if (temperature[j] > T.mortality[7]) { +# temperature[j] * slopes[7] + intercepts[7]; +# } +# # If mortality > 100, make it equal to 100. +# mortality.probability[mortality.probability>100] = 100; +# # If mortality <0, make equal to 0. +# mortality.probability[mortality.probability<0] = 0; +# } +# # Make mortality adjustments based on adj parameter. +# mortality.probability = (100 - mortality.probability) * adj + mortality.probability; +# # if mortality > 100, make it equal to 100. +# mortality.probability[mortality.probability>100] = 100; +# # If mortality <0, make equal to 0. +# mortality.probability[mortality.probability<0] = 0; +# # Change percent to proportion. +# mortality.probability = mortality.probability / 100; +# return(mortality.probability) +#} + +mortality.egg = function(temperature) { + if (temperature < 12.7) { + mortality.probability = 0.8; + } else { + mortality.probability = 0.8 - temperature / 40.0; + if (mortality.probability < 0) { + mortality.probability = 0.01; + } } - # Calculates mortality based on temperature. - mortality.probability = NULL; - for (j in 1:length(temperature)) { - mortality.probability[j] = if(temperature[j] <= T.mortality[2]) { - temperature[j] * slopes[1] + intercepts[1]; - } else if (temperature[j] > T.mortality[2] && temperature[j] <= T.mortality[3]) { - temperature[j] * slopes[2] + intercepts[2]; - } else if (temperature[j] > T.mortality[3] && temperature[j] <= T.mortality[4]) { - temperature[j] * slopes[3] + intercepts[3]; - } else if (temperature[j] > T.mortality[4] && temperature[j] <= T.mortality[5]) { - temperature[j] * slopes[4] + intercepts[4]; - } else if (temperature[j] > T.mortality[5] && temperature[j] <= T.mortality[6]) { - temperature[j] * slopes[5] + intercepts[5]; - } else if (temperature[j] > T.mortality[6] && temperature[j] <= T.mortality[7]) { - temperature[j] * slopes[6] + intercepts[6]; - } else if (temperature[j] > T.mortality[7]) { - temperature[j] * slopes[7] + intercepts[7]; - } - # If mortality > 100, make it equal to 100. - mortality.probability[mortality.probability>100] = 100; - # If mortality <0, make equal to 0. - mortality.probability[mortality.probability<0] = 0; - } - # Make mortality adjustments based on adj parameter. - mortality.probability = (100 - mortality.probability) * adj + mortality.probability; - # if mortality > 100, make it equal to 100. - mortality.probability[mortality.probability>100] = 100; - # If mortality <0, make equal to 0. - mortality.probability[mortality.probability<0] = 0; - # Change percent to proportion. - mortality.probability = mortality.probability / 100; - return(mortality.probability) + return (mortality.probability); } mortality.nymph = function(temperature) { @@ -311,6 +323,8 @@ prepend_end_doy_norm = 0; # The start DOY for norm data appended to ytd data. append_start_doy_norm = 0; + cat("start_date: ", start_date, "\n"); + cat("end_date: ", end_date, "\n"); if (is.null(start_date) && is.null(end_date)) { # We're not dealing with a date interval. date_interval = FALSE; @@ -326,6 +340,7 @@ start_date_doy = as.integer(strftime(start_date, format="%j")); end_date_doy = as.integer(strftime(end_date, format="%j")); } + cat("date_interval: ", date_interval, "\n"); if (is.null(input_ytd)) { # We're processing only the 30 year normals data. processing_year_to_date_data = FALSE; @@ -363,16 +378,13 @@ } end_date_ytd_row = which(temperature_data_frame$DATE==end_date); if (length(end_date_ytd_row) > 0) { + cat("I'm here...\n"); end_date_ytd_row = end_date_ytd_row[1]; # The end date is contained within the input_ytd data. end_doy_ytd = as.integer(temperature_data_frame$DOY[end_date_ytd_row]); - if (end_doy_ytd > end_date_ytd_row + 1) { - # The input year-to-date dataset is missing 1 or more - # days of data. - days_missing = end_doy_ytd - end_date_ytd_row; - msg = cat("The year-to-date dataset is missing ", days_missing, " days of data.\n"); - stop_err(msg); - } + cat("end_doy_ytd: ", end_doy_ytd, "\n"); + cat("end_date_ytd_row: ", end_date_ytd_row, "\n"); + cat("start_date_ytd_row: ", start_date_ytd_row, "\n"); } else { end_date_ytd_row = 0; } @@ -396,13 +408,10 @@ # Save the first DOY to later check if start_date is Jan 1. start_doy_ytd = as.integer(temperature_data_frame$DOY[1]); end_doy_ytd = as.integer(temperature_data_frame$DOY[num_ytd_rows]); - if (end_doy_ytd > end_date_ytd_row + 1) { - # The input year-to-date dataset is missing 1 or more - # days of data. - days_missing = end_doy_ytd - end_date_ytd_row; - msg = cat("The year-to-date dataset is missing ", days_missing, " days of data.\n"); - stop_err(msg); - } + cat("I'm here 2...\n"); + cat("end_doy_ytd: ", end_doy_ytd, "\n"); + cat("end_date_ytd_row: ", end_date_ytd_row, "\n"); + cat("start_date_ytd_row: ", start_date_ytd_row, "\n"); } } else { # We're processing only the 30 year normals data, so create an empty @@ -557,6 +566,8 @@ } } } + # Ensure all DOY values are consectuive integers. + validate_doys(temperature_data_frame); # Add a column containing the daylight length for each day. temperature_data_frame = add_daylight_length(temperature_data_frame); return(list(temperature_data_frame, start_date, end_date, prepend_end_doy_norm, append_start_doy_norm, is_leap_year, location)); @@ -856,7 +867,8 @@ } if (vector.individual[2] == 0) { # Egg. - death.probability = opt$egg_mortality * mortality.egg(mean.temp, adj=opt$egg_mortality); + # death.probability = opt$egg_mortality * mortality.egg(mean.temp, adj=opt$egg_mortality); + death.probability = opt$egg_mortality * mortality.egg(mean.temp); } else if (vector.individual[2] == 1 | vector.individual[2] == 2) { # Nymph.