Mercurial > repos > greg > insect_phenology_model
comparison insect_phenology_model.R @ 43:fd3c00392fce draft
Uploaded
author | greg |
---|---|
date | Mon, 23 Apr 2018 09:49:08 -0400 |
parents | f4d683709b7f |
children | 315c5e1bc44a |
comparison
equal
deleted
inserted
replaced
42:64132300c62e | 43:fd3c00392fce |
---|---|
5 option_list <- list( | 5 option_list <- list( |
6 make_option(c("--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"), | 6 make_option(c("--adult_mortality"), action="store", dest="adult_mortality", type="integer", help="Adjustment rate for adult mortality"), |
7 make_option(c("--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"), | 7 make_option(c("--adult_accumulation"), action="store", dest="adult_accumulation", type="integer", help="Adjustment of degree-days accumulation (old nymph->adult)"), |
8 make_option(c("--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"), | 8 make_option(c("--egg_mortality"), action="store", dest="egg_mortality", type="integer", help="Adjustment rate for egg mortality"), |
9 make_option(c("--input_norm"), action="store", dest="input_norm", help="30 year normals temperature data for selected station"), | 9 make_option(c("--input_norm"), action="store", dest="input_norm", help="30 year normals temperature data for selected station"), |
10 make_option(c("--input_ytd"), action="store", dest="input_ytd", help="Year-to-date temperature data for selected location"), | 10 make_option(c("--input_ytd"), action="store", dest="input_ytd", default=NULL, help="Year-to-date temperature data for selected location"), |
11 make_option(c("--insect"), action="store", dest="insect", help="Insect name"), | 11 make_option(c("--insect"), action="store", dest="insect", help="Insect name"), |
12 make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"), | 12 make_option(c("--insects_per_replication"), action="store", dest="insects_per_replication", type="integer", help="Number of insects with which to start each replication"), |
13 make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"), | 13 make_option(c("--life_stages"), action="store", dest="life_stages", help="Selected life stages for plotting"), |
14 make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"), | 14 make_option(c("--life_stages_adult"), action="store", dest="life_stages_adult", default=NULL, help="Adult life stages for plotting"), |
15 make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"), | 15 make_option(c("--life_stages_nymph"), action="store", dest="life_stages_nymph", default=NULL, help="Nymph life stages for plotting"), |
16 make_option(c("--location"), action="store", dest="location", help="Selected location"), | 16 make_option(c("--location"), action="store", dest="location", help="Selected location"), |
17 make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), | 17 make_option(c("--min_clutch_size"), action="store", dest="min_clutch_size", type="integer", help="Adjustment of minimum clutch size"), |
18 make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), | 18 make_option(c("--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), |
19 make_option(c("--num_days_ytd"), action="store", dest="num_days_ytd", type="integer", help="Total number of days in the temperature dataset"), | 19 make_option(c("--num_days_ytd"), action="store", dest="num_days_ytd", default=NULL, type="integer", help="Total number of days in the year-to-date temperature dataset"), |
20 make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"), | 20 make_option(c("--nymph_mortality"), action="store", dest="nymph_mortality", type="integer", help="Adjustment rate for nymph mortality"), |
21 make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"), | 21 make_option(c("--old_nymph_accumulation"), action="store", dest="old_nymph_accumulation", type="integer", help="Adjustment of degree-days accumulation (young nymph->old nymph)"), |
22 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), | 22 make_option(c("--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), |
23 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), | 23 make_option(c("--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), |
24 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"), | 24 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"), |
260 # that was prepended to the year-to-date data. | 260 # that was prepended to the year-to-date data. |
261 tick_index = get_tick_index(i, last_tick, ticks, month_labels) | 261 tick_index = get_tick_index(i, last_tick, ticks, month_labels) |
262 ticks[tick_index] = i; | 262 ticks[tick_index] = i; |
263 month_labels[tick_index] = "End prepended 30 year normals"; | 263 month_labels[tick_index] = "End prepended 30 year normals"; |
264 last_tick = i; | 264 last_tick = i; |
265 } else if (i==end_doy_ytd+1) { | 265 } else if (end_doy_ytd > 0 & i==end_doy_ytd+1) { |
266 # Add a tick for the start of the 30 year normnals data | 266 # Add a tick for the start of the 30 year normnals data |
267 # that was appended to the year-to-date data. | 267 # that was appended to the year-to-date data. |
268 tick_index = get_tick_index(i, last_tick, ticks, month_labels) | 268 tick_index = get_tick_index(i, last_tick, ticks, month_labels) |
269 ticks[tick_index] = i; | 269 ticks[tick_index] = i; |
270 month_labels[tick_index] = "Start appended 30 year normals"; | 270 month_labels[tick_index] = "Start appended 30 year normals"; |
352 } | 352 } |
353 return(mortality.probability); | 353 return(mortality.probability); |
354 } | 354 } |
355 | 355 |
356 parse_input_data = function(input_ytd, input_norm, num_days_ytd) { | 356 parse_input_data = function(input_ytd, input_norm, num_days_ytd) { |
357 # Read the input_ytd temperature datafile into a data frame. | 357 if (is.null(input_ytd)) { |
358 # The input_ytd data has the following 6 columns: | 358 # We're analysing only the 30 year normals data, so create an empty |
359 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX | 359 # data frame for containing temperature data after it is converted |
360 temperature_data_frame = read.csv(file=input_ytd, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=","); | 360 # from the 30 year normals format to the year-to-date format. |
361 # Set the temperature_data_frame column names for access. | 361 temperature_data_frame = data.frame(matrix(ncol=6, nrow=0)); |
362 colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX"); | 362 colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX"); |
363 # Get the start date. | 363 # Base all dates on the current date since 30 year |
364 start_date = temperature_data_frame$DATE[1]; | 364 # normals data does not include any dates. |
365 end_date = temperature_data_frame$DATE[num_days_ytd]; | 365 year = format(Sys.Date(), "%Y"); |
366 start_date = paste(year, "01", "01", sep="-"); | |
367 end_date = paste(year, "12", "31", sep="-"); | |
368 # Set invalid start and end DOY. | |
369 start_doy_ytd = 0; | |
370 end_doy_ytd = 0; | |
371 } else { | |
372 # Read the input_ytd temperature datafile into a data frame. | |
373 # The input_ytd data has the following 6 columns: | |
374 # LATITUDE, LONGITUDE, DATE, DOY, TMIN, TMAX | |
375 temperature_data_frame = read.csv(file=input_ytd, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=","); | |
376 # Set the temperature_data_frame column names for access. | |
377 colnames(temperature_data_frame) = c("LATITUDE", "LONGITUDE", "DATE", "DOY", "TMIN", "TMAX"); | |
378 # Get the start date. | |
379 start_date = temperature_data_frame$DATE[1]; | |
380 end_date = temperature_data_frame$DATE[num_days_ytd]; | |
381 # Extract the year from the start date. | |
382 date_str = format(start_date); | |
383 date_str_items = strsplit(date_str, "-")[[1]]; | |
384 year = date_str_items[1]; | |
385 # Save the first DOY to later check if start_date is Jan 1. | |
386 start_doy_ytd = as.integer(temperature_data_frame$DOY[1]); | |
387 end_doy_ytd = as.integer(temperature_data_frame$DOY[num_days_ytd]); | |
388 } | |
366 # See if we're in a leap year. | 389 # See if we're in a leap year. |
367 is_leap_year = is_leap_year(start_date); | 390 is_leap_year = is_leap_year(start_date); |
368 # Get the number of days in the year. | 391 # Get the number of days in the year. |
369 total_days = get_total_days(is_leap_year); | 392 total_days = get_total_days(is_leap_year); |
370 # Extract the year from the start date. | |
371 date_str = format(start_date); | |
372 date_str_items = strsplit(date_str, "-")[[1]]; | |
373 year = date_str_items[1]; | |
374 # Save the first DOY to later check if start_date is Jan 1. | |
375 start_doy_ytd = as.integer(temperature_data_frame$DOY[1]); | |
376 end_doy_ytd = as.integer(temperature_data_frame$DOY[num_days_ytd]); | |
377 # Read the input_norm temperature datafile into a data frame. | 393 # Read the input_norm temperature datafile into a data frame. |
378 # The input_norm data has the following 10 columns: | 394 # The input_norm data has the following 10 columns: |
379 # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX | 395 # STATIONID, LATITUDE, LONGITUDE, ELEV_M, NAME, ST, MMDD, DOY, TMIN, TMAX |
380 norm_data_frame = read.csv(file=input_norm, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=","); | 396 norm_data_frame = read.csv(file=input_norm, header=T, strip.white=TRUE, stringsAsFactors=FALSE, sep=","); |
381 # Set the norm_data_frame column names for access. | 397 # Set the norm_data_frame column names for access. |
383 # All normals data includes Feb 29 which is row 60 in | 399 # All normals data includes Feb 29 which is row 60 in |
384 # the data, so delete that row if we're not in a leap year. | 400 # the data, so delete that row if we're not in a leap year. |
385 if (!is_leap_year) { | 401 if (!is_leap_year) { |
386 norm_data_frame = norm_data_frame[-c(60),]; | 402 norm_data_frame = norm_data_frame[-c(60),]; |
387 } | 403 } |
388 if (start_doy_ytd > 1) { | 404 if (is.null(input_ytd)) { |
389 # The year-to-date data starts after Jan 1, so create a | 405 # Convert the 30 year normals data to the year-to-date format. |
390 # temporary data frame to contain the 30 year normals data | 406 for (i in 1:total_days) { |
391 # from Jan 1 to the date immediately prior to start_date. | 407 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i); |
392 tmp_data_frame = temperature_data_frame[FALSE,]; | 408 } |
393 for (i in 1:start_doy_ytd-1) { | 409 } else { |
394 tmp_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i); | 410 # Merge the year-to-date data with the 30 year normals data. |
395 } | 411 if (start_doy_ytd > 1) { |
396 # Next merge the temporary data frame with the year-to-date data frame. | 412 # The year-to-date data starts after Jan 1, so create a |
397 temperature_data_frame = rbind(tmp_data_frame, temperature_data_frame); | 413 # temporary data frame to contain the 30 year normals data |
398 } | 414 # from Jan 1 to the date immediately prior to start_date. |
399 # Define the next row for the year-to-date data from the 30 year normals data. | 415 tmp_data_frame = temperature_data_frame[FALSE,]; |
400 first_normals_append_row = end_doy_ytd + 1; | 416 for (i in 1:start_doy_ytd-1) { |
401 # Append the 30 year normals data to the year-to-date data. | 417 tmp_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i); |
402 for (i in first_normals_append_row:total_days) { | 418 } |
403 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i); | 419 # Next merge the temporary data frame with the year-to-date data frame. |
420 temperature_data_frame = rbind(tmp_data_frame, temperature_data_frame); | |
421 } | |
422 # Define the next row for the year-to-date data from the 30 year normals data. | |
423 first_normals_append_row = end_doy_ytd + 1; | |
424 # Append the 30 year normals data to the year-to-date data. | |
425 for (i in first_normals_append_row:total_days) { | |
426 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i); | |
427 } | |
404 } | 428 } |
405 # Add a column containing the daylight length for each day. | 429 # Add a column containing the daylight length for each day. |
406 temperature_data_frame = add_daylight_length(temperature_data_frame, total_days); | 430 temperature_data_frame = add_daylight_length(temperature_data_frame, total_days); |
407 return(list(temperature_data_frame, start_date, end_date, start_doy_ytd, end_doy_ytd, is_leap_year, total_days)); | 431 return(list(temperature_data_frame, start_date, end_date, start_doy_ytd, end_doy_ytd, is_leap_year, total_days)); |
408 } | 432 } |