1- # ' @title Plot matching heatmaps for modeled and observed temp
2- # ' @param nc_file Netcdf model output file
3- # ' @param field_file CSV or TSV field data file (see \link{resample_to_field} for format)
4- # ' @param var_name a character vector of valid variable names (see \code{\link{sim_vars}})
5- # ' @param fig_path F if plot to screen, string path if save plot as .png
6- # ' @param resample sample the model output to the same time points as the observations?
7- # ' @param \dots additional arguments passed to \code{\link{resample_to_field}}
8- # '
9- # ' @seealso Internally uses \link{get_var} and \link{resample_to_field}
10- # '
11- # '
12- # '@examples
13- # 'sim_folder <- run_example_sim(verbose = FALSE)
14- # 'nc_file <- file.path(sim_folder, 'output.nc')
15- # 'nml_file <- file.path(sim_folder, 'glm2.nml')
16- # 'field_file <- file.path(sim_folder, 'field_data.tsv')
17- # '
18- # 'run_glm(sim_folder)
19- # '
20- # 'plot_var_compare(nc_file, field_file, 'temp', resample=FALSE) ##makes a plot!
21- # '
22- # '@importFrom akima interp
23- # '@export
24- plot_var_compare = function (nc_file , field_file , var_name , fig_path = FALSE , resample = TRUE , ... ){
25-
26- heatmaps <- .is_heatmap(nc_file , var_name )
27- if (! heatmaps ){
28- warning(' plot_var_compare not implemented for 1D variables'
29- return ()
30- }
31-
32- start_par = par(no.readonly = TRUE )
33- # Create layout
34-
35- mod_temp = get_var(nc_file , var_name , reference = ' surface'
36- mod_depths = get.offsets(mod_temp )
37-
38-
39- data = resample_to_field(nc_file , field_file , var_name = var_name , ... )
40- if (resample ){
41- model_df <- resample_sim(mod_temp , t_out = unique(data $ DateTime ))
42- }else {
43- model_df = mod_temp
44- }
45-
46- # Pivot observed into table
47- x = as.numeric(as.POSIXct(data $ DateTime ))
48- y = data $ Depth
49- z = data [,paste0(' Observed_' var_name )]
50- x_out = sort(unique(x ))
51- y_out = sort(unique(mod_depths ))
52-
53- # remove any NA values before the 2D interp
54- x = x [! is.na(z )]
55- y = y [! is.na(z )]
56- z = z [! is.na(z )]
57-
58- # Added a scaling factor to Y. Interp won't interpolate if X and Y are on vastly different scales.
59- # I don't use Y from here later, so it doesn't matter what the mangitude of the values is.
60- interped = interp(x , y * 1e6 , z , x_out , y_out * 1e6 )
61-
62- gen_default_fig(filename = fig_path , num_divs = 2 )# , omi = c(0.1, 0.5, 0, 0))
63- .stacked_layout(heatmaps , num_divs = 2 )
64- obs_df <- data.frame (interped $ z )
65- names(obs_df ) <- paste(' var_' y_out , sep = ' '
66- obs_df <- cbind(data.frame (DateTime = as.POSIXct(x_out , origin = ' 1970-01-01' obs_df )
67-
68- # Use model to define X-axis plotting extent for both graphs
69- xaxis <- get_xaxis(model_df [,1 ])
70-
71- y.text = y_out [1 ]+ diff(range(y_out ))* 0.05 # note, reference will ALWAYS be surface for compare to field data
72- .plot_df_heatmap(obs_df , bar_title = .unit_label(nc_file ,var_name ), overlays = c(points(x = x ,y = y ),text(x_out [1 ],y = y.text ,' Observed' pos = 4 , offset = 1 )), xaxis = xaxis )
73-
74- .plot_df_heatmap(model_df , bar_title = .unit_label(nc_file ,var_name ), overlays = text(x_out [1 ],y = y.text ,' Modeled' pos = 4 , offset = 1 ), xaxis = xaxis )
75-
76- par(start_par )# set PAR back to what it started at
77- if (is.character(fig_path ))
78- dev.off()
1+ # ' @title Plot matching heatmaps for modeled and observed temp
2+ # ' @param nc_file Netcdf model output file
3+ # ' @param field_file CSV or TSV field data file (see \link{resample_to_field} for format)
4+ # ' @param var_name a character vector of valid variable names (see \code{\link{sim_vars}})
5+ # ' @param fig_path F if plot to screen, string path if save plot as .png
6+ # ' @param resample sample the model output to the same time points as the observations?
7+ # ' @param col_lim range for heatmap (in units of the variable)
8+ # ' @param \dots additional arguments passed to \code{\link{resample_to_field}}
9+ # '
10+ # ' @seealso Internally uses \link{get_var} and \link{resample_to_field}
11+ # '
12+ # '
13+ # '@examples
14+ # 'sim_folder <- run_example_sim(verbose = FALSE)
15+ # 'nc_file <- file.path(sim_folder, 'output.nc')
16+ # 'nml_file <- file.path(sim_folder, 'glm2.nml')
17+ # 'field_file <- file.path(sim_folder, 'field_data.tsv')
18+ # '
19+ # 'run_glm(sim_folder)
20+ # '
21+ # 'plot_var_compare(nc_file, field_file, 'temp', resample=FALSE) ##makes a plot!
22+ # '
23+ # '@importFrom akima interp
24+ # '@export
25+ plot_var_compare = function (nc_file , field_file , var_name , fig_path = FALSE , resample = TRUE , col_lim , ... ){
26+
27+ heatmaps <- .is_heatmap(nc_file , var_name )
28+ if (! heatmaps ){
29+ warning(' plot_var_compare not implemented for 1D variables'
30+ return ()
31+ }
32+
33+ start_par = par(no.readonly = TRUE )
34+ # Create layout
35+
36+ mod_temp = get_var(nc_file , var_name , reference = ' surface'
37+ mod_depths = get.offsets(mod_temp )
38+
39+
40+ data = resample_to_field(nc_file , field_file , var_name = var_name , ... )
41+ if (resample ){
42+ model_df <- resample_sim(mod_temp , t_out = unique(data $ DateTime ))
43+ }else {
44+ model_df = mod_temp
45+ }
46+
47+ # Pivot observed into table
48+ x = as.numeric(as.POSIXct(data $ DateTime ))
49+ y = data $ Depth
50+ z = data [,paste0(' Observed_' var_name )]
51+ x_out = sort(unique(x ))
52+ y_out = sort(unique(mod_depths ))
53+
54+ # remove any NA values before the 2D interp
55+ x = x [! is.na(z )]
56+ y = y [! is.na(z )]
57+ z = z [! is.na(z )]
58+
59+ # Added a scaling factor to Y. Interp won't interpolate if X and Y are on vastly different scales.
60+ # I don't use Y from here later, so it doesn't matter what the mangitude of the values is.
61+ interped = interp(x , y * 1e6 , z , x_out , y_out * 1e6 )
62+
63+ gen_default_fig(filename = fig_path , num_divs = 2 )# , omi = c(0.1, 0.5, 0, 0))
64+ .stacked_layout(heatmaps , num_divs = 2 )
65+ obs_df <- data.frame (interped $ z )
66+ names(obs_df ) <- paste(' var_' y_out , sep = ' '
67+ obs_df <- cbind(data.frame (DateTime = as.POSIXct(x_out , origin = ' 1970-01-01' obs_df )
68+
69+ # Use model to define X-axis plotting extent for both graphs
70+ xaxis <- get_xaxis(model_df [,1 ])
71+
72+ y.text = y_out [1 ]+ diff(range(y_out ))* 0.05 # note, reference will ALWAYS be surface for compare to field data
73+ if (missing(col_lim ))
74+ col_lim = range(data [, - 1 ], na.rm = TRUE )
75+
76+ .plot_df_heatmap(obs_df , bar_title = .unit_label(nc_file ,var_name ), overlays = c(points(x = x ,y = y ),text(x_out [1 ],y = y.text ,' Observed' pos = 4 , offset = 1 )), xaxis = xaxis , col_lim = col_lim )
77+
78+ .plot_df_heatmap(model_df , bar_title = .unit_label(nc_file ,var_name ), overlays = text(x_out [1 ],y = y.text ,' Modeled' pos = 4 , offset = 1 ), xaxis = xaxis , col_lim = col_lim )
79+
80+ par(start_par )# set PAR back to what it started at
81+ if (is.character(fig_path ))
82+ dev.off()
7983}
0 commit comments