Skip to content

Commit

Permalink
Merge pull request #212 from lawinslow/master
Browse files Browse the repository at this point in the history
Switch over to Sparkling example
  • Loading branch information
Jordan S Read committed May 26, 2016
2 parents 1523d9c + 1e208a8 commit 1a16b2e
Show file tree
Hide file tree
Showing 21 changed files with 40,280 additions and 23,333 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: glmtools
Type: Package
Title: glmtools
Version: 0.14.2
Version: 0.14.3
Date: 2016-05-22
Authors@R: c( person("Jordan", "Read", role = c("aut","cre"),
email = "jread@usgs.gov"),
Expand Down
2 changes: 1 addition & 1 deletion R/compare_to_field.R
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,7 @@
#' metric = 'water.temperature', as_value = FALSE)
#'print(paste(temp_rmse,'deg C RMSE'))
#'# function in development
#'buoy_file <- file.path(sim_folder, 'buoy_data.csv')
#'buoy_file <- file.path(sim_folder, 'field_data.tsv')
#'temp_rmse <- compare_to_field(nc_file, buoy_file,
#' metric = 'water.temperature', as_value = FALSE,
#' method = 'interp',precision = 'hours')
Expand Down
2 changes: 1 addition & 1 deletion R/get_temp.R
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@
#'temp_bot <- get_temp(nc_file)
#'
#'#-- get temporal subset--
#'t_out <- seq(as.POSIXct("2011-04-04"), as.POSIXct("2011-06-01"), by = 86400)
#'t_out <- seq(as.POSIXct("2010-04-15"), as.POSIXct("2010-06-01"), by = 86400)
#'temp_surf <- get_temp(nc_file, reference = 'surface', z_out = 0, t_out = t_out)
#'plot(temp_surf)
#'@export
Expand Down
2 changes: 1 addition & 1 deletion R/resample_to_field.R
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
#'
#'temps <- resample_to_field(nc_file, field_file)
#'buoy_file <- system.file('extdata', 'buoy_data.csv', package = 'glmtools')
#'temps <- resample_to_field(nc_file, buoy_file, precision = 'hours')
#'temps <- resample_to_field(nc_file, buoy_file, precision = 'mins')
#'@export
resample_to_field <- function(nc_file, field_file, method = 'match', precision = 'days', var_name='temp'){

Expand Down
42 changes: 27 additions & 15 deletions R/run_example_sim.R
Original file line number Diff line number Diff line change
Expand Up @@ -20,44 +20,56 @@ run_example_sim = function(sim_folder, verbose = TRUE){
if(verbose){cat('sim_folder argument missing, using temp directory: ', sim_folder,'\n\n')}
}
nml_file <- file.path(sim_folder, 'glm2.nml')
driver_file <- file.path(sim_folder, 'Anvil_driver.csv')
driver_file <- file.path(sim_folder, 'nldas_driver.csv')
calibration_tsv <- file.path(sim_folder, 'field_data.tsv')
calibration_csv <- file.path(sim_folder, 'field_data.csv')
stage_csv <- file.path(sim_folder, 'field_stage.csv')
ice_snow_csv <- file.path(sim_folder, 'snow_ice_depth_obs.csv')
ice_dur_csv <- file.path(sim_folder, 'ice_duration_obs.csv')

buoy_csv <- file.path(sim_folder, 'buoy_data.csv')
nc_file <- file.path(sim_folder, paste0(nc_out, '.nc'))
# move glm2.nml to sim_folder

file.copy(from = system.file('extdata', 'Anvil_driver.csv', package = 'glmtools'),
to = driver_file)
file.copy(from = system.file('extdata', 'nldas_driver.csv', package = 'glmtools'), to = driver_file)
if(verbose){cat('driver data file copied to ', driver_file,'\n')}

file.copy(from = system.file('extdata', 'field_data.tsv', package = 'glmtools'), to = calibration_tsv)
file.copy(from = system.file('extdata', 'field_data.csv', package = 'glmtools'), to = calibration_csv)
file.copy(from = system.file('extdata', 'buoy_data.csv', package = 'glmtools'), to = buoy_csv)
file.copy(from = system.file('extdata', 'field_stage.csv', package = 'glmtools'), to = stage_csv)
file.copy(from = system.file('extdata', basename(ice_snow_csv), package = 'glmtools'), to = ice_snow_csv)
file.copy(from = system.file('extdata', basename(ice_dur_csv), package = 'glmtools'), to = ice_dur_csv)

nml <- read_nml() # read in default nml from GLMr
nml <- set_nml(nml, arg_list = list('Kw'=0.55, 'lake_name'='Anvil',
nml <- set_nml(nml, arg_list = list('Kw'=0.331, 'lake_name'='Sparkling',
'bsn_vals' = 15,
'H' = c(510.5363, 511.23299, 511.92967, 512.62636, 513.32304, 514.01973, 514.71641, 515.4131, 516.10979, 516.80647, 517.50316, 518.19984, 518.89653, 519.59321, 520.2899),
'A' = c(0, 108964, 217929, 326893, 435858, 544822.5, 653787, 762751, 871716, 980680, 1089645, 1198609, 1307574, 1416538, 1525503),
'start' = '2011-04-01 00:00:00',
'stop' = '2011-09-02 00:00:00',
'H' = c(301.712, 303.018285714286, 304.324571428571, 305.630857142857, 306.937142857143, 308.243428571429, 309.549714285714, 310.856, 312.162285714286, 313.468571428571, 314.774857142857, 316.081142857143, 317.387428571429, 318.693714285714, 320),
'A' = c(0, 45545.8263571429, 91091.6527142857, 136637.479071429, 182183.305428571, 227729.131785714, 273274.958142857, 318820.7845, 364366.610857143, 409912.437214286, 455458.263571429, 501004.089928571, 546549.916285714, 592095.742642857, 637641.569),
'start' = '2010-04-15 00:00:00',
'stop' = '2010-12-30 00:00:00',
'dt' = 3600,
'out_fn' = nc_out,
'nsave' = 24,
'csv_point_nlevs' = 0,
'num_depths' = 3,
'lake_depth' = 9.7536,
'the_depths' = c(0, 1.2, 9.7536),
'the_temps' = c(12, 10, 7),
'lake_depth' = 18.288,
'the_depths' = c(0, 0.2, 18.288),
'the_temps' = c(3, 4, 4),
'the_sals' = c(0, 0, 0),
'subdaily' = FALSE,
'meteo_fl' = 'Anvil_driver.csv',
'cd' = 0.00108,
'meteo_fl' = 'nldas_driver.csv',
'max_layer_thick' = 3,
'sw_factor' = 1.08,
'coef_wind_stir' = 0.402,
'coef_mix_hyp' = 0.2,
'coef_mix_KH' = 0.1,
'cd' = 0.0013,
'ce' = 0.00132,
'sed_temp_mean' = 4.5,
'sed_temp_amplitude' = 0.25,
'min_layer_thick' = 0.1,
'coef_mix_conv' = 0.20,
'num_outlet' = 0))

if(verbose){cat('writing nml file to ', nml_file,'\n')}
write_nml(glm_nml = nml, file = nml_file)

Expand Down
8 changes: 4 additions & 4 deletions R/time_helpers.R
Original file line number Diff line number Diff line change
Expand Up @@ -87,9 +87,9 @@ get_prec_time <- function(time_secs){
if (time_secs >= 3600 & time_secs < 86400){
prec = 'hours'
} else if (time_secs >= 60 & time_secs < 3600){
prec = 'minutes'
prec = 'mins'
} else if (time_secs < 60){
prec = 'seconds'
prec = 'secs'
} else {
prec = 'days'
}
Expand All @@ -100,9 +100,9 @@ get_sec_unit <- function(unit){
# gotta be a POSIXct method for this...
if (unit == 'hours'){
secs = 3600
} else if (unit == 'minutes'){
} else if (unit == 'mins'){
secs = 60
} else if (unit == 'seconds'){
} else if (unit == 'secs'){
secs = 1
} else if (unit == 'days'){
secs = 86400
Expand Down
165 changes: 165 additions & 0 deletions demo/test_glm_ice_sparkling.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,165 @@


library(glmtools)
sim_dir = run_example_sim()

Sys.setenv(tz='UTC')

#compare temps
ncfile = file.path(sim_dir, 'output.nc')
tempfield = file.path(sim_dir, 'field_data.csv')
plot_temp_compare(ncfile, tempfield, resample = FALSE)

rmse = validate_sim(ncfile, tempfield, 'water.temperature', fig_path=FALSE, report = TRUE)

################################################################################
#Ice compare
################################################################################
ice_obs = read.table(file.path(sim_dir, 'snow_ice_depth_obs.csv'), sep=',', header=TRUE)
ice_obs$DateTime = as.POSIXct(ice_obs$sampledate, tz="GMT")
ice_obs$avsnow = ice_obs$avsnow/100
ice_obs$whiteice = ice_obs$whiteice/100
ice_obs$blueice = ice_obs$blueice/100


snow = get_var(ncfile, 'hsnow')
attributes(snow$DateTime)$tzone = 'GMT'
wice = get_var(ncfile, 'hwice')
attributes(wice$DateTime)$tzone = 'GMT'
bice = get_var(ncfile, 'hice')
attributes(bice$DateTime)$tzone = 'GMT'
wtr = get_temp(ncfile, reference='surface')
attributes(wtr$DateTime)$tzone = 'GMT'

all_ice = merge(snow, ice_obs, all.x=TRUE, by='DateTime')
all_ice = merge(wice, all_ice, all.x=TRUE)
all_ice = merge(bice, all_ice, all.x=TRUE)
all_ice$yday = as.POSIXlt(all_ice$DateTime)$yday
all_ice$year = as.POSIXlt(all_ice$DateTime)$year + 1900

ice = get_ice(ncfile)
attributes(ice$DateTime)$tzone = 'GMT'
ice_on_off = mda.lakes::get_ice_onoff(ice, wtr)

lter_on_off = read.table(file.path(sim_dir,'ice_duration_obs.csv'), sep=',', header=TRUE)
all_ice = merge(all_ice, ice_on_off, all.x=TRUE)
all_ice = merge(all_ice, lter_on_off, all.x=TRUE)

all_ice = all_ice[order(all_ice$DateTime), ]

################################################################################
## Now plots
################################################################################
plot(all_ice$avsnow, all_ice$hsnow)
abline(0,1)

plot(all_ice$whiteice, all_ice$hwice)
abline(0,1)

plot(all_ice$blueice, all_ice$hice)
abline(0,1)

plot(hice~DateTime, subset(all_ice, year < 1990), type='l')
points(all_ice$DateTime, all_ice$blueice)

#png('figures/sp_bice.png', res=300, width=2000, height=2000)
par(mfrow=c(3,1))
plot(hice~DateTime, subset(all_ice, year < 1990), type='l')
title('Blue Ice')
points(all_ice$DateTime, all_ice$blueice)
plot(hice~DateTime, subset(all_ice, year > 1990 & year < 2000), type='l')
points(all_ice$DateTime, all_ice$blueice)
plot(hice~DateTime, subset(all_ice, year > 2000 & year < 2010), type='l')
points(all_ice$DateTime, all_ice$blueice)
#dev.off()

#png('figures/sp_wice.png', res=300, width=2000, height=2000)
par(mfrow=c(3,1))
plot(hwice~DateTime, subset(all_ice, year < 1990), type='l')
title('White Ice')
points(all_ice$DateTime, all_ice$whiteice)
plot(hwice~DateTime, subset(all_ice, year > 1990 & year < 2000), type='l')
points(all_ice$DateTime, all_ice$whiteice)
plot(hwice~DateTime, subset(all_ice, year > 2000 & year < 2010), type='l')
points(all_ice$DateTime, all_ice$whiteice)
#dev.off()

#png('figures/sp_snow.png', res=300, width=2000, height=2000)
par(mfrow=c(3,1))
plot(hsnow~DateTime, subset(all_ice, year < 1990), type='l')
title('Snow')
points(all_ice$DateTime, all_ice$avsnow)
plot(hsnow~DateTime, subset(all_ice, year > 1990 & year < 2000), type='l')
points(all_ice$DateTime, all_ice$avsnow)
plot(hsnow~DateTime, subset(all_ice, year > 2000 & year < 2010), type='l')
points(all_ice$DateTime, all_ice$avsnow)
#dev.off()



####
#ICE

mean(abs(all_ice$firstopen - as.POSIXlt(all_ice$off)$yday+1), na.rm=TRUE)
mean((all_ice$firstopen - as.POSIXlt(all_ice$off)$yday+1), na.rm=TRUE)

plot(all_ice$firstopen, as.POSIXlt(all_ice$off)$yday+1)
abline(0,1)

mean(abs(all_ice$lastopen - as.POSIXlt(all_ice$on)$yday+1), na.rm=TRUE)
mean((all_ice$lastopen - as.POSIXlt(all_ice$on)$yday+1), na.rm=TRUE)

plot(all_ice$lastopen, as.POSIXlt(all_ice$on)$yday+1)
abline(0,1)


##figure for Hipsey paper

#png('figures/sp_bws_plots.png', res=300, width=2000, height=2000)
par(mfrow=c(3,1), oma=c(5,4,0,0), mar=c(0,0,1,1))

plot(hice~DateTime, subset(all_ice, year > 1990 & year < 2001), type='l', ylab='Blue Ice (m)', xaxt='n')
lines(hice~DateTime, all_ice)
points(all_ice$DateTime, all_ice$blueice, col=rgb(0, 0, 0, 0.5), pch=16)
mtext(side=2, text='Blue Ice (m)', line=2.5)
legend('topright', legend = 'A', bty = 'n', inset = c(-0.01,-0.06), cex=2)

plot(hwice~DateTime, subset(all_ice, year > 1990 & year < 2001), type='l', ylab='White Ice (m)', xaxt='n')
lines(hwice~DateTime, all_ice)
points(all_ice$DateTime, all_ice$whiteice, col=rgb(0, 0, 0, 0.5), pch=16)
mtext(side=2, text='White Ice (m)', line=2.5)
legend('topright', legend = 'B', bty = 'n', inset = c(-0.01,-0.06), cex=2)

legend('topleft', legend = c('Modeled', 'Observed'), pch = c(NA, 16), lty=c(1,NA), col=c('black', rgb(0, 0, 0, 0.5)), inset=c(0.01, 0.05))

plot(hsnow~DateTime, subset(all_ice, year > 1990 & year < 2001), type='l', ylab='Snow (m)')
lines(hsnow~DateTime, all_ice)
points(all_ice$DateTime, all_ice$avsnow, col=rgb(0, 0, 0, 0.5), pch=16)
mtext(side=2, text='Snow (m)', line=2.5)
mtext(side=1, text='Year', line=2.5)
legend('topright', legend = 'C', bty = 'n', inset = c(-0.01,-0.06), cex=2)

#dev.off()


### Compare ice on and off for Hipsey paper

library(mda.lakes)
library(lubridate)
onoff = lter_on_off
mod_onoff = get_ice_onoff(get_ice(ncfile), get_temp(ncfile, reference = 'surface'))

obsmod = merge(onoff, mod_onoff, by='year')


#png('figures/sp_on_off_plots.png', res=300, width=2400, height=1500)
par(mfrow=c(1,2), mar=c(5,4,1,1))

plot(yday(obsmod$on), yday(as.POSIXct(obsmod$datelastopen)), pch=16, col=rgb(0, 0, 0, 0.5),
xlab='Modeled Ice On', ylab='Observed Ice On')
abline(0,1)
plot(yday(obsmod$off[-1]), yday(as.POSIXct(obsmod$datefirstopen[-1])), pch=16, col=rgb(0, 0, 0, 0.5),
xlab='Modeled Ice Off', ylab='Observed Ice Off')
abline(0,1)

#dev.off()
Loading

0 comments on commit 1a16b2e

Please sign in to comment.