insect_phenology_model: insect_phenology

comparison insect_phenology_model.R @ 39:169c8180205a draft

Uploaded

author	greg
date	Wed, 11 Apr 2018 11:26:53 -0400
parents	c0f76f4f84fc
children	d8e6304dc5e4

comparison

equal deleted inserted replaced

-:c0f76f4f84fc
+:169c8180205a
 # F2 mean.
 f2_m = apply(f2_replications, 1, mean);
 # F2 standard error.
 f2_se = apply(f2_replications, 1, sd) / sqrt(opt$replications);
 return(list(p_m, p_se, f1_m, f1_se, f2_m, f2_se))
+}
+get_next_normals_row = function(norm_data_frame, year, is_leap_year, index) {
+# Return the next 30 year normals row formatted
+# appropriately for the year-to-date data.
+latitude = norm_data_frame[index,"LATITUDE"][1];
+longitude = norm_data_frame[index,"LONGITUDE"][1];
+# Format the date.
+mmdd = norm_data_frame[index,"MMDD"][1];
+date_str = paste(year, mmdd, sep="-");
+doy = norm_data_frame[index,"DOY"][1];
+if (!is_leap_year) {
+# Since all normals data includes Feb 29, we have to
+# subtract 1 from DOY if we're not in a leap year since
+# we removed the Feb 29 row from the data frame above.
+doy = as.integer(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));
 }
 get_temperature_at_hour = function(latitude, temperature_data_frame, row, num_days) {
 # Base development threshold for Brown Marmorated Stink Bug
 # insect phenology model.
 }
 get_total_days = function(is_leap_year) {
 # Get the total number of days in the current year.
 if (is_leap_year) {
-return (366);
+return(366);
 } else {
-return (365);
+return(365);
 }
 }
-get_x_axis_ticks_and_labels = function(temperature_data_frame, num_rows, num_days_ytd) {
+get_x_axis_ticks_and_labels = function(temperature_data_frame, num_rows, start_doy_ytd, end_doy_ytd) {
 # Keep track of the years to see if spanning years.
 month_labels = list();
 ticks = list();
 current_month_label = NULL;
 last_tick = 0;
 for (i in 1:num_rows) {
-if (i==num_days_ytd) {
+if (start_doy_ytd > 1 & i==start_doy_ytd-1) {
-# Add a tick for the start of the 30 year normnals data.
+# Add a tick for the end of the 30 year normnals data
+# that was prepended to the year-to-date data.
 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
 ticks[tick_index] = i;
-month_labels[tick_index] = "Start 30 year normals";
+month_labels[tick_index] = "End prepended 30 year normals";
 last_tick = i;
-}
+} else  if (i==end_doy_ytd+1) {
-# Get the year and month from the date which
+# Add a tick for the start of the 30 year normnals data
-# has the format YYYY-MM-DD.
+# that was appended to the year-to-date data.
-date = format(temperature_data_frame$DATE[i]);
+tick_index = get_tick_index(i, last_tick, ticks, month_labels)
-# Get the month label.
-items = strsplit(date, "-")[[1]];
-month = items[2];
-month_label = month.abb[as.integer(month)];
-tick_index = get_tick_index(i, last_tick, ticks, month_labels)
-if (!identical(current_month_label, month_label)) {
-# Add an x-axis tick for the month.
 ticks[tick_index] = i;
-month_labels[tick_index] = month_label;
+month_labels[tick_index] = "Start appended 30 year normals";
-current_month_label = month_label;
 last_tick = i;
-}
+} else if (i==num_rows) {
-tick_index = get_tick_index(i, last_tick, ticks, month_labels)
-if (!is.null(tick_index)) {
-# Get the day.
-day = weekdays(as.Date(date));
-if (day=="Sunday") {
-# Add an x-axis tick if we're on a Sunday.
-ticks[tick_index] = i;
-# Add a blank month label so it is not displayed.
-month_labels[tick_index] = "";
-last_tick = i;
-}
-}
-if (i==num_rows) {
 # Add a tick for the last day of the year.
 tick_index = get_tick_index(i, last_tick, ticks, month_labels)
 ticks[tick_index] = i;
 month_labels[tick_index] = "";
 last_tick = i;
+} else {
+# Get the year and month from the date which
+# has the format YYYY-MM-DD.
+date = format(temperature_data_frame$DATE[i]);
+# Get the month label.
+items = strsplit(date, "-")[[1]];
+month = items[2];
+month_label = month.abb[as.integer(month)];
+if (!identical(current_month_label, month_label)) {
+# Add an x-axis tick for the month.
+tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+ticks[tick_index] = i;
+month_labels[tick_index] = month_label;
+current_month_label = month_label;
+last_tick = i;
+}
+tick_index = get_tick_index(i, last_tick, ticks, month_labels)
+if (!is.null(tick_index)) {
+# Get the day.
+day = weekdays(as.Date(date));
+if (day=="Sunday") {
+# Add an x-axis tick if we're on a Sunday.
+ticks[tick_index] = i;
+# Add a blank month label so it is not displayed.
+month_labels[tick_index] = "";
+last_tick = i;
+}
+}
 }
 }
 return(list(ticks, month_labels));
 }
 # Extract the year from the date_str.
 date = format(date_str);
 items = strsplit(date, "-")[[1]];
 year = as.integer(items[1]);
 if (((year %% 4 == 0) & (year %% 100 != 0)) | (year %% 400 == 0)) {
-return (TRUE);
+return(TRUE);
 } else {
-return (FALSE);
+return(FALSE);
 }
 }
 mortality.adult = function(temperature) {
 if (temperature < 12.7) {
 colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX");
 # Get the start date.
 start_date = temperature_data_frame$DATE[1];
 # See if we're in a leap year.
 is_leap_year = is_leap_year(start_date);
-# get the number of days in the year.
+# Get the number of days in the year.
 total_days = get_total_days(is_leap_year);
 # Extract the year from the start date.
 date_str = format(start_date);
 date_str_items = strsplit(date_str, "-")[[1]];
 year = date_str_items[1];
+# 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_days_ytd]);
 # Read the input_norm temperature datafile into a data frame.
 # The input_norm data has the following 10 columns:
 # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX
 norm_data_frame = read.csv(file=input_norm, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=",");
 # Set the norm_data_frame column names for access.
 # All normals data includes Feb 29 which is row 60 in
 # the data, so delete that row if we're not in a leap year.
 if (!is_leap_year) {
 norm_data_frame = norm_data_frame[-c(60),];
 }
-# Define the next row for the temperature_data_frame from the 30 year normals data.
+if (start_doy_ytd > 1) {
-first_normals_row = num_days_ytd + 1;
+# The year-to-date data starts after Jan 1, so create a
-for (i in first_normals_row:total_days) {
+# temporary data frame to contain the 30 year normals data
-latitude = norm_data_frame[i,"LATITUDE"][1];
+# from Jan 1 to the date immediately prior to start_date.
-longitude = norm_data_frame[i,"LONGITUDE"][1];
+tmp_data_frame = temperature_data_frame[FALSE,];
-# Format the date.
+for (i in 1:start_doy_ytd-1) {
-mmdd = norm_data_frame[i,"MMDD"][1];
+tmp_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
-date_str = paste(year, mmdd, sep="-");
+}
-doy = norm_data_frame[i,"DOY"][1];
+# Next merge the temporary data frame with the year-to-date data frame.
-if (!is_leap_year) {
+temperature_data_frame = rbind(tmp_data_frame, temperature_data_frame);
-# Since all normals data includes Feb 29, we have to
+}
-# subtract 1 from DOY if we're not in a leap year since
+# Define the next row for the year-to-date data from the 30 year normals data.
-# we removed the Feb 29 row from the data frame above.
+first_normals_append_row = end_doy_ytd + 1;
-doy = as.integer(doy) - 1;
+# Append the 30 year normals data to the year-to-date data.
-}
+for (i in first_normals_append_row:total_days) {
-tmin = norm_data_frame[i,"TMIN"][1];
+temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
-tmax = norm_data_frame[i,"TMAX"][1];
-# Append the row to temperature_data_frame.
-new_row = list(latitude, longitude, date_str, doy, tmin, tmax);
-temperature_data_frame[i,] = new_row;
 }
 # Add a column containing the daylight length for each day.
 temperature_data_frame = add_daylight_length(temperature_data_frame, total_days);
-return(temperature_data_frame);
+return(list(temperature_data_frame, start_date, start_doy_ytd, end_doy_ytd, is_leap_year, total_days));
 }
 render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval,
 replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL,
 life_stages_adult=NULL, life_stages_nymph=NULL) {
 if (chart_type=="pop_size_by_life_stage") {
 if (life_stage=="Total") {
 title = paste(insect, ": Reps", replications, ":", life_stage, "Pop :", location, ": Lat", latitude, ":", start_date, "-", end_date, sep=" ");
 legend_text = c("Egg", "Nymph", "Adult");
 columns = c(4, 2, 1);
 plot_generations_separately = FALSE;
 }
 # Display the total number of days in the Galaxy history item blurb.
 cat("Year-to-date number of days: ", opt$num_days_ytd, "\n");
-# Read the temperature data into a data frame.
+# Parse the inputs.
-temperature_data_frame = parse_input_data(opt$input_ytd, opt$input_norm, opt$num_days_ytd);
+data_list = parse_input_data(opt$input_ytd, opt$input_norm, opt$num_days_ytd);
+temperature_data_frame = data_list[[1]];
+# Information needed for plots.
+start_date = data_list[[2]];
+end_date = temperature_data_frame$DATE[opt$num_days_ytd];
+start_doy_ytd = data_list[[3]];
+end_doy_ytd = data_list[[4]];
+is_leap_year = data_list[[5]];
+total_days = data_list[[6]];
+total_days_vector = c(1:total_days);
 # 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);
 }
-# Information needed for plots.
-start_date = temperature_data_frame$DATE[1];
-end_date = temperature_data_frame$DATE[opt$num_days_ytd];
-# See if we're in a leap year.
-is_leap_year = is_leap_year(start_date);
-total_days = get_total_days(is_leap_year);
-total_days_vector = c(1:total_days);
 # Get the ticks date labels for plots.
-ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, total_days, opt$num_days_ytd);
+ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, total_days, start_doy_ytd, end_doy_ytd);
 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.
 latitude = temperature_data_frame$LATITUDE[1];

Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 39:169c8180205a draft