| Previous changeset 12:e82ea5f1a179 (2017-08-14) Next changeset 14:6c78e618a7c0 (2017-09-28) |
|
Commit message:
Uploaded |
|
modified:
insect_phenology_model.R insect_phenology_model.xml |
|
added:
tool-data/ghcnd_stations.txt.sample |
|
removed:
tool-data/locations.txt.sample tool-data/years.txt.sample |
| b |
| diff -r e82ea5f1a179 -r 80d38a26b2e8 insect_phenology_model.R --- a/insect_phenology_model.R Mon Aug 14 13:44:45 2017 -0400 +++ b/insect_phenology_model.R Thu Sep 28 10:45:53 2017 -0400 |
| [ |
| @@ -12,6 +12,7 @@ make_option(c("-i", "--max_clutch_size"), action="store", dest="max_clutch_size", type="integer", help="Adjustment of maximum clutch size"), make_option(c("-j", "--nymph_mort"), action="store", dest="nymph_mort", type="integer", help="Adjustment rate for nymph mortality"), make_option(c("-k", "--old_nymph_accum"), action="store", dest="old_nymph_accum", type="integer", help="Adjustment of DD accumulation (young nymph->old nymph)"), + make_option(c("-n", "--num_days"), action="store", dest="num_days", type="integer", help="Total number of days in the temperature dataset"), make_option(c("-o", "--output"), action="store", dest="output", help="Output dataset"), make_option(c("-p", "--oviposition"), action="store", dest="oviposition", type="integer", help="Adjustment for oviposition rate"), make_option(c("-q", "--photoperiod"), action="store", dest="photoperiod", type="double", help="Critical photoperiod for diapause induction/termination"), @@ -28,36 +29,39 @@ data.input=function(loc, year, temperature.dataset) { - expdata <- matrix(rep(0, 365 * 3), nrow=365) + expdata <- matrix(rep(0, opt$num_days * 3), nrow=opt$num_days) namedat <- paste(loc, year, ".Rdat", sep="") temp.data <- read.csv(file=temperature.dataset, header=T) - expdata[,1] <- c(1:365) + expdata[,1] <- c(1:opt$num_days) # Minimum - expdata[,2] <- temp.data[c(1:365), 3] + expdata[,2] <- temp.data[c(1:opt$num_days), 3] # Maximum - expdata[,3] <- temp.data[c(1:365), 2] + expdata[,3] <- temp.data[c(1:opt$num_days), 2] save(expdata, file=namedat) namedat } -daylength=function(latitude) +daylength=function(latitude, num_days) { # from Forsythe 1995 p=0.8333 dl <- NULL - for (i in 1:365) { + for (i in 1:num_days) { theta <- 0.2163108 + 2 * atan(0.9671396 * tan(0.00860 * (i - 186))) phi <- asin(0.39795 * cos(theta)) dl[i] <- 24 - 24 / pi * acos((sin(p * pi / 180) + sin(latitude * pi / 180) * sin(phi)) / (cos(latitude * pi / 180) * cos(phi))) } - dl # return a vector of daylength in 365 days + # return a vector of daylength for the number of + # days specifie din the input temperature data + dl } -hourtemp=function(latitude, date, temperature_file_path) +hourtemp=function(latitude, date, temperature_file_path, num_days) { load(temperature_file_path) - threshold <- 14.17 # base development threshold for BMSB + # base development threshold for BMSB + threshold <- 14.17 dnp <- expdata[date, 2] # daily minimum dxp <- expdata[date, 3] # daily maximum dmean <- 0.5 * (dnp + dxp) @@ -67,22 +71,28 @@ dd <- 0 } else { - dlprofile <- daylength(latitude) # extract daylength data for entire year + # extract daylength data for the number of + # days specified in the input temperature data + dlprofile <- daylength(latitude, num_days) T <- NULL # initialize hourly temperature dh <- NULL #initialize degree hour vector - # date <- 200 - y <- dlprofile[date] # calculate daylength in given date - z <- 24 - y # night length - a <- 1.86 # lag coefficient - b <- 2.20 # night coefficient - #tempdata <- read.csv("tempdata.csv") #import raw data set - # Should be outside function otherwise its redundant - risetime <- 12 - y / 2 # sunrise time - settime <- 12 + y / 2 # sunset time + # calculate daylength in given date + y <- dlprofile[date] + # night length + z <- 24 - y + # lag coefficient + a <- 1.86 + # night coefficient + b <- 2.20 + # sunrise time + risetime <- 12 - y / 2 + # sunset time + settime <- 12 + y / 2 ts <- (dxp - dnp) * sin(pi * (settime - 5) / (y + 2 * a)) + dnp for (i in 1:24) { if (i > risetime && i<settime) { - m <- i - 5 # number of hours after Tmin until sunset + # number of hours after Tmin until sunset + m <- i - 5 T[i]=(dxp - dnp) * sin(pi * m / (y + 2 * a)) + dnp if (T[i]<8.4) { dh[i] <- 0 @@ -205,8 +215,8 @@ temperature_file_path <- data.input(opt$location, opt$year, opt$temperature_dataset) # Initialize matrix for results from all replications -S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, 365 * opt$replications), ncol = opt$replications) -newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, 365 * opt$replications), ncol=opt$replications) +S0.rep <- S1.rep <- S2.rep <- S3.rep <- S4.rep <- S5.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol = opt$replications) +newborn.rep <- death.rep <- adult.rep <- pop.rep <- g0.rep <- g1.rep <- g2.rep <- g0a.rep <- g1a.rep <- g2a.rep <- matrix(rep(0, opt$num_days * opt$replications), ncol=opt$replications) # loop through replications for (N.rep in 1:opt$replications) { @@ -220,27 +230,27 @@ # complete matrix for the population vec.mat <- t(matrix(vec.mat, nrow=5)) # complete photoperiod profile in a year, requires daylength function - ph.p <- daylength(opt$latitude) + ph.p <- daylength(opt$latitude, opt$num_days) # time series of population size tot.pop <- NULL # gen.0 pop size - gen0.pop <- rep(0, 365) - gen1.pop <- rep(0, 365) - gen2.pop <- rep(0, 365) - S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, 365) - g0.adult <- g1.adult <- g2.adult <- rep(0, 365) - N.newborn <- N.death <- N.adult <- rep(0, 365) - dd.day <- rep(0, 365) + gen0.pop <- rep(0, opt$num_days) + gen1.pop <- rep(0, opt$num_days) + gen2.pop <- rep(0, opt$num_days) + S0 <- S1 <- S2 <- S3 <- S4 <- S5 <- rep(0, opt$num_days) + g0.adult <- g1.adult <- g2.adult <- rep(0, opt$num_days) + N.newborn <- N.death <- N.adult <- rep(0, opt$num_days) + dd.day <- rep(0, opt$num_days) # start tick ptm <- proc.time() # all the days - for (day in 1:365) { + for (day in 1:opt$num_days) { # photoperiod in the day photoperiod <- ph.p[day] - temp.profile <- hourtemp(opt$latitude, day, temperature_file_path) + temp.profile <- hourtemp(opt$latitude, day, temperature_file_path, opt$num_days) mean.temp <- temp.profile[1] dd.temp <- temp.profile[2] dd.day[day] <- dd.temp @@ -514,8 +524,7 @@ N.newborn[day] <- n.newborn N.death[day] <- n.death N.adult[day] <- n.adult - #print(c(N.rep, day, n, n.adult)) - } # end of 365 days + } # end of days specified in the input temperature data dd.cum <- cumsum(dd.day) # collect all the outputs @@ -537,15 +546,11 @@ g2a.rep[,N.rep] <- g2.adult } -# save(dd.day, dd.cum, S0.rep, S1.rep, S2.rep, S3.rep, S4.rep, S5.rep, newborn.rep, death.rep, adult.rep, pop.rep, g0.rep, g1.rep, g2.rep, g0a.rep, g1a.rep, g2a.rep, file=opt$output) -# maybe do not need to export this bit, but for now just leave it as-is -# do we need to export this Rdat file? - # Data analysis and visualization # default: plot 1 year of result # but can be expanded to accommodate multiple years n.yr <- 1 -day.all <- c(1:365 * n.yr) +day.all <- c(1:opt$num_days * n.yr) # mean value for adults sa <- apply((S3.rep + S4.rep + S5.rep), 1, mean) |
| b |
| diff -r e82ea5f1a179 -r 80d38a26b2e8 insect_phenology_model.xml --- a/insect_phenology_model.xml Mon Aug 14 13:44:45 2017 -0400 +++ b/insect_phenology_model.xml Thu Sep 28 10:45:53 2017 -0400 |
| b |
| @@ -9,12 +9,13 @@ -a $adult_mort -b $adult_accum -c $egg_mort --d $latitude --e '$location' +-d $location.fields.latitude +-e '$location.fields.name' -f $max_clutch_size -i $min_clutch_size -j $nymph_mort -k $old_nymph_accum +-n $temperature_data.metadata.data_lines -o '$output' -p $oviposition -q $photoperiod @@ -25,8 +26,13 @@ -y $young_nymph_accum ]]></command> <inputs> - <param name="location" type="text" value="" optional="false" label="Location" /> - <param name="latitude" type="float" value="0.0" label="Latitude of selected location" /> + <param name="location" type="select" label="Global Historical Climatology Network Station"> + <options from_file="ghcnd_stations.txt"> + <column name="name" index="4"/> + <column name="station" index="0"/> + <column name="latitude" index="1"/> + </options> + </param> <param name="temperature_data" type="data" format="csv" label="Temperature data" /> <param name="year" type="integer" value="2017" min="1995" label="Temperature data year" /> <param name="replications" type="integer" value="10" min="1" label="Number of replications" /> @@ -46,7 +52,7 @@ </param> </inputs> <outputs> - <data name="output" format="pdf" label="${tool.name} ${location}, ${year} lat:${latitude} on ${on_string}" /> + <data name="output" format="pdf" label="${tool.name} ${location}, ${year} on ${on_string}" /> </outputs> <tests> <test> @@ -85,5 +91,7 @@ </help> <citations> <citation type="doi">10.3389/fphys.2016.00165</citation> + <citation type="doi">10.1175/JTECH-D-11-00103.22</citation> + <citation type="doi">10.7289/V5D21VHZ</citation> </citations> </tool> |
| b |
| diff -r e82ea5f1a179 -r 80d38a26b2e8 tool-data/ghcnd_stations.txt.sample --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool-data/ghcnd_stations.txt.sample Thu Sep 28 10:45:53 2017 -0400 |
| b |
| b'@@ -0,0 +1,103746 @@\n+IV. FORMAT OF "ghcnd-stations.txt"\n+\n+------------------------------\n+Variable Columns Type\n+------------------------------\n+ID 1-11 Character\n+LATITUDE 13-20 Real\n+LONGITUDE 22-30 Real\n+ELEVATION 32-37 Real\n+STATE 39-40 Character\n+NAME 42-71 Character\n+GSN FLAG 73-75 Character\n+HCN/CRN FLAG 77-79 Character\n+WMO ID 81-85 Character\n+------------------------------\n+\n+These variables have the following definitions:\n+\n+ID is the station identification code. Note that the first two\n+ characters denote the FIPS country code, the third character \n+ is a network code that identifies the station numbering system \n+ used, and the remaining eight characters contain the actual \n+ station ID. \n+\n+ See "ghcnd-countries.txt" for a complete list of country codes.\n+\t See "ghcnd-states.txt" for a list of state/province/territory codes.\n+\n+ The network code has the following five values:\n+\n+ 0 = unspecified (station identified by up to eight \n+\t alphanumeric characters)\n+\t 1 = Community Collaborative Rain, Hail,and Snow (CoCoRaHS)\n+\t based identification number. To ensure consistency with\n+\t with GHCN Daily, all numbers in the original CoCoRaHS IDs\n+\t have been left-filled to make them all four digits long. \n+\t In addition, the characters "-" and "_" have been removed \n+\t to ensure that the IDs do not exceed 11 characters when \n+\t preceded by "US1". For example, the CoCoRaHS ID \n+\t "AZ-MR-156" becomes "US1AZMR0156" in GHCN-Daily\n+ C = U.S. Cooperative Network identification number (last six \n+ characters of the GHCN-Daily ID)\n+\t E = Identification number used in the ECA&D non-blended\n+\t dataset\n+\t M = World Meteorological Organization ID (last five\n+\t characters of the GHCN-Daily ID)\n+\t N = Identification number used in data supplied by a \n+\t National Meteorological or Hydrological Center\n+\t R = U.S. Interagency Remote Automatic Weather Station (RAWS)\n+\t identifier\n+\t S = U.S. Natural Resources Conservation Service SNOwpack\n+\t TELemtry (SNOTEL) station identifier\n+ W = WBAN identification number (last five characters of the \n+ GHCN-Daily ID)\n+\n+LATITUDE is latitude of the station (in decimal degrees).\n+\n+LONGITUDE is the longitude of the station (in decimal degrees).\n+\n+ELEVATION is the elevation of the station (in meters, missing = -999.9).\n+\n+\n+STATE is the U.S. postal code for the state (for U.S. stations only).\n+\n+NAME is the name of the station.\n+\n+GSN FLAG is a flag that indicates whether the station is part of the GCOS\n+ Surface Network (GSN). The flag is assigned by cross-referencing \n+ the number in the WMOID field with the official list of GSN \n+ stations. There are two possible values:\n+\n+ Blank = non-GSN station or WMO Station number not available\n+ GSN = GSN station \n+\n+HCN/ is a flag that indicates whether the station is part of the U.S.\n+CRN FLAG Historical Climatology Network (HCN). There are three possible \n+ values:\n+\n+ Blank = Not a member of the U.S. Historical Climatology \n+\t or U.S. Climate Reference Networks\n+ HCN = U.S. Historical Climatology Network station\n+\t CRN = U.S. Climate Reference Network or U.S. Regional Climate \n+\t Network Station\n+\n+WMO ID is the World Meteorological Organization (WMO) number for the\n+ station. If the station has no WMO number (or one has not yet \n+\t been matched to this station), then the field is blank.\n+\n+--------------------------------------------------------------------------------\n+--------------------------------------------------------------------------------\n+\n+ACW00011604 17.1167 -61.7833 10.1 ST JOHNS COOLIDGE FLD '..b'W00041606 19.2833 166.6500 4.3 UM WAKE ISLAND 91246\n+WZ004094600 -27.1700 31.2700 983.0 DWALENI \n+WZ004451000 -26.6700 31.0700 1030.0 MANKAYANE \n+WZ004455110 -26.5330 31.3000 641.0 MANZINI/MATSAPA AIR 68396\n+WZ004467410 -26.8500 31.9200 100.0 BIG BEND(WISSEL RODE) \n+WZ004822290 -26.3200 31.1300 1219.0 MBABANE \n+WZ004834260 -26.1000 31.7500 250.0 SWAZILAND RANCH-HOMEST \n+ZA000067403 -9.8000 29.0830 1324.0 KAWAMBWA 67403\n+ZA000067441 -11.7500 24.4330 1363.0 MWINILUNGA GSN 67441\n+ZA000067475 -10.2170 31.1330 1384.0 KASAMA GSN 67475\n+ZA000067531 -13.5330 23.1170 1078.0 ZAMBEZI 67531\n+ZA000067541 -13.5330 25.8500 1234.0 KASEMPA 67541\n+ZA000067543 -13.6000 24.2000 1075.0 KABOMPO 67543\n+ZA000067551 -12.1830 26.3830 1386.0 SOLWEZI 67551\n+ZA000067561 -13.0000 28.6500 1270.0 NDOLA 67561\n+ZA000067571 -13.2330 30.2170 1384.0 SERENJE 67571\n+ZA000067581 -13.5500 32.5830 1032.0 CHIPATA GSN 67581\n+ZA000067583 -12.2830 33.2000 1143.0 LUNDAZI 67583\n+ZA000067633 -15.2500 23.1500 1053.0 MONGU GSN 67633\n+ZA000067641 -14.8000 24.8000 1213.0 KAOMA 67641\n+ZA000067659 -15.7670 27.9170 978.0 KAFUE POLDER 67659\n+ZA000067667 -15.5500 28.2500 1213.0 MOUNT MAKULU 67667\n+ZA000067741 -17.4670 24.3000 951.0 SESHEKE 67741\n+ZA000067743 -17.8170 25.8170 986.0 LIVINGSTONE GSN 67743\n+ZA000067753 -16.8330 27.0670 1278.0 CHOMA 67753\n+ZAM00067663 -14.4500 28.4670 1207.0 KABWE/MILLIKEN 67663\n+ZI000067755 -17.6170 27.3330 617.0 BINGA 67755\n+ZI000067761 -16.5170 28.8830 518.0 KARIBA 67761\n+ZI000067765 -16.8330 29.6170 1344.0 KAROI 67765\n+ZI000067775 -17.9170 31.1330 1480.0 HARARE (KUTSAGA) GSN 67775\n+ZI000067779 -16.7830 31.5830 966.0 MOUNT DARWIN 67779\n+ZI000067781 -17.4170 32.2170 1244.0 MUTOKO 67781\n+ZI000067789 -17.0330 30.8500 1480.0 MVURWI 67789\n+ZI000067843 -18.1000 25.8500 1062.0 VICTORIA FALLS 67843\n+ZI000067853 -18.6330 27.0000 1077.0 HWANGE NATIONAL PAR 67853\n+ZI000067861 -18.2170 28.9330 1282.0 GOKWE 67861\n+ZI000067865 -18.9330 29.8330 1215.0 KWEKWE 67865\n+ZI000067867 -19.4500 29.8500 1429.0 GWERU 67867\n+ZI000067889 -18.2830 32.7500 1880.0 WYANGA 67889\n+ZI000067964 -20.1500 28.6170 1344.0 BULAWAYO (GOETZ OBS GSN 67964\n+ZI000067965 -20.0170 28.6170 1326.0 BULAWAYO AIRPORT 67965\n+ZI000067969 -21.0500 29.3670 861.0 WEST NICHOLSON 67969\n+ZI000067975 -20.0670 30.8670 1095.0 MASVINGO 67975\n+ZI000067977 -21.0170 31.5830 430.0 BUFFALO RANGE 67977\n+ZI000067983 -20.2000 32.6160 1132.0 CHIPINGE GSN 67983\n+ZI000067991 -22.2170 30.0000 457.0 BEITBRIDGE 67991\n' |
| b |
| diff -r e82ea5f1a179 -r 80d38a26b2e8 tool-data/locations.txt.sample --- a/tool-data/locations.txt.sample Mon Aug 14 13:44:45 2017 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -1,16 +0,0 @@ -Asheville NC asheville_nc:35.58 -Bridgeton NJ bridgeton_nj:39.43 -Davis CA davis_ca:38.55 -Geneva NY geneva_ny:42.88 -Homestead FL homestead_fl:25.47 -NJ_companion nj_companion:39.62 -NJ_management nj_management:39.66 -PA_companion pa_companion:39.98 -PA_management pa_management:39.97 -Riverside CA riverside_ca:33.95 -Salem OR salem_or:44.93 -VA_companion va_companion:39.15 -VA_management va_management:39.13 -Wenatchee WA wneatchee_wa:47.42 -WV_companion wv_companion:39.38 -WV_management wv_management:39.52 |
| b |
| diff -r e82ea5f1a179 -r 80d38a26b2e8 tool-data/years.txt.sample --- a/tool-data/years.txt.sample Mon Aug 14 13:44:45 2017 -0400 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 |
| b |
| @@ -1,31 +0,0 @@ -1995 1995 -1996 1996 -1997 1997 -1998 1998 -1999 1999 -2000 2000 -2001 2001 -2002 2002 -2003 2003 -2004 2004 -2005 2005 -2006 2006 -2007 2007 -2008 2008 -2009 2009 -2010 2010 -2011 2011 -2012 2012 -2013 2013 -2014 2014 -2015 2015 -2016 2016 -2017 2017 -2018 2018 -2019 2019 -2020 2020 -2021 2021 -2022 2022 -2023 2023 -2024 2024 -2025 2025 |