Repository 'insect_phenology_model'
hg clone https://toolshed.g2.bx.psu.edu/repos/greg/insect_phenology_model

Changeset 49:1b6864c5b50a (2018-06-05)
Previous changeset 48:99e1c1300fcd (2018-06-05) Next changeset 50:927321ed0322 (2018-08-07)
Commit message:
Uploaded
modified:
insect_phenology_model.R
b
diff -r 99e1c1300fcd -r 1b6864c5b50a insect_phenology_model.R
--- a/insect_phenology_model.R Tue Jun 05 07:52:51 2018 -0400
+++ b/insect_phenology_model.R Tue Jun 05 07:52:59 2018 -0400
[
b'@@ -6,6 +6,7 @@\n     make_option(c("--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"),\n     make_option(c("--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"),\n     make_option(c("--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"),\n+    make_option(c("--end_date"), action="store", dest="end_date", default=NULL, help="End date for custom date interval"),\n     make_option(c("--input_norm"), action="store", dest="input_norm", help="30 year normals temperature data for selected station"),\n     make_option(c("--input_ytd"), action="store", dest="input_ytd", default=NULL, help="Year-to-date temperature data for selected location"),\n     make_option(c("--insect"), action="store", dest="insect", help="Insect name"),\n@@ -24,6 +25,7 @@\n     make_option(c("--plot_generations_separately"), action="store", dest="plot_generations_separately", help="Plot Plot P, F1 and F2 as separate lines or pool across them"),\n     make_option(c("--plot_std_error"), action="store", dest="plot_std_error", help="Plot Standard error"),\n     make_option(c("--replications"), action="store", dest="replications", type="integer", help="Number of replications"),\n+    make_option(c("--start_date"), action="store", dest="start_date", default=NULL, help="Start date for custom date interval"),\n     make_option(c("--young_nymph_accumulation"), action="store", dest="young_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (egg->young nymph)")\n )\n \n@@ -31,10 +33,11 @@\n args <- parse_args(parser, positional_arguments=TRUE);\n opt <- args$options;\n \n-add_daylight_length = function(temperature_data_frame, num_rows) {\n-    # Return a vector of daylight length (photoperido profile) for\n-    # the number of days specified in the input_ytd temperature data\n-    # (from Forsythe 1995).\n+add_daylight_length = function(temperature_data_frame) {\n+    # Return temperature_data_frame with an added column\n+    # of daylight length (photoperido profile).\n+    num_rows = dim(temperature_data_frame)[1];\n+    # From Forsythe 1995.\n     p = 0.8333;\n     latitude = temperature_data_frame$LATITUDE[1];\n     daylight_length_vector = NULL;\n@@ -63,6 +66,26 @@\n     return(data_frame);\n }\n \n+extract_date_interval_rows = function(df, start_date, end_date) {\n+    date_interval_rows = df[df$DATE >= start_date & df$DATE <= end_date];\n+    return(date_interval_rows);\n+}\n+\n+from_30_year_normals = function(norm_data_frame, start_date_doy, end_date_doy, year) {\n+    # The data we want is fully contained within the 30 year normals data.\n+    first_norm_row = which(norm_data_frame$DOY==start_date_doy);\n+    last_norm_row = which(norm_data_frame$DOY==end_date_doy);\n+    # Add 1 to the number of rows to ensure that the end date is included.\n+    tmp_data_frame_rows = last_norm_row - first_norm_row + 1;\n+    tmp_data_frame = get_new_temperature_data_frame(nrow=tmp_data_frame_rows);\n+    j = 0;\n+    for (i in first_norm_row:last_norm_row) {\n+        j = j + 1;\n+        tmp_data_frame[j,] = get_next_normals_row(norm_data_frame, year, i);\n+    }\n+    return (tmp_data_frame);\n+}\n+\n get_file_path = function(life_stage, base_name, life_stage_nymph=NULL, life_stage_adult=NULL) {\n     if (!is.null(life_stage_nymph)) {\n         lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph);\n@@ -123,7 +146,49 @@\n     return(list(p_m, p_se, f1_m, f1_se, f2_m, f2_se))\n }\n \n-get_next_normals_row = function(norm_data_frame, year, is_leap_year, index) {\n+get_new_norm_data_frame = function(is_leap_year, input_norm=NULL, nrow=0) {\n+    # The input_norm data has the following 10 columns:\n+    # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX\n+    column_names = c("STATIONID", "LATITUDE","LONGITUDE", "ELEV_M", "NAME", "ST", "MMDD", "DOY"'..b'eded for plots, some of these values are\n+# being reset here since in some case they were set above.\n start_date = data_list[[2]];\n end_date = data_list[[3]];\n-start_doy_ytd = data_list[[4]];\n-end_doy_ytd = data_list[[5]];\n+prepend_end_doy_norm = data_list[[4]];\n+append_start_doy_norm = data_list[[5]];\n is_leap_year = data_list[[6]];\n-total_days = data_list[[7]];\n-total_days_vector = c(1:total_days);\n-location =  data_list[[8]];\n+location = data_list[[7]];\n \n+if (is.null(opt$start_date) && is.null(opt$end_date)) {\n+    # We\'re plotting an entire year.\n+    date_interval = FALSE;\n+    # Display the total number of days in the Galaxy history item blurb.\n+    if (processing_year_to_date_data) {\n+        cat("Number of days year-to-date: ", opt$num_days_ytd, "\\n");\n+    } else {\n+        if (is_leap_year) {\n+            num_days = 366;\n+        } else {\n+            num_days = 365;\n+        }\n+        cat("Number of days in year: ", num_days, "\\n");\n+    }\n+} else {\n+    # FIXME: currently custom date fields are free text, but\n+    # Galaxy should soon include support for a date selector\n+    # at which point this tool should be enhanced to use it.\n+    # Validate start_date.\n+    date_interval = TRUE;\n+    # Calaculate the number of days in the date interval rather\n+    # than using the number of rows in the input temperature data.\n+    start_date = validate_date(opt$start_date);\n+    # Validate end_date.\n+    end_date = validate_date(opt$end_date);\n+    if (start_date >= end_date) {\n+        stop_err("The start date must be between 1 and 50 days before the end date when setting date intervals for plots.");\n+    }\n+    # Calculate the number of days in the date interval.\n+    num_days = difftime(end_date, start_date, units=c("days"));\n+    # Add 1 to the number of days to make the dates inclusive.  For\n+    # example, if the user enters a date range of 2018-01-01 to\n+    # 2018-01-31, they likely expect the end date to be included.\n+    num_days = num_days + 1;\n+    if (num_days > 50) {\n+        # We need to restrict date intervals since\n+        # plots render tick marks for each day.\n+        stop_err("Date intervals for plotting cannot exceed 50 days.");\n+    }\n+    # Display the total number of days in the Galaxy history item blurb.\n+    cat("Number of days in date interval: ", num_days, "\\n");\n+}\n # Create copies of the temperature data for generations P, F1 and F2 if we\'re plotting generations separately.\n if (plot_generations_separately) {\n     temperature_data_frame_P = data.frame(temperature_data_frame);\n@@ -549,7 +892,7 @@\n }\n \n # Get the ticks date labels for plots.\n-ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, total_days, start_doy_ytd, end_doy_ytd);\n+ticks_and_labels = get_x_axis_ticks_and_labels(temperature_data_frame, prepend_end_doy_norm, append_start_doy_norm, date_interval);\n ticks = c(unlist(ticks_and_labels[1]));\n date_labels = c(unlist(ticks_and_labels[2]));\n # All latitude values are the same, so get the value for plots from the first row.\n@@ -615,6 +958,7 @@\n     }\n }\n # Initialize matrices.\n+total_days = dim(temperature_data_frame)[1];\n if (process_eggs) {\n     Eggs.replications = matrix(rep(0, total_days*opt$replications), ncol=opt$replications);\n }\n@@ -778,7 +1122,7 @@\n         doy = temperature_data_frame$DOY[row];\n         # Photoperiod in the day.\n         photoperiod = temperature_data_frame$DAYLEN[row];\n-        temp.profile = get_temperature_at_hour(latitude, temperature_data_frame, row, total_days);\n+        temp.profile = get_temperature_at_hour(latitude, temperature_data_frame, row);\n         mean.temp = temp.profile[1];\n         averages.temp = temp.profile[2];\n         averages.day[row] = averages.temp;\n@@ -1459,6 +1803,7 @@\n     write.csv(temperature_data_frame_F2, file=file_path, row.names=F);\n }\n \n+total_days_vector = c(1:dim(temperature_data_frame)[1]);\n if (plot_generations_separately) {\n     for (life_stage in life_stages) {\n         if (life_stage == "Egg") {\n'