Mercurial > repos > greg > insect_phenology_model

comparison insect_phenology_model.R @ 45:315c5e1bc44a draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author greg
date Mon, 23 Apr 2018 10:19:06 -0400
parents fd3c00392fce
children 0791cca1fc5c
comparison
equal deleted inserted replaced
44:c61d3d9d44db 45:315c5e1bc44a
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", default=NULL, 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", default=NULL, type="integer", help="Total number of days in the year-to-date 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)"),
351 mortality.probability = temperature * 0.0008 + 0.03; 351 mortality.probability = temperature * 0.0008 + 0.03;
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, location) {
357 if (is.null(input_ytd)) { 357 if (is.null(input_ytd)) {
358 # We're analysing only the 30 year normals data, so create an empty 358 # We're analysing only the 30 year normals data, so create an empty
359 # data frame for containing temperature data after it is converted 359 # data frame for containing temperature data after it is converted
360 # from the 30 year normals format to the year-to-date format. 360 # from the 30 year normals format to the year-to-date format.
361 temperature_data_frame = data.frame(matrix(ncol=6, nrow=0)); 361 temperature_data_frame = data.frame(matrix(ncol=6, nrow=0));
399 # All normals data includes Feb 29 which is row 60 in 399 # All normals data includes Feb 29 which is row 60 in
400 # 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.
401 if (!is_leap_year) { 401 if (!is_leap_year) {
402 norm_data_frame = norm_data_frame[-c(60),]; 402 norm_data_frame = norm_data_frame[-c(60),];
403 } 403 }
404 # Set the location to be the station name if the user elected no to enter it.
405 if (is.null(location)) {
406 location = norm_data_frame$NAME[1];
407 }
404 if (is.null(input_ytd)) { 408 if (is.null(input_ytd)) {
405 # Convert the 30 year normals data to the year-to-date format. 409 # Convert the 30 year normals data to the year-to-date format.
406 for (i in 1:total_days) { 410 for (i in 1:total_days) {
407 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i); 411 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
408 } 412 }
426 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i); 430 temperature_data_frame[i,] = get_next_normals_row(norm_data_frame, year, is_leap_year, i);
427 } 431 }
428 } 432 }
429 # Add a column containing the daylight length for each day. 433 # Add a column containing the daylight length for each day.
430 temperature_data_frame = add_daylight_length(temperature_data_frame, total_days); 434 temperature_data_frame = add_daylight_length(temperature_data_frame, total_days);
431 return(list(temperature_data_frame, start_date, end_date, start_doy_ytd, end_doy_ytd, is_leap_year, total_days)); 435 return(list(temperature_data_frame, start_date, end_date, start_doy_ytd, end_doy_ytd, is_leap_year, total_days, location));
432 } 436 }
433 437
434 render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval, 438 render_chart = function(ticks, date_labels, chart_type, plot_std_error, insect, location, latitude, start_date, end_date, days, maxval,
435 replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL, 439 replications, life_stage, group, group_std_error, group2=NULL, group2_std_error=NULL, group3=NULL, group3_std_error=NULL,
436 life_stages_adult=NULL, life_stages_nymph=NULL) { 440 life_stages_adult=NULL, life_stages_nymph=NULL) {
523 } 527 }
524 # Display the total number of days in the Galaxy history item blurb. 528 # Display the total number of days in the Galaxy history item blurb.
525 cat("Year-to-date number of days: ", opt$num_days_ytd, "\n"); 529 cat("Year-to-date number of days: ", opt$num_days_ytd, "\n");
526 530
527 # Parse the inputs. 531 # Parse the inputs.
528 data_list = parse_input_data(opt$input_ytd, opt$input_norm, opt$num_days_ytd); 532 data_list = parse_input_data(opt$input_ytd, opt$input_norm, opt$num_days_ytd, opt$location);
529 temperature_data_frame = data_list[[1]]; 533 temperature_data_frame = data_list[[1]];
530 # Information needed for plots. 534 # Information needed for plots.
531 start_date = data_list[[2]]; 535 start_date = data_list[[2]];
532 end_date = data_list[[3]]; 536 end_date = data_list[[3]];
533 start_doy_ytd = data_list[[4]]; 537 start_doy_ytd = data_list[[4]];
534 end_doy_ytd = data_list[[5]]; 538 end_doy_ytd = data_list[[5]];
535 is_leap_year = data_list[[6]]; 539 is_leap_year = data_list[[6]];
536 total_days = data_list[[7]]; 540 total_days = data_list[[7]];
537 total_days_vector = c(1:total_days); 541 total_days_vector = c(1:total_days);
542 location = data_list[[8]];
538 543
539 # Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately. 544 # Create copies of the temperature data for generations P, F1 and F2 if we're plotting generations separately.
540 if (plot_generations_separately) { 545 if (plot_generations_separately) {
541 temperature_data_frame_P = data.frame(temperature_data_frame); 546 temperature_data_frame_P = data.frame(temperature_data_frame);
542 temperature_data_frame_F1 = data.frame(temperature_data_frame); 547 temperature_data_frame_F1 = data.frame(temperature_data_frame);
1462 file_path = get_file_path(life_stage, "egg_pop_by_generation.pdf") 1467 file_path = get_file_path(life_stage, "egg_pop_by_generation.pdf")
1463 pdf(file=file_path, width=20, height=30, bg="white"); 1468 pdf(file=file_path, width=20, height=30, bg="white");
1464 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1469 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1465 # Egg population size by generation. 1470 # Egg population size by generation.
1466 maxval = max(P_eggs+F1_eggs+F2_eggs) + 100; 1471 maxval = max(P_eggs+F1_eggs+F2_eggs) + 100;
1467 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, 1472 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
1468 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P_eggs, group_std_error=P_eggs.std_error, 1473 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P_eggs, group_std_error=P_eggs.std_error,
1469 group2=F1_eggs, group2_std_error=F1_eggs.std_error, group3=F2_eggs, group3_std_error=F2_eggs.std_error); 1474 group2=F1_eggs, group2_std_error=F1_eggs.std_error, group3=F2_eggs, group3_std_error=F2_eggs.std_error);
1470 # Turn off device driver to flush output. 1475 # Turn off device driver to flush output.
1471 dev.off(); 1476 dev.off();
1472 } else if (life_stage == "Nymph") { 1477 } else if (life_stage == "Nymph") {
1502 group2 = F1_total_nymphs; 1507 group2 = F1_total_nymphs;
1503 group2_std_error = F1_total_nymphs.std_error; 1508 group2_std_error = F1_total_nymphs.std_error;
1504 group3 = F2_total_nymphs; 1509 group3 = F2_total_nymphs;
1505 group3_std_error = F2_total_nymphs.std_error; 1510 group3_std_error = F2_total_nymphs.std_error;
1506 } 1511 }
1507 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, 1512 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
1508 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error, 1513 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1509 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_nymph=life_stage_nymph); 1514 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_nymph=life_stage_nymph);
1510 # Turn off device driver to flush output. 1515 # Turn off device driver to flush output.
1511 dev.off(); 1516 dev.off();
1512 } 1517 }
1552 group2 = F1_total_adults; 1557 group2 = F1_total_adults;
1553 group2_std_error = F1_total_adults.std_error; 1558 group2_std_error = F1_total_adults.std_error;
1554 group3 = F2_total_adults; 1559 group3 = F2_total_adults;
1555 group3_std_error = F2_total_adults.std_error; 1560 group3_std_error = F2_total_adults.std_error;
1556 } 1561 }
1557 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, 1562 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
1558 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error, 1563 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1559 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_adult=life_stage_adult); 1564 group2=group2, group2_std_error=group2_std_error, group3=group3, group3_std_error=group3_std_error, life_stages_adult=life_stage_adult);
1560 # Turn off device driver to flush output. 1565 # Turn off device driver to flush output.
1561 dev.off(); 1566 dev.off();
1562 } 1567 }
1568 file_path = get_file_path(life_stage, "total_pop_by_generation.pdf") 1573 file_path = get_file_path(life_stage, "total_pop_by_generation.pdf")
1569 pdf(file=file_path, width=20, height=30, bg="white"); 1574 pdf(file=file_path, width=20, height=30, bg="white");
1570 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1575 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1571 # Total population size by generation. 1576 # Total population size by generation.
1572 maxval = max(P+F1+F2) + 100; 1577 maxval = max(P+F1+F2) + 100;
1573 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, 1578 render_chart(ticks, date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, location, latitude,
1574 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P, group_std_error=P.std_error, 1579 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=P, group_std_error=P.std_error,
1575 group2=F1, group2_std_error=F1.std_error, group3=F2, group3_std_error=F2.std_error); 1580 group2=F1, group2_std_error=F1.std_error, group3=F2, group3_std_error=F2.std_error);
1576 # Turn off device driver to flush output. 1581 # Turn off device driver to flush output.
1577 dev.off(); 1582 dev.off();
1578 } 1583 }
1585 file_path = get_file_path(life_stage, "egg_pop.pdf") 1590 file_path = get_file_path(life_stage, "egg_pop.pdf")
1586 pdf(file=file_path, width=20, height=30, bg="white"); 1591 pdf(file=file_path, width=20, height=30, bg="white");
1587 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1592 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1588 # Egg population size. 1593 # Egg population size.
1589 maxval = max(eggs+eggs.std_error) + 100; 1594 maxval = max(eggs+eggs.std_error) + 100;
1590 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, 1595 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
1591 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error); 1596 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error);
1592 # Turn off device driver to flush output. 1597 # Turn off device driver to flush output.
1593 dev.off(); 1598 dev.off();
1594 } else if (life_stage == "Nymph") { 1599 } else if (life_stage == "Nymph") {
1595 for (life_stage_nymph in life_stages_nymph) { 1600 for (life_stage_nymph in life_stages_nymph) {
1610 # Old nymph population size. 1615 # Old nymph population size.
1611 group = old_nymphs; 1616 group = old_nymphs;
1612 group_std_error = old_nymphs.std_error; 1617 group_std_error = old_nymphs.std_error;
1613 } 1618 }
1614 maxval = max(group+group_std_error) + 100; 1619 maxval = max(group+group_std_error) + 100;
1615 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, 1620 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
1616 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error, 1621 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1617 life_stages_nymph=life_stage_nymph); 1622 life_stages_nymph=life_stage_nymph);
1618 # Turn off device driver to flush output. 1623 # Turn off device driver to flush output.
1619 dev.off(); 1624 dev.off();
1620 } 1625 }
1641 # Diapausing adult population size. 1646 # Diapausing adult population size.
1642 group = diapausing_adults; 1647 group = diapausing_adults;
1643 group_std_error = diapausing_adults.std_error 1648 group_std_error = diapausing_adults.std_error
1644 } 1649 }
1645 maxval = max(group+group_std_error) + 100; 1650 maxval = max(group+group_std_error) + 100;
1646 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, 1651 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
1647 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error, 1652 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error,
1648 life_stages_adult=life_stage_adult); 1653 life_stages_adult=life_stage_adult);
1649 # Turn off device driver to flush output. 1654 # Turn off device driver to flush output.
1650 dev.off(); 1655 dev.off();
1651 } 1656 }
1655 file_path = get_file_path(life_stage, "total_pop.pdf") 1660 file_path = get_file_path(life_stage, "total_pop.pdf")
1656 pdf(file=file_path, width=20, height=30, bg="white"); 1661 pdf(file=file_path, width=20, height=30, bg="white");
1657 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); 1662 par(mar=c(5, 6, 4, 4), mfrow=c(3, 1));
1658 # Total population size. 1663 # Total population size.
1659 maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100; 1664 maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100;
1660 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, 1665 render_chart(ticks, date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, location, latitude,
1661 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=total_adults, group_std_error=total_adults.std_error, 1666 start_date, end_date, total_days_vector, maxval, opt$replications, life_stage, group=total_adults, group_std_error=total_adults.std_error,
1662 group2=total_nymphs, group2_std_error=total_nymphs.std_error, group3=eggs, group3_std_error=eggs.std_error); 1667 group2=total_nymphs, group2_std_error=total_nymphs.std_error, group3=eggs, group3_std_error=eggs.std_error);
1663 # Turn off device driver to flush output. 1668 # Turn off device driver to flush output.
1664 dev.off(); 1669 dev.off();
1665 } 1670 }