Mercurial > repos > greg > insect_phenology_model
changeset 60:393085589438 draft
Uploaded
| author | greg |
|---|---|
| date | Thu, 20 Dec 2018 09:15:52 -0500 |
| parents | 892cf703be62 |
| children | 734f8e4dfbe5 |
| files | insect_phenology_model.R insect_phenology_model.xml test-data/output_combined1.csv test-data/output_combined2.csv test-data/output_combined3.csv test-data/output_combined4.csv test-data/output_f1_3.csv test-data/output_f1_4.csv test-data/output_f2_3.csv test-data/output_f2_4.csv test-data/output_p_3.csv test-data/output_p_4.csv utils.R |
| diffstat | 13 files changed, 92 insertions(+), 73 deletions(-) [+] |
line wrap: on
line diff
--- a/insect_phenology_model.R Wed Nov 21 11:42:37 2018 -0500 +++ b/insect_phenology_model.R Thu Dec 20 09:15:52 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) { @@ -366,13 +378,6 @@ 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); - } } else { end_date_ytd_row = 0; } @@ -396,13 +401,6 @@ # 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); - } } } else { # We're processing only the 30 year normals data, so create an empty @@ -557,6 +555,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 +856,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.
--- a/insect_phenology_model.xml Wed Nov 21 11:42:37 2018 -0500 +++ b/insect_phenology_model.xml Thu Dec 20 09:15:52 2018 -0500 @@ -66,13 +66,13 @@ <validator type="expression" message="30 year normals temperature data must have 10 columns and 366 rows">value is not None and value.metadata.columns==10 and value.metadata.data_lines==366</validator> </param> <conditional name="merge_ytd_temperature_data_cond"> - <param name="merge_ytd_temperature_data" type="select" label="Merge year-to-date temperature data with 30 year normals temperature data?"> + <param name="merge_ytd_temperature_data" type="select" label="Merge daily actuals temperature data with 30 year normals temperature data?"> <option value="yes" selected="true">Yes</option> <option value="no">No</option> </param> <when value="yes"> - <param name="input_ytd" type="data" format="csv" label="Year-to-date temperature data"> - <validator type="expression" message="Year-to-date temperature data must have 6 columns">value is not None and value.metadata.columns==6</validator> + <param name="input_ytd" type="data" format="csv" label="Daily actuals temperature data"> + <validator type="expression" message="Daily actuals temperature data must have 6 columns">value is not None and value.metadata.columns==6</validator> </param> <param name="location" type="text" value="" optional="true" label="Location" help="Enter the location or leave blank to use the station name from 30 year normals data."/> </when> @@ -81,7 +81,7 @@ <param name="insect" type="select" label="Select insect"> <option value="BMSB" selected="True">Brown Marmorated Stink Bug</option> </param> - <param name="replications" type="integer" value="10" min="2" label="Number of replications"/> + <param name="replications" type="integer" value="5" min="2" label="Number of replications"/> <param name="insects_per_replication" type="integer" value="1000" min="100" label="Number of insects with which to start each replication"/> <param name="photoperiod" type="float" value="13.5" min="0" label="Critical photoperiod for diapause induction/termination"/> <param name="egg_mortality" type="float" value="1" min="0" max="1" label="Adjustment rate for egg mortality" help="Floating point value between 0 and 1"/> @@ -242,7 +242,7 @@ * **30 year normals temperature data** - the dataset from your history containing the 30-year normals temperature data (available at http://pestwatch.psu.edu/ghcn). * **Merge year-to-date temperature data with 30 year normals temperature data** - select Yes to merge a year-to-date temperature dataset from your history into the selected 30 year normals temperature data. - * **Year-to-date temperature data** - the dataset from your history containing the year-to-date temperature data (available at http://pestwatch.psu.edu/minmax). + * **Daily actuals temperature data** - the dataset from your history containing the daily actuals temperature data (available at http://pestwatch.psu.edu/minmax). * **Location** - the location associated with the selected temperature data. * **Select insect** - currently only the Brown Marmorated Stink Bug can be analyzed.
--- a/test-data/output_combined1.csv Wed Nov 21 11:42:37 2018 -0500 +++ b/test-data/output_combined1.csv Thu Dec 20 09:15:52 2018 -0500 @@ -1,1 +1,1 @@ -"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","YOUNGNYMPH","YOUNGNYMPHSE","PRE.VITADULT","PRE.VITADULTSE" +"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","DEGREE.DAYS","YOUNGNYMPH","YOUNGNYMPHSE","PRE.VITADULT","PRE.VITADULTSE"
--- a/test-data/output_combined2.csv Wed Nov 21 11:42:37 2018 -0500 +++ b/test-data/output_combined2.csv Thu Dec 20 09:15:52 2018 -0500 @@ -1,1 +1,1 @@ -"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","EGG","EGGSE","TOTALNYMPH","TOTALNYMPHSE","TOTALADULT","TOTALADULTSE" +"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","DEGREE.DAYS","EGG","EGGSE","TOTALNYMPH","TOTALNYMPHSE","TOTALADULT","TOTALADULTSE"
--- a/test-data/output_combined3.csv Wed Nov 21 11:42:37 2018 -0500 +++ b/test-data/output_combined3.csv Thu Dec 20 09:15:52 2018 -0500 @@ -1,1 +1,1 @@ -"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","OLDNYMPH","OLDNYMPHSE" +"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","DEGREE.DAYS","OLDNYMPH","OLDNYMPHSE"
--- a/test-data/output_combined4.csv Wed Nov 21 11:42:37 2018 -0500 +++ b/test-data/output_combined4.csv Thu Dec 20 09:15:52 2018 -0500 @@ -1,1 +1,1 @@ -"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","EGG","EGGSE","TOTALNYMPH","TOTALNYMPHSE","TOTALADULT","TOTALADULTSE" +"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","DEGREE.DAYS","EGG","EGGSE","TOTALNYMPH","TOTALNYMPHSE","TOTALADULT","TOTALADULTSE"
--- a/test-data/output_f1_3.csv Wed Nov 21 11:42:37 2018 -0500 +++ b/test-data/output_f1_3.csv Thu Dec 20 09:15:52 2018 -0500 @@ -1,1 +1,1 @@ -"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","OLDNYMPH.F1","OLDNYMPH.F1.SE" +"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","DEGREE.DAYS","OLDNYMPH","OLDNYMPHSE","OLDNYMPH.F1","OLDNYMPH.F1.SE"
--- a/test-data/output_f1_4.csv Wed Nov 21 11:42:37 2018 -0500 +++ b/test-data/output_f1_4.csv Thu Dec 20 09:15:52 2018 -0500 @@ -1,1 +1,1 @@ -"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","EGG.F1","EGG.F1.SE","TOTALNYMPH.F1","TOTALNYMPH.F1.SE","TOTALADULT.F1","TOTALADULT.F1.SE" +"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","DEGREE.DAYS","EGG","EGGSE","TOTALNYMPH","TOTALNYMPHSE","TOTALADULT","TOTALADULTSE","EGG.F1","EGG.F1.SE","TOTALNYMPH.F1","TOTALNYMPH.F1.SE","TOTALADULT.F1","TOTALADULT.F1.SE","ALL.TOTAL.F1","ALL.TOTAL.F1.SE"
--- a/test-data/output_f2_3.csv Wed Nov 21 11:42:37 2018 -0500 +++ b/test-data/output_f2_3.csv Thu Dec 20 09:15:52 2018 -0500 @@ -1,1 +1,1 @@ -"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","OLDNYMPH.F2","OLDNYMPH.F2.SE" +"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","DEGREE.DAYS","OLDNYMPH","OLDNYMPHSE","OLDNYMPH.F2","OLDNYMPH.F2.SE"
--- a/test-data/output_f2_4.csv Wed Nov 21 11:42:37 2018 -0500 +++ b/test-data/output_f2_4.csv Thu Dec 20 09:15:52 2018 -0500 @@ -1,1 +1,1 @@ -"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","EGG.F2","EGG.F2.SE","TOTALNYMPH.F2","TOTALNYMPH.F2.SE","TOTALADULT.F2","TOTALADULT.F2.SE" +"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","DEGREE.DAYS","EGG","EGGSE","TOTALNYMPH","TOTALNYMPHSE","TOTALADULT","TOTALADULTSE","EGG.F2","EGG.F2.SE","TOTALNYMPH.F2","TOTALNYMPH.F2.SE","TOTALADULT.F2","TOTALADULT.F2.SE","ALL.TOTAL.F2","ALL.TOTAL.F2.SE"
--- a/test-data/output_p_3.csv Wed Nov 21 11:42:37 2018 -0500 +++ b/test-data/output_p_3.csv Thu Dec 20 09:15:52 2018 -0500 @@ -1,1 +1,1 @@ -"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","OLDNYMPH.P","OLDNYMPH.P.SE" +"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","DEGREE.DAYS","OLDNYMPH","OLDNYMPHSE","OLDNYMPH.P","OLDNYMPH.P.SE"
--- a/test-data/output_p_4.csv Wed Nov 21 11:42:37 2018 -0500 +++ b/test-data/output_p_4.csv Thu Dec 20 09:15:52 2018 -0500 @@ -1,1 +1,1 @@ -"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","EGG.P","EGG.P.SE","TOTALNYMPH.P","TOTALNYMPH.P.SE","TOTALADULT.P","TOTALADULT.P.SE","ALL.TOTAL.P","ALL.TOTAL.P.SE" +"LATITUDE","LONGITUDE","DATE","DOY","TMIN","TMAX","DAYLEN","DEGREE.DAYS","EGG","EGGSE","TOTALNYMPH","TOTALNYMPHSE","TOTALADULT","TOTALADULTSE","EGG.P","EGG.P.SE","TOTALNYMPH.P","TOTALNYMPH.P.SE","TOTALADULT.P","TOTALADULT.P.SE","ALL.TOTAL.P","ALL.TOTAL.P.SE"
--- a/utils.R Wed Nov 21 11:42:37 2018 -0500 +++ b/utils.R Thu Dec 20 09:15:52 2018 -0500 @@ -333,3 +333,21 @@ } return(valid_date); } + +validate_doys = function(temperature_data_frame) { + # Ensure all DOY values are consecutive integers. + last_doy = 0; + num_rows = dim(temperature_data_frame)[1]; + for (i in 1:num_rows) { + doy = as.integer(temperature_data_frame$DOY[i]); + if (last_doy == 0) { + last_doy = doy; + } else { + if (doy != last_doy + 1) { + msg = paste("DOY ", doy, " is not consecutive (previous DOY is ", last_doy, ".", sep=""); + stop_err(msg); + } + last_doy = doy; + } + } +}