insect_phenology_model: insect_phenology

comparison insect_phenology_model.R @ 56:829518206949 draft

Uploaded

author	greg
date	Mon, 05 Nov 2018 14:25:19 -0500
parents	5bb1d76c29ca
children	20c1451d8869

comparison

equal deleted inserted replaced

-:8c2b9eb9f9af
+:829518206949
 add_daylight_length = function(temperature_data_frame) {
 # Return temperature_data_frame with an added column
 # of daylight length (photoperiod profile).
 num_rows = dim(temperature_data_frame)[1];
-# From Forsythe 1995.
-p = 0.8333;
 latitude = temperature_data_frame$LATITUDE[1];
 daylight_length_vector = NULL;
 for (i in 1:num_rows) {
 # Get the day of the year from the current row
 # of the temperature data for computation.
 doy = temperature_data_frame$DOY[i];
-theta = 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (doy - 186)));
+daylight_length_vector[i] = get_daylen(doy, latitude);
-phi = asin(0.39795 * cos(theta));
-# Compute the length of daylight for the day of the year.
-darkness_length = 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)));
-daylight_length_vector[i] = 24 - darkness_length;
 }
 # Append daylight_length_vector as a new column to temperature_data_frame.
 temperature_data_frame = append_vector(temperature_data_frame, daylight_length_vector, "DAYLEN");
 return(temperature_data_frame);
 }
 for (i in first_norm_row:last_norm_row) {
 j = j + 1;
 tmp_data_frame[j,] = get_next_normals_row(norm_data_frame, year, i);
 }
 return (tmp_data_frame);
+}
+get_daylen = function(doy, latitude, p=0.8333) {
+# The default value for p is from Forsythe 1995.
+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.
+darkness_length = 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi)));
+daylight_length = 24 - darkness_length;
+return (daylight_length);
 }
 get_new_norm_data_frame = function(is_leap_year, input_norm=NULL, nrow=0) {
 # The input_norm data has the following 10 columns:
 # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX
 mmdd = norm_data_frame[index,"MMDD"][1];
 date_str = paste(year, mmdd, sep="-");
 doy = norm_data_frame[index,"DOY"][1];
 tmin = norm_data_frame[index,"TMIN"][1];
 tmax = norm_data_frame[index,"TMAX"][1];
-return(list(latitude, longitude, date_str, doy, tmin, tmax));
+daylen = get_daylen(doy, latitude);
-}
+# Calculate the average temperature based on the
+# tmin and tmax values from the 30 year normals data.
-get_temperature_at_hour = function(latitude, temperature_data_frame, row) {
+temp.profile = get_temperature_at_hour(latitude, row, tmin=tmin, tmax=tmax, daylen=daylen);
+mean_temp = temp.profile[1];
+avg_temp = temp.profile[2];
+return(list(latitude, longitude, date_str, doy, tmin, tmax, daylen, avg_temp));
+}
+get_temperature_at_hour = function(latitude, row, temperature_data_frame=NULL, tmin=NULL, tmax=NULL, daylen=NULL) {
 # Base development threshold for Brown Marmorated Stink Bug
 # insect phenology model.
 threshold = 14.17;
-# Minimum temperature for current row.
+if (is.null(temperature_data_frame)) {
-curr_min_temp = temperature_data_frame$TMIN[row];
+# The values of tmin and tmax cannot be NULL.
-# Maximum temperature for current row.
+curr_min_temp = tmin;
-curr_max_temp = temperature_data_frame$TMAX[row];
+curr_max_temp = tmax;
+} else {
+# Minimum temperature for current row.
+curr_min_temp = temperature_data_frame$TMIN[row];
+# Maximum temperature for current row.
+curr_max_temp = temperature_data_frame$TMAX[row];
+}
 # Mean temperature for current row.
 curr_mean_temp = 0.5 * (curr_min_temp + curr_max_temp);
 # Initialize degree day accumulation
 averages = 0;
 if (curr_max_temp < threshold) {
 else {
 # Initialize hourly temperature.
 T = NULL;
 # Initialize degree hour vector.
 dh = NULL;
-# Daylight length for current row.
+if (is.null(temperature_data_frame)) {
-y = temperature_data_frame$DAYLEN[row];
+# The value of daylen cannot be NULL.
+y = daylen;
+} else {
+# Daylight length for current row.
+y = temperature_data_frame$DAYLEN[row];
+}
 # Darkness length.
 z = 24 - y;
 # Lag coefficient.
 a = 1.86;
 # Darkness coefficient.
 num_days = 365;
 }
 cat("Number of days in year: ", num_days, "\n");
 }
-# Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately.
-if (plot_generations_separately) {
-temperature_data_frame_P = data.frame(temperature_data_frame);
-temperature_data_frame_F1 = data.frame(temperature_data_frame);
-temperature_data_frame_F2 = data.frame(temperature_data_frame);
-}
 # Get the ticks date labels for plots.
 ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, prepend_end_doy_norm, append_start_doy_norm);
 ticks = c(unlist(ticks_and_labels[1]));
 date_labels = c(unlist(ticks_and_labels[2]));
 # All latitude values are the same, so get the value for plots from the first row.
 for (row in 1:total_days) {
 # Get the integer day of the year for the current row.
 doy = temperature_data_frame$DOY[row];
 # Photoperiod in the day.
 photoperiod = temperature_data_frame$DAYLEN[row];
-temp.profile = get_temperature_at_hour(latitude, temperature_data_frame, row);
+temp.profile = get_temperature_at_hour(latitude, row, temperature_data_frame=temperature_data_frame);
 mean.temp = temp.profile[1];
 averages.temp = temp.profile[2];
 averages.day[row] = averages.temp;
 # Trash bin for death.
 death.vector = NULL;
 F2.total_adult[row] = sum((vector.matrix[,1]==2 & vector.matrix[,2]==3) | (vector.matrix[,1]==2 & vector.matrix[,2]==4) | (vector.matrix[,1]==2 & vector.matrix[,2]==5));
 }
 }
 }   # End of days specified in the input_ytd temperature data.
+# Set the cumulative average temperature (this is never used).
 averages.cum = cumsum(averages.day);
 # Define the output values.
 if (process_eggs) {
 Eggs.replications[,current_replication] = Eggs;
 }
 }
 population.replications[,current_replication] = total.population;
 # End processing replications.
 }
+# Append the averages.day vector (i.e., degree-days)
+# to the various temperature_data_frames.
+temperature_data_frame = append_vector(temperature_data_frame, averages.day, "DEGREE.DAYS");
 if (process_eggs) {
 # Mean value for eggs.
 eggs = apply(Eggs.replications, 1, mean);
 temperature_data_frame = append_vector(temperature_data_frame, eggs, "EGG");
 }
 }
 }
 if (plot_generations_separately) {
+# Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately.
+temperature_data_frame_P = data.frame(temperature_data_frame);
+temperature_data_frame_F1 = data.frame(temperature_data_frame);
+temperature_data_frame_F2 = data.frame(temperature_data_frame);
 m_se = get_mean_and_std_error(P.replications, F1.replications, F2.replications);
 P = m_se[[1]];
 P.std_error = m_se[[2]];
 F1 = m_se[[3]];
 F1.std_error = m_se[[4]];

Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 56:829518206949 draft