Mercurial > repos > greg > insect_phenology_model
changeset 18:f5ecff4800f2 draft
Uploaded
author | greg |
---|---|
date | Tue, 06 Mar 2018 13:39:14 -0500 |
parents | 6f31ade9a0f4 |
children | 3c6e94e477cb |
files | insect_phenology_model.R insect_phenology_model.xml test-data/pop_by_generation.pdf test-data/total_pop_by_generation.pdf |
diffstat | 4 files changed, 161 insertions(+), 105 deletions(-) [+] |
line wrap: on
line diff
--- a/insect_phenology_model.R Tue Mar 06 08:26:50 2018 -0500 +++ b/insect_phenology_model.R Tue Mar 06 13:39:14 2018 -0500 @@ -53,16 +53,22 @@ } dev.egg = function(temperature) { + # This function is not used, but was + # in the original, so keep it for now. dev.rate = -0.9843 * temperature + 33.438; return(dev.rate); } dev.emerg = function(temperature) { + # This function is not used, but was + # in the original, so keep it for now. emerg.rate = -0.5332 * temperature + 24.147; return(emerg.rate); } dev.old = function(temperature) { + # This function is not used, but was + # in the original, so keep it for now. n34 = -0.6119 * temperature + 17.602; n45 = -0.4408 * temperature + 19.036; dev.rate = mean(n34 + n45); @@ -70,6 +76,8 @@ } dev.young = function(temperature) { + # This function is not used, but was + # in the original, so keep it for now. n12 = -0.3728 * temperature + 14.68; n23 = -0.6119 * temperature + 25.249; dev.rate = mean(n12 + n23); @@ -95,6 +103,35 @@ return(c(unlist(month_labels))); } +get_life_stage_index = function(life_stage, life_stage_nymph=NULL, life_stage_adult=NULL) { + # Name collection elements so that they + # are displayed in logical order. + if (life_stage=="Egg") { + lsi = "01"; + } else if (life_stage=="Nymph") { + if (life_stage_nymph=="Young") { + lsi = "02"; + } else if (life_stage_nymph=="Old") { + lsi = "03"; + } else if (life_stage_nymph=="Total") { + lsi="04"; + } + } else if (life_stage=="Adult") { + if (life_stage_adult=="Pre-vittelogenic") { + lsi = "05"; + } else if (life_stage_adult=="Vittelogenic") { + lsi = "06"; + } else if (life_stage_adult=="Diapausing") { + lsi = "07"; + } else if (life_stage_adult=="Total") { + lsi = "08"; + } + } else if (life_stage=="Total") { + lsi = "09"; + } + return(lsi); +} + get_temperature_at_hour = function(latitude, temperature_data_frame, row, num_days) { # Base development threshold for Brown Marmorated Stink Bug # insect phenology model. @@ -218,8 +255,8 @@ render_chart = function(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) { + 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=" "); @@ -389,7 +426,7 @@ population.replications = matrix(rep(0, opt$num_days*opt$replications), ncol=opt$replications); # Process replications. -for (N.replications in 1:opt$replications) { +for (current_replication in 1:opt$replications) { # Start with the user-defined number of insects per replication. num_insects = opt$insects_per_replication; # Generation, Stage, degree-days, T, Diapause. @@ -474,6 +511,7 @@ death.probability = opt$egg_mortality * mortality.egg(mean.temp); } else if (vector.individual[2] == 1 | vector.individual[2] == 2) { + # Nymph. death.probability = opt$nymph_mortality * mortality.nymph(mean.temp); } else if (vector.individual[2] == 3 | vector.individual[2] == 4 | vector.individual[2] == 5) { @@ -711,7 +749,7 @@ if (plot_generations_separately) { if (process_eggs) { - # For egg life stage of generation F1 population size, + # For egg life stage of generation P population size, # column 1 (generation) is 0 and column 2 (Stage) is 0. P.egg[row] = sum(vector.matrix[,1]==0 & vector.matrix[,2]==0); # For egg life stage of generation F1 population size, @@ -722,7 +760,7 @@ F2.egg[row] = sum(vector.matrix[,1]==2 & vector.matrix[,2]==0); } if (process_nymphs) { - # For nymph life stage of generation F1 population + # For nymph life stage of generation P population # size, one of the following combinations is required: # - column 1 (Generation) is 0 and column 2 (Stage) is 1 (Young nymph) # - column 1 (Generation) is 0 and column 2 (Stage) is 2 (Old nymph) @@ -765,44 +803,45 @@ # Define the output values. if (process_eggs) { - Eggs.replications[,N.replications] = Eggs; + Eggs.replications[,current_replication] = Eggs; } if (process_nymphs) { - YoungNymphs.replications[,N.replications] = YoungNymphs; - OldNymphs.replications[,N.replications] = OldNymphs; + YoungNymphs.replications[,current_replication] = YoungNymphs; + OldNymphs.replications[,current_replication] = OldNymphs; } if (process_adults) { - Previtellogenic.replications[,N.replications] = Previtellogenic; - Vitellogenic.replications[,N.replications] = Vitellogenic; - Diapausing.replications[,N.replications] = Diapausing; + Previtellogenic.replications[,current_replication] = Previtellogenic; + Vitellogenic.replications[,current_replication] = Vitellogenic; + Diapausing.replications[,current_replication] = Diapausing; } - newborn.replications[,N.replications] = N.newborn; - adult.replications[,N.replications] = N.adult; - death.replications[,N.replications] = N.death; + newborn.replications[,current_replication] = N.newborn; + adult.replications[,current_replication] = N.adult; + death.replications[,current_replication] = N.death; if (plot_generations_separately) { # P is Parental, or overwintered adults. - P.replications[,N.replications] = overwintering_adult.population; + P.replications[,current_replication] = overwintering_adult.population; # F1 is the first field-produced generation. - F1.replications[,N.replications] = first_generation.population; + F1.replications[,current_replication] = first_generation.population; # F2 is the second field-produced generation. - F2.replications[,N.replications] = second_generation.population; + F2.replications[,current_replication] = second_generation.population; if (process_eggs) { - P_eggs.replications[,N.replications] = P.egg; - F1_eggs.replications[,N.replications] = F1.egg; - F2_eggs.replications[,N.replications] = F2.egg; + P_eggs.replications[,current_replication] = P.egg; + F1_eggs.replications[,current_replication] = F1.egg; + F2_eggs.replications[,current_replication] = F2.egg; } if (process_nymphs) { - P_nymphs.replications[,N.replications] = P.nymph; - F1_nymphs.replications[,N.replications] = F1.nymph; - F2_nymphs.replications[,N.replications] = F2.nymph; + P_nymphs.replications[,current_replication] = P.nymph; + F1_nymphs.replications[,current_replication] = F1.nymph; + F2_nymphs.replications[,current_replication] = F2.nymph; } if (process_adults) { - P_adults.replications[,N.replications] = P.adult; - F1_adults.replications[,N.replications] = F1.adult; - F2_adults.replications[,N.replications] = F2.adult; + P_adults.replications[,current_replication] = P.adult; + F1_adults.replications[,current_replication] = F1.adult; + F2_adults.replications[,current_replication] = F2.adult; } } - population.replications[,N.replications] = total.population; + population.replications[,current_replication] = total.population; + # End processing replications. } if (process_eggs) { @@ -814,17 +853,17 @@ if (process_nymphs) { # Calculate nymph populations for selected life stage. for (life_stage_nymph in life_stages_nymph) { - if (life_stages_nymph=="Total") { + if (life_stage_nymph=="Total") { # Mean value for all nymphs. total_nymphs = apply((YoungNymphs.replications+OldNymphs.replications), 1, mean); # Standard error for all nymphs. total_nymphs.std_error = apply((YoungNymphs.replications+OldNymphs.replications) / sqrt(opt$replications), 1, sd); - } else if (life_stages_nymph=="Young") { + } else if (life_stage_nymph=="Young") { # Mean value for young nymphs. young_nymphs = apply(YoungNymphs.replications, 1, mean); # Standard error for young nymphs. young_nymphs.std_error = apply(YoungNymphs.replications / sqrt(opt$replications), 1, sd); - } else if (life_stages_nymph=="Old") { + } else if (life_stage_nymph=="Old") { # Mean value for old nymphs. old_nymphs = apply(OldNymphs.replications, 1, mean); # Standard error for old nymphs. @@ -835,22 +874,22 @@ if (process_adults) { # Calculate adult populations for selected life stage. for (life_stage_adult in life_stages_adult) { - if (life_stages_adult=="Total") { + if (life_stage_adult=="Total") { # Mean value for all adults. total_adults = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, mean); # Standard error for all adults. total_adults.std_error = apply((Previtellogenic.replications+Vitellogenic.replications+Diapausing.replications), 1, sd) / sqrt(opt$replications); - } else if (life_stages_adult == "Pre-vittelogenic") { + } else if (life_stage_adult == "Pre-vittelogenic") { # Mean value for previtellogenic adults. previttelogenic_adults = apply(Previtellogenic.replications, 1, mean); # Standard error for previtellogenic adults. previttelogenic_adults.std_error = apply(Previtellogenic.replications, 1, sd) / sqrt(opt$replications); - } else if (life_stages_adult == "Vittelogenic") { + } else if (life_stage_adult == "Vittelogenic") { # Mean value for vitellogenic adults. vittelogenic_adults = apply(Vitellogenic.replications, 1, mean); # Standard error for vitellogenic adults. vittelogenic_adults.std_error = apply(Vitellogenic.replications, 1, sd) / sqrt(opt$replications); - } else if (life_stages_adult == "Diapausing") { + } else if (life_stage_adult == "Diapausing") { # Mean value for vitellogenic adults. diapausing_adults = apply(Diapausing.replications, 1, mean); # Standard error for vitellogenic adults. @@ -929,11 +968,13 @@ if (life_stage == "Egg") { # Start PDF device driver. dev.new(width=20, height=30); - file_path = paste(output_dir, "egg_pop_by_generation.pdf", sep="/"); + lsi = get_life_stage_index(life_stage); + file_name = paste(lsi, "egg_pop_by_generation.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); pdf(file=file_path, width=20, height=30, bg="white"); par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); # Egg population size by generation. - maxval = max(F2_eggs) + 100; + maxval = max(P_eggs+F1_eggs+F2_eggs) + 100; render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, opt$replications, life_stage, group=P_eggs, group_std_error=P_eggs.std_error, group2=F1_eggs, group2_std_error=F1_eggs.std_error, group3=F2_eggs, group3_std_error=F2_eggs.std_error); @@ -943,15 +984,16 @@ for (life_stage_nymph in life_stages_nymph) { # Start PDF device driver. dev.new(width=20, height=30); - file_name = paste(tolower(life_stage_nymph), "nymph_pop_by_generation.pdf", sep="_"); + lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph); + file_name = paste(lsi, tolower(life_stage_nymph), "nymph_pop_by_generation.pdf", sep="_"); file_path = paste(output_dir, file_name, sep="/"); pdf(file=file_path, width=20, height=30, bg="white"); par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); # Nymph population size by generation. - maxval = max(F2_nymphs) + 100; + maxval = max(P_nymphs+F1_nymphs+F2_nymphs) + 100; render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, - opt$replications, life_stage, group=P_nymphs, group_std_error=P_nymphs.std_error, group2=F1_nymphs, group2_std_error=F1_nymphs.std_error, - group3=F2_nymphs, group3_std_error=F2_nymphs.std_error, life_stages_nymph=life_stage_nymph); + opt$replications, life_stage, group=P_nymphs, group_std_error=P_nymphs.std_error, group2=F1_nymphs, group2_std_error=F1_nymphs.std_error, + group3=F2_nymphs, group3_std_error=F2_nymphs.std_error, life_stages_nymph=life_stage_nymph); # Turn off device driver to flush output. dev.off(); } @@ -959,12 +1001,13 @@ for (life_stage_adult in life_stages_adult) { # Start PDF device driver. dev.new(width=20, height=30); - file_name = paste(tolower(life_stage_adult), "adult_pop_by_generation.pdf", sep="_"); + lsi = get_life_stage_index(life_stage, life_stage_adult=life_stage_adult); + file_name = paste(lsi, tolower(life_stage_adult), "adult_pop_by_generation.pdf", sep="_"); file_path = paste(output_dir, file_name, sep="/"); pdf(file=file_path, width=20, height=30, bg="white"); par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); # Adult population size by generation. - maxval = max(F2_adults) + 100; + maxval = max(P_adults+F1_adults+F2_adults) + 100; render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, opt$replications, life_stage, group=P_adults, group_std_error=P_adults.std_error, group2=F1_adults, group2_std_error=F1_adults.std_error, group3=F2_adults, group3_std_error=F2_adults.std_error, life_stages_adult=life_stage_adult); @@ -973,12 +1016,16 @@ } } else if (life_stage == "Total") { # Start PDF device driver. + # Name collection elements so that they + # are displayed in logical order. dev.new(width=20, height=30); - file_path = paste(output_dir, "total_pop_by_generation.pdf", sep="/"); + lsi = get_life_stage_index(life_stage); + file_name = paste(lsi, "total_pop_by_generation.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); pdf(file=file_path, width=20, height=30, bg="white"); par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); # Total population size by generation. - maxval = max(F2); + maxval = max(P+F1+F2) + 100; render_chart(date_labels, "pop_size_by_generation", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, opt$replications, life_stage, group=P, group_std_error=P.std_error, group2=F1, group2_std_error=F1.std_error, group3=F2, group3_std_error=F2.std_error); # Turn off device driver to flush output. @@ -990,11 +1037,13 @@ if (life_stage == "Egg") { # Start PDF device driver. dev.new(width=20, height=30); - file_path = paste(output_dir, "egg_pop.pdf", sep="/"); + lsi = get_life_stage_index(life_stage); + file_name = paste(lsi, "egg_pop.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); pdf(file=file_path, width=20, height=30, bg="white"); par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); # Egg population size. - maxval = max(eggs+eggs.std_error); + maxval = max(eggs+eggs.std_error) + 100; render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, opt$replications, life_stage, group=eggs, group_std_error=eggs.std_error); # Turn off device driver to flush output. @@ -1003,7 +1052,8 @@ for (life_stage_nymph in life_stages_nymph) { # Start PDF device driver. dev.new(width=20, height=30); - file_name = paste(tolower(life_stage_nymph), "nymph_pop.pdf", sep="_"); + lsi = get_life_stage_index(life_stage, life_stage_nymph=life_stage_nymph); + file_name = paste(lsi, tolower(life_stage_nymph), "nymph_pop.pdf", sep="_"); file_path = paste(output_dir, file_name, sep="/"); pdf(file=file_path, width=20, height=30, bg="white"); par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); @@ -1020,7 +1070,7 @@ group = old_nymphs; group_std_error = old_nymphs.std_error; } - maxval = max(group+group.std_error); + maxval = max(group+group_std_error) + 100; render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error, life_stages_nymph=life_stage_nymph); # Turn off device driver to flush output. @@ -1030,7 +1080,8 @@ for (life_stage_adult in life_stages_adult) { # Start PDF device driver. dev.new(width=20, height=30); - file_name = paste(tolower(life_stages_adult), "adult_pop.pdf", sep="_"); + lsi = get_life_stage_index(life_stage, life_stage_adult=life_stage_adult); + file_name = paste(lsi, tolower(life_stage_adult), "adult_pop.pdf", sep="_"); file_path = paste(output_dir, file_name, sep="/"); pdf(file=file_path, width=20, height=30, bg="white"); par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); @@ -1051,7 +1102,7 @@ group = diapausing_adults; group_std_error = diapausing_adults.std_error } - maxval = max(group+group_std_error); + maxval = max(group+group_std_error) + 100; render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, opt$replications, life_stage, group=group, group_std_error=group_std_error, life_stages_adult=life_stage_adult); # Turn off device driver to flush output. @@ -1060,11 +1111,13 @@ } else if (life_stage == "Total") { # Start PDF device driver. dev.new(width=20, height=30); - file_path = paste(output_dir, "total_pop.pdf", sep="/"); + lsi = get_life_stage_index(life_stage); + file_name = paste(lsi, "total_pop.pdf", sep="_"); + file_path = paste(output_dir, file_name, sep="/"); pdf(file=file_path, width=20, height=30, bg="white"); par(mar=c(5, 6, 4, 4), mfrow=c(3, 1)); # Total population size. - maxval = max(eggs+eggs.std_error, nymphs+nymphs.std_error, adults+adults.std_error); + maxval = max(eggs+eggs.std_error, total_nymphs+total_nymphs.std_error, total_adults+total_adults.std_error) + 100; render_chart(date_labels, "pop_size_by_life_stage", opt$plot_std_error, opt$insect, opt$location, latitude, start_date, end_date, days, maxval, opt$replications, life_stage, group=total_adults, group_std_error=total_adults.std_error, group2=total_nymphs, group2_std_error=total_nymphs.std_error, group3=eggs, group3_std_error=eggs.std_error);
--- a/insect_phenology_model.xml Tue Mar 06 08:26:50 2018 -0500 +++ b/insect_phenology_model.xml Tue Mar 06 13:39:14 2018 -0500 @@ -118,7 +118,10 @@ <param name="location" value="State College PA" /> <param name="replications" value="2" /> <output_collection name="output_collection" type="list"> - <element name="total_pop_by_generation.pdf" file="total_pop_by_generation.pdf" ftype="pdf" compare="contains"/> + <element name="01_egg_pop_by_generation.pdf" file="pop_by_generation.pdf" ftype="pdf" compare="contains"/> + <element name="04_total_nymph_pop_by_generation.pdf" file="pop_by_generation.pdf" ftype="pdf" compare="contains"/> + <element name="08_total_adult_pop_by_generation.pdf" file="pop_by_generation.pdf" ftype="pdf" compare="contains"/> + <element name="09_total_pop_by_generation.pdf" file="pop_by_generation.pdf" ftype="pdf" compare="contains"/> </output_collection> </test> </tests>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/pop_by_generation.pdf Tue Mar 06 13:39:14 2018 -0500 @@ -0,0 +1,52 @@ +%PDF +1 0 obj +<< +/CreationDate +/ModDate +/Title (R Graphics Output) +/Producer +/Creator (R) +<< /Type /Catalog /Pages 3 0 R >> +endobj +7 0 obj +<< /Type /Page /Parent 3 0 R /Contents 8 0 R /Resources 4 0 R >> +endobj +3 0 obj +<< /Type /Pages /Kids [ 7 0 R ] /Count 1 /MediaBox [0 0 1440 2160] >> +endobj +4 0 obj +<< +/ProcSet [/PDF /Text] +/Font <</F2 10 0 R /F3 11 0 R >> +/ExtGState << >> +/ColorSpace << /sRGB 5 0 R >> +>> +endobj +5 0 obj +[/ICCBased 6 0 R] +endobj +6 0 obj +<< /Alternate /DeviceRGB /N 3 /Length 2596 /Filter /FlateDecode >> +stream +endobj +9 0 obj +<< +/Type /Encoding /BaseEncoding /WinAnsiEncoding +/Differences [ 45/minus 96/quoteleft +144/dotlessi /grave /acute /circumflex /tilde /macron /breve /dotaccent +/dieresis /.notdef /ring /cedilla /.notdef /hungarumlaut /ogonek /caron /space] +>> +endobj +10 0 obj +<< /Type /Font /Subtype /Type1 /Name /F2 /BaseFont /Helvetica +/Encoding 9 0 R >> +endobj +11 0 obj +<< /Type /Font /Subtype /Type1 /Name /F3 /BaseFont /Helvetica-Bold +/Encoding 9 0 R >> +endobj +xref +trailer +<< /Size 12 /Info 1 0 R /Root 2 0 R >> +startxref +%%EOF
--- a/test-data/total_pop_by_generation.pdf Tue Mar 06 08:26:50 2018 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,52 +0,0 @@ -%PDF -1 0 obj -<< -/CreationDate -/ModDate -/Title (R Graphics Output) -/Producer -/Creator (R) -<< /Type /Catalog /Pages 3 0 R >> -endobj -7 0 obj -<< /Type /Page /Parent 3 0 R /Contents 8 0 R /Resources 4 0 R >> -endobj -3 0 obj -<< /Type /Pages /Kids [ 7 0 R ] /Count 1 /MediaBox [0 0 1440 2160] >> -endobj -4 0 obj -<< -/ProcSet [/PDF /Text] -/Font <</F2 10 0 R /F3 11 0 R >> -/ExtGState << >> -/ColorSpace << /sRGB 5 0 R >> ->> -endobj -5 0 obj -[/ICCBased 6 0 R] -endobj -6 0 obj -<< /Alternate /DeviceRGB /N 3 /Length 2596 /Filter /FlateDecode >> -stream -endobj -9 0 obj -<< -/Type /Encoding /BaseEncoding /WinAnsiEncoding -/Differences [ 45/minus 96/quoteleft -144/dotlessi /grave /acute /circumflex /tilde /macron /breve /dotaccent -/dieresis /.notdef /ring /cedilla /.notdef /hungarumlaut /ogonek /caron /space] ->> -endobj -10 0 obj -<< /Type /Font /Subtype /Type1 /Name /F2 /BaseFont /Helvetica -/Encoding 9 0 R >> -endobj -11 0 obj -<< /Type /Font /Subtype /Type1 /Name /F3 /BaseFont /Helvetica-Bold -/Encoding 9 0 R >> -endobj -xref -trailer -<< /Size 12 /Info 1 0 R /Root 2 0 R >> -startxref -%%EOF