Skip to content

Commit

Permalink
putting plot of most recent data against forecast on website (#111)
Browse files Browse the repository at this point in the history
* changed name of plot_data function to plot_species_forecast for clarity

* added code to plot the most recent data from the repo against the most recent forecast prior to that data collection.  there are still some issues, right now this will only work when June 2017 is the most recent data

* removed old code from report.Rmd
  • Loading branch information
emchristensen authored and sdtaylor committed Sep 27, 2017
1 parent 2fbb0f5 commit 7fc9d13
Show file tree
Hide file tree
Showing 3 changed files with 39 additions and 32 deletions.
23 changes: 10 additions & 13 deletions forecast_tools.R
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,13 @@ get_sp_predicts = function(data, lvl, lead_time) {
data2 = filter(data1, newmoonnumber == target_moon)
}

plot_data = function(data) {
#' this is the 'second plot on the 'Species-level Forecast' on the 'Current Forecast' page on the website
#'
#'
#' @param data
#' @return sp_predict is a plot object -- plot(sp_predict) displays it
#'
plot_species_forecast = function(data) {
newmoons_table = read.csv(
text = getURL(
"https://raw.githubusercontent.com/weecology/PortalData/master/Rodents/moon_dates.csv"))
Expand Down Expand Up @@ -111,21 +117,12 @@ plot_data = function(data) {
ylab("Species") +
xlab("Abundance") +
scale_y_discrete(breaks = reorder(data$species,data$estimate),labels = reorder(species_names$scientificname,species_names$estimate))
if(!is.na(period_code)){
rodents = abundance("repo", shape='flat')
observed = dplyr::filter(rodents, period == period_code)
joined_data = left_join(data, observed, by = "species")
joined_data[is.na(joined_data)] = 0
joined_data[joined_data$species=='total','abundance'] = sum(joined_data$abundance,na.rm=T)
joined_data[joined_data$species=='total','period'] = period_code
sp_predict = sp_predict +
geom_point(data = joined_data, mapping = aes(x = abundance, y = species),
color = "blue")
}
plot(sp_predict)

return(sp_predict)
}



#' Compares forecasts to observations over different lead times.
#' Error can be any function. The level, species, and currency columns from
#' observations and forecasts must have matching values.
Expand Down
3 changes: 2 additions & 1 deletion index.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,8 @@ source("forecast_tools.R")
data = compile_forecasts()
ensemble = dplyr::filter(data, level == 'All', model == 'Ensemble', date == max(date))
sp_predictions = get_sp_predicts(ensemble, 'All', lead_time = 1)
plot_data(sp_predictions)
sp_predict = plot_species_forecast(sp_predictions)
plot(sp_predict)
```
```{r highabund_ts_forecasts, echo=FALSE, message=FALSE, warning=FALSE}
most_abund_sp = sp_predictions %>%
Expand Down
45 changes: 27 additions & 18 deletions report.Rmd
Original file line number Diff line number Diff line change
@@ -1,27 +1,36 @@
---
title: "Report - June 2017"
title: "Report - September 2017"
---
Here's the most recent data (blue) compared to the forecast performed right before data collection (black, with error bars)
```{r, echo=FALSE, message=FALSE, warning=FALSE}
source("forecast_tools.R")
library(portalr)
Note the forecast for D. spectabalis this month.
# get most recent rodent data
rodents = abundance("repo", shape='flat')
period_code = max(rodents$period)
observed = dplyr::filter(rodents, period == period_code)
```{r speciesforecasts, echo=FALSE, message=FALSE, warning=FALSE}
source("forecast_tools.R")
# get new moon number associated with this period
newmoons_table = read.csv(text = getURL("https://raw.githubusercontent.com/weecology/PortalData/master/Rodents/moon_dates.csv"))
moon_code = dplyr::filter(newmoons_table,period==period_code) %>% select(newmoonnumber) %>% as.integer()
# get most recent forecast prior to data collection -- WIP!!!!!!!! there's an issue up
data = read.csv('predictions/2017-06-23Allforecasts.csv',na.strings = '')
data = plyr::rename(data,c('NewMoonNumber' = 'newmoonnumber'))
data = compile_forecasts()
ensemble = dplyr::filter(data, level == 'All', model == 'Ensemble', date == max(date))
ensemble = dplyr::filter(data, level == 'All', model == 'Ensemble', newmoonnumber==moon_code)
sp_predictions = get_sp_predicts(ensemble, 'All', lead_time = 1)
ds_forecast = forecast_viz(obs_data = obs_data_newmoon,
obs_date_col_name = "newmoonnumber",
obs_val_col_name = 'DS',
for_data = for_data,
for_date_col_name = "newmoonnumber",
for_val_col_name = "estimate",
for_model_name = "ensemble",
for_lowerpi_col_name = "LowerPI",
for_upperpi_col_name = "UpperPI",
start_newmoon = 300)
plot(ds_forecast)
joined_data = left_join(sp_predictions, observed, by = "species")
joined_data[is.na(joined_data)] = 0
joined_data[joined_data$species=='total','abundance'] = sum(joined_data$abundance,na.rm=T)
joined_data[joined_data$species=='total','period'] = period_code
sp_predict = plot_species_forecast(sp_predictions)
sp_predict = sp_predict + geom_point(data = joined_data, mapping = aes(x = abundance, y = species),
color = "blue")
plot(sp_predict)
```

There is an increasing change that we will see a spectab in the coming months. Does this mean Stephanie will return?

0 comments on commit 7fc9d13

Please sign in to comment.