Formalize a small example

[mdodrill-usgs]

Formalize a small example which illustrates the full pipeline to run a streamMetabolizer model. This could then be added to the package for users to run/diagnose install problems.

Bob may have an example of this which we can build off.

[robohall]

Here is minimal code and data.

library(tidyverse)
library(lubridate)
library(streamMetabolizer)

#Here is minimal code for estimating metabolism from using streamMetabolizer and data from a PME sensor.
#Use it to test if sM is working on your computer and to help with building your own code
#Bob Hall
#update, Aug 2026


setwd("/Users/bob.hall/streampulse")
#set your own wd

##constants you will need.  Enter your own
lat<- 48.5
long<- -113.6 #(note negative for us folks in western hemisphere)
altitude<- 914 # m asl
bpstd <- 1013 #Do not change. ignoring weather relate changes in bp for this minimal code.  1013 mb = 760 mm Hg.  
K600_guess <- 10 #your best guess of K600 in units of 1/d.  Big slow rivers, 0.5 to 2.  Medium slow rivers, 2-10. 
#Fun rivers to float 10-20.  Non bubbly streams 20-40.  Bubbly streams 50-200.
mean_depth <- 0.25 # knowing mean river depth (m) is not important for code testing, but it sure is important
#if you want to publish meaningful estimates of metabolism.  Measure your streams's depth at many cross sections.

####function returns pressure in mb  which is what sM wants
bpcalc_atm<- function(bpst, alt) {
  bpst*exp((-9.80665*0.0289644*alt)/(8.31447*(273.15+15)))
}
bpcalc_atm(bpst=bpstd, alt=altitude)



###read in a .csv file that has unix time, oxygen conc and temp from a PME file
#You want 5 to 50 days of data here

##Make your .csv file look exactly like this.  Nothing else!  
# unixtime   temp   oxy
# 1721418480 18.048 8.467
# 1721419080 16.720 8.869
# 1721419680 16.490 8.965

#where unixtime is "Unix" , temp is  "Temperature_degC", and oxy is  "Dissolved Oxygen_mgL" in the PME file
#And god forbid never use PME % sat data


oxy<-read.csv("beaver_oxy.csv")
head(oxy)  #look ok??


####
##run the below code without messing with it

oxy$time<-as_datetime(oxy$unixtime) #converts unit time to a R time object in UTC

oxy$solar.time<-convert_UTC_to_solartime(oxy$time, longitude= long, time.type="mean solar")


oxy$light<- calc_light(oxy$solar.time, latitude=lat, longitude=long, max.PAR =2326, attach.units = F)

oxy$osat<- calc_DO_sat(temp.water=oxy$temp, pressure.air=bpcalc_atm(bpst=bpstd, alt=altitude))


#set sM model and specs.  Ignoring Q variation here with a simpler pooling model
#model is state space.

model_name <- mm_name(type='bayes', pool_K600='normal', err_obs_iid = T, err_proc_iid =T)
model_specs <- specs(model_name, K600_daily_meanlog_meanlog=log(K600_guess), K600_daily_meanlog_sdlog=0.7, K600_daily_sdlog_sigma=0.1, burnin_steps=1000, 
                  saved_steps=1000) 

### Make your data frame

data_sm<-data.frame(DO.obs=oxy$oxy, DO.sat=oxy$osat, 
                    temp.water=oxy$temp, depth=rep(mean_depth,length(oxy$oxy)), 
                    light=oxy$light, solar.time=oxy$solar.time)

head(data_sm) #always take a look at the dataframe!


fit <- metab(model_specs, data=data_sm, info=c(site='Beaver Brook', source='HT Odum'))


plot_DO_preds(predict_DO(fit))
plot_metab_preds(predict_metab(fit))
params<- get_params(fit , uncertainty='ci')

plot(params$K600.daily,params$ER.daily) #if these covary, then equifinality warning!  Indeed my data have that.


mean(params$GPP.daily, na.rm=T)
mean(params$ER.daily, na.rm=T)

[beaver_oxy.csv](https://github.com/user-attachments/files/31046584/beaver_oxy.csv)