Skip to content

Commit 1a16b2e

Browse files
author
Jordan S Read
committed
Merge pull request #212 from lawinslow/master
Switch over to Sparkling example
2 parents 1523d9c + 1e208a8 commit 1a16b2e

21 files changed

+40280
-23333
lines changed

‎DESCRIPTION‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
Package: glmtools
22
Type: Package
33
Title: glmtools
4-
Version: 0.14.2
4+
Version: 0.14.3
55
Date: 2016-05-22
66
Authors@R: c( person("Jordan", "Read", role = c("aut","cre"),
77
email = "[email protected]"),

‎R/compare_to_field.R‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@
2727
#' metric = 'water.temperature', as_value = FALSE)
2828
#'print(paste(temp_rmse,'deg C RMSE'))
2929
#'# function in development
30-
#'buoy_file <- file.path(sim_folder, 'buoy_data.csv')
30+
#'buoy_file <- file.path(sim_folder, 'field_data.tsv')
3131
#'temp_rmse <- compare_to_field(nc_file, buoy_file,
3232
#' metric = 'water.temperature', as_value = FALSE,
3333
#' method = 'interp',precision = 'hours')

‎R/get_temp.R‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,7 @@
2424
#'temp_bot <- get_temp(nc_file)
2525
#'
2626
#'#-- get temporal subset--
27-
#'t_out <- seq(as.POSIXct("2011-04-04"), as.POSIXct("2011-06-01"), by = 86400)
27+
#'t_out <- seq(as.POSIXct("2010-04-15"), as.POSIXct("2010-06-01"), by = 86400)
2828
#'temp_surf <- get_temp(nc_file, reference = 'surface', z_out = 0, t_out = t_out)
2929
#'plot(temp_surf)
3030
#'@export

‎R/resample_to_field.R‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
#'
1919
#'temps <- resample_to_field(nc_file, field_file)
2020
#'buoy_file <- system.file('extdata', 'buoy_data.csv', package = 'glmtools')
21-
#'temps <- resample_to_field(nc_file, buoy_file, precision = 'hours')
21+
#'temps <- resample_to_field(nc_file, buoy_file, precision = 'mins')
2222
#'@export
2323
resample_to_field <- function(nc_file, field_file, method = 'match', precision = 'days', var_name='temp'){
2424

‎R/run_example_sim.R‎

Lines changed: 27 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -20,44 +20,56 @@ run_example_sim = function(sim_folder, verbose = TRUE){
2020
if(verbose){cat('sim_folder argument missing, using temp directory: ', sim_folder,'\n\n')}
2121
}
2222
nml_file <- file.path(sim_folder, 'glm2.nml')
23-
driver_file <- file.path(sim_folder, 'Anvil_driver.csv')
23+
driver_file <- file.path(sim_folder, 'nldas_driver.csv')
2424
calibration_tsv <- file.path(sim_folder, 'field_data.tsv')
2525
calibration_csv <- file.path(sim_folder, 'field_data.csv')
2626
stage_csv <- file.path(sim_folder, 'field_stage.csv')
27+
ice_snow_csv <- file.path(sim_folder, 'snow_ice_depth_obs.csv')
28+
ice_dur_csv <- file.path(sim_folder, 'ice_duration_obs.csv')
2729

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

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

3636
file.copy(from = system.file('extdata', 'field_data.tsv', package = 'glmtools'), to = calibration_tsv)
3737
file.copy(from = system.file('extdata', 'field_data.csv', package = 'glmtools'), to = calibration_csv)
38-
file.copy(from = system.file('extdata', 'buoy_data.csv', package = 'glmtools'), to = buoy_csv)
3938
file.copy(from = system.file('extdata', 'field_stage.csv', package = 'glmtools'), to = stage_csv)
39+
file.copy(from = system.file('extdata', basename(ice_snow_csv), package = 'glmtools'), to = ice_snow_csv)
40+
file.copy(from = system.file('extdata', basename(ice_dur_csv), package = 'glmtools'), to = ice_dur_csv)
4041

4142
nml <- read_nml() # read in default nml from GLMr
42-
nml <- set_nml(nml, arg_list = list('Kw'=0.55, 'lake_name'='Anvil',
43+
nml <- set_nml(nml, arg_list = list('Kw'=0.331, 'lake_name'='Sparkling',
4344
'bsn_vals' = 15,
44-
'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),
45-
'A' = c(0, 108964, 217929, 326893, 435858, 544822.5, 653787, 762751, 871716, 980680, 1089645, 1198609, 1307574, 1416538, 1525503),
46-
'start' = '2011-04-01 00:00:00',
47-
'stop' = '2011-09-02 00:00:00',
45+
'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),
46+
'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),
47+
'start' = '2010-04-15 00:00:00',
48+
'stop' = '2010-12-30 00:00:00',
4849
'dt' = 3600,
4950
'out_fn' = nc_out,
5051
'nsave' = 24,
5152
'csv_point_nlevs' = 0,
5253
'num_depths' = 3,
53-
'lake_depth' = 9.7536,
54-
'the_depths' = c(0, 1.2, 9.7536),
55-
'the_temps' = c(12, 10, 7),
54+
'lake_depth' = 18.288,
55+
'the_depths' = c(0, 0.2, 18.288),
56+
'the_temps' = c(3, 4, 4),
5657
'the_sals' = c(0, 0, 0),
5758
'subdaily' = FALSE,
58-
'meteo_fl' = 'Anvil_driver.csv',
59-
'cd' = 0.00108,
59+
'meteo_fl' = 'nldas_driver.csv',
60+
'max_layer_thick' = 3,
61+
'sw_factor' = 1.08,
62+
'coef_wind_stir' = 0.402,
63+
'coef_mix_hyp' = 0.2,
64+
'coef_mix_KH' = 0.1,
65+
'cd' = 0.0013,
66+
'ce' = 0.00132,
67+
'sed_temp_mean' = 4.5,
68+
'sed_temp_amplitude' = 0.25,
69+
'min_layer_thick' = 0.1,
70+
'coef_mix_conv' = 0.20,
6071
'num_outlet' = 0))
72+
6173
if(verbose){cat('writing nml file to ', nml_file,'\n')}
6274
write_nml(glm_nml = nml, file = nml_file)
6375

‎R/time_helpers.R‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -87,9 +87,9 @@ get_prec_time <- function(time_secs){
8787
if (time_secs >= 3600 & time_secs < 86400){
8888
prec = 'hours'
8989
} else if (time_secs >= 60 & time_secs < 3600){
90-
prec = 'minutes'
90+
prec = 'mins'
9191
} else if (time_secs < 60){
92-
prec = 'seconds'
92+
prec = 'secs'
9393
} else {
9494
prec = 'days'
9595
}
@@ -100,9 +100,9 @@ get_sec_unit <- function(unit){
100100
# gotta be a POSIXct method for this...
101101
if (unit == 'hours'){
102102
secs = 3600
103-
} else if (unit == 'minutes'){
103+
} else if (unit == 'mins'){
104104
secs = 60
105-
} else if (unit == 'seconds'){
105+
} else if (unit == 'secs'){
106106
secs = 1
107107
} else if (unit == 'days'){
108108
secs = 86400

‎demo/test_glm_ice_sparkling.R‎

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
1+
2+
3+
library(glmtools)
4+
sim_dir = run_example_sim()
5+
6+
Sys.setenv(tz='UTC')
7+
8+
#compare temps
9+
ncfile = file.path(sim_dir, 'output.nc')
10+
tempfield = file.path(sim_dir, 'field_data.csv')
11+
plot_temp_compare(ncfile, tempfield, resample = FALSE)
12+
13+
rmse = validate_sim(ncfile, tempfield, 'water.temperature', fig_path=FALSE, report = TRUE)
14+
15+
################################################################################
16+
#Ice compare
17+
################################################################################
18+
ice_obs = read.table(file.path(sim_dir, 'snow_ice_depth_obs.csv'), sep=',', header=TRUE)
19+
ice_obs$DateTime = as.POSIXct(ice_obs$sampledate, tz="GMT")
20+
ice_obs$avsnow = ice_obs$avsnow/100
21+
ice_obs$whiteice = ice_obs$whiteice/100
22+
ice_obs$blueice = ice_obs$blueice/100
23+
24+
25+
snow = get_var(ncfile, 'hsnow')
26+
attributes(snow$DateTime)$tzone = 'GMT'
27+
wice = get_var(ncfile, 'hwice')
28+
attributes(wice$DateTime)$tzone = 'GMT'
29+
bice = get_var(ncfile, 'hice')
30+
attributes(bice$DateTime)$tzone = 'GMT'
31+
wtr = get_temp(ncfile, reference='surface')
32+
attributes(wtr$DateTime)$tzone = 'GMT'
33+
34+
all_ice = merge(snow, ice_obs, all.x=TRUE, by='DateTime')
35+
all_ice = merge(wice, all_ice, all.x=TRUE)
36+
all_ice = merge(bice, all_ice, all.x=TRUE)
37+
all_ice$yday = as.POSIXlt(all_ice$DateTime)$yday
38+
all_ice$year = as.POSIXlt(all_ice$DateTime)$year + 1900
39+
40+
ice = get_ice(ncfile)
41+
attributes(ice$DateTime)$tzone = 'GMT'
42+
ice_on_off = mda.lakes::get_ice_onoff(ice, wtr)
43+
44+
lter_on_off = read.table(file.path(sim_dir,'ice_duration_obs.csv'), sep=',', header=TRUE)
45+
all_ice = merge(all_ice, ice_on_off, all.x=TRUE)
46+
all_ice = merge(all_ice, lter_on_off, all.x=TRUE)
47+
48+
all_ice = all_ice[order(all_ice$DateTime), ]
49+
50+
################################################################################
51+
## Now plots
52+
################################################################################
53+
plot(all_ice$avsnow, all_ice$hsnow)
54+
abline(0,1)
55+
56+
plot(all_ice$whiteice, all_ice$hwice)
57+
abline(0,1)
58+
59+
plot(all_ice$blueice, all_ice$hice)
60+
abline(0,1)
61+
62+
plot(hice~DateTime, subset(all_ice, year < 1990), type='l')
63+
points(all_ice$DateTime, all_ice$blueice)
64+
65+
#png('figures/sp_bice.png', res=300, width=2000, height=2000)
66+
par(mfrow=c(3,1))
67+
plot(hice~DateTime, subset(all_ice, year < 1990), type='l')
68+
title('Blue Ice')
69+
points(all_ice$DateTime, all_ice$blueice)
70+
plot(hice~DateTime, subset(all_ice, year > 1990 & year < 2000), type='l')
71+
points(all_ice$DateTime, all_ice$blueice)
72+
plot(hice~DateTime, subset(all_ice, year > 2000 & year < 2010), type='l')
73+
points(all_ice$DateTime, all_ice$blueice)
74+
#dev.off()
75+
76+
#png('figures/sp_wice.png', res=300, width=2000, height=2000)
77+
par(mfrow=c(3,1))
78+
plot(hwice~DateTime, subset(all_ice, year < 1990), type='l')
79+
title('White Ice')
80+
points(all_ice$DateTime, all_ice$whiteice)
81+
plot(hwice~DateTime, subset(all_ice, year > 1990 & year < 2000), type='l')
82+
points(all_ice$DateTime, all_ice$whiteice)
83+
plot(hwice~DateTime, subset(all_ice, year > 2000 & year < 2010), type='l')
84+
points(all_ice$DateTime, all_ice$whiteice)
85+
#dev.off()
86+
87+
#png('figures/sp_snow.png', res=300, width=2000, height=2000)
88+
par(mfrow=c(3,1))
89+
plot(hsnow~DateTime, subset(all_ice, year < 1990), type='l')
90+
title('Snow')
91+
points(all_ice$DateTime, all_ice$avsnow)
92+
plot(hsnow~DateTime, subset(all_ice, year > 1990 & year < 2000), type='l')
93+
points(all_ice$DateTime, all_ice$avsnow)
94+
plot(hsnow~DateTime, subset(all_ice, year > 2000 & year < 2010), type='l')
95+
points(all_ice$DateTime, all_ice$avsnow)
96+
#dev.off()
97+
98+
99+
100+
####
101+
#ICE
102+
103+
mean(abs(all_ice$firstopen - as.POSIXlt(all_ice$off)$yday+1), na.rm=TRUE)
104+
mean((all_ice$firstopen - as.POSIXlt(all_ice$off)$yday+1), na.rm=TRUE)
105+
106+
plot(all_ice$firstopen, as.POSIXlt(all_ice$off)$yday+1)
107+
abline(0,1)
108+
109+
mean(abs(all_ice$lastopen - as.POSIXlt(all_ice$on)$yday+1), na.rm=TRUE)
110+
mean((all_ice$lastopen - as.POSIXlt(all_ice$on)$yday+1), na.rm=TRUE)
111+
112+
plot(all_ice$lastopen, as.POSIXlt(all_ice$on)$yday+1)
113+
abline(0,1)
114+
115+
116+
##figure for Hipsey paper
117+
118+
#png('figures/sp_bws_plots.png', res=300, width=2000, height=2000)
119+
par(mfrow=c(3,1), oma=c(5,4,0,0), mar=c(0,0,1,1))
120+
121+
plot(hice~DateTime, subset(all_ice, year > 1990 & year < 2001), type='l', ylab='Blue Ice (m)', xaxt='n')
122+
lines(hice~DateTime, all_ice)
123+
points(all_ice$DateTime, all_ice$blueice, col=rgb(0, 0, 0, 0.5), pch=16)
124+
mtext(side=2, text='Blue Ice (m)', line=2.5)
125+
legend('topright', legend = 'A', bty = 'n', inset = c(-0.01,-0.06), cex=2)
126+
127+
plot(hwice~DateTime, subset(all_ice, year > 1990 & year < 2001), type='l', ylab='White Ice (m)', xaxt='n')
128+
lines(hwice~DateTime, all_ice)
129+
points(all_ice$DateTime, all_ice$whiteice, col=rgb(0, 0, 0, 0.5), pch=16)
130+
mtext(side=2, text='White Ice (m)', line=2.5)
131+
legend('topright', legend = 'B', bty = 'n', inset = c(-0.01,-0.06), cex=2)
132+
133+
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))
134+
135+
plot(hsnow~DateTime, subset(all_ice, year > 1990 & year < 2001), type='l', ylab='Snow (m)')
136+
lines(hsnow~DateTime, all_ice)
137+
points(all_ice$DateTime, all_ice$avsnow, col=rgb(0, 0, 0, 0.5), pch=16)
138+
mtext(side=2, text='Snow (m)', line=2.5)
139+
mtext(side=1, text='Year', line=2.5)
140+
legend('topright', legend = 'C', bty = 'n', inset = c(-0.01,-0.06), cex=2)
141+
142+
#dev.off()
143+
144+
145+
### Compare ice on and off for Hipsey paper
146+
147+
library(mda.lakes)
148+
library(lubridate)
149+
onoff = lter_on_off
150+
mod_onoff = get_ice_onoff(get_ice(ncfile), get_temp(ncfile, reference = 'surface'))
151+
152+
obsmod = merge(onoff, mod_onoff, by='year')
153+
154+
155+
#png('figures/sp_on_off_plots.png', res=300, width=2400, height=1500)
156+
par(mfrow=c(1,2), mar=c(5,4,1,1))
157+
158+
plot(yday(obsmod$on), yday(as.POSIXct(obsmod$datelastopen)), pch=16, col=rgb(0, 0, 0, 0.5),
159+
xlab='Modeled Ice On', ylab='Observed Ice On')
160+
abline(0,1)
161+
plot(yday(obsmod$off[-1]), yday(as.POSIXct(obsmod$datefirstopen[-1])), pch=16, col=rgb(0, 0, 0, 0.5),
162+
xlab='Modeled Ice Off', ylab='Observed Ice Off')
163+
abline(0,1)
164+
165+
#dev.off()

0 commit comments

Comments
 (0)