Appalachicola exploration
Below is an exploration of DO and metabolism data from the Cat Point station at Appalachicola, using data from 2006 to 2019. It incorporates some changes to the WtRegDO package on the sine branch. Specifically, the meteval function was updated to include anomalous estimates by month and changes the required data inputs if using different gas exchanges methods (i.e., air temperature is no longer a requirement if using the Wanninkhof method).
knitr::opts_chunk$set(echo = TRUE, message = F, echo = T, warning = F, fig.path = '')
# libraries
box::use(
WtRegDO[ecometab, wtreg, meteval],
here[...],
doParallel[registerDoParallel],
parallel[detectCores],
dplyr[...],
tidyr[...],
SWMPr[import_local, qaqc, comb],
ggplot2[...],
tibble[enframe],
lubridate[...],
ggforce[facet_zoom],
oce[...],
patchwork[...]
)
# base ggplot theme
thm <- theme_minimal(base_size = 14) +
theme(
legend.title = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
strip.background = element_blank(),
legend.position = 'top'
)Data prep
This code chunk does the following:
- Imports water quality data from cat point and weather data from east bay.
- Combines the files at hourly time step and preps for detiding.
- Detides observed DO using window widths from Jill’s presentation.
- Calculates open water metabolism estimates on observed and detided DO.
- Saves the data files use below.
The script is not run in this Rmd file since it takes a while, but it is shown for reproducibility.
# import apacp data, clean up with qaqc
apacpwq <- import_local(station_code = 'apacpwq', path = here('data-raw/284523.zip'))
apacpwq <- qaqc(apacpwq, qaqc_keep = c('0', '1', '2', '3', '4', '5'))
apaebmet <- import_local(station_code = 'apaebmet', path = here('data-raw/284523.zip'))
apaebmet <- qaqc(apaebmet, qaqc_keep = c('0', '1', '2', '3', '4', '5'))
# combine wx and wq, select/rename relevant columns
apa <- comb(apacpwq, apaebmet, timestep = 60, method = 'union') %>%
select(
DateTimeStamp = datetimestamp,
Temp = temp,
Sal = sal,
DO_obs = do_mgl,
WSpd = wspd,
BP = bp,
Tide = depth,
PAR = totpar
)
# site metadata
locs <- SWMPr::stat_locs %>%
filter(station_code == 'apacp')
lat <- locs$latitude
long <- locs$longitude
tz <- attr(apa$DateTimeStamp, 'tzone')
# setup parallel backend
ncores <- detectCores()
registerDoParallel(cores = ncores - 1)
# get evalcor results to check sun angle/tidal height correlations
cordat <- evalcor(apa, tz, lat, long, progress = T, plot = F)
apa$cordat <- cordat
apaobs <- apa %>%
filter(!is.na(Tide))
# detiding
apadtd <- wtreg(apaobs, tz = tz, lat = lat, long = long, sine = T, gasex = 'Wanninkhof', wins = list(3, 12, 0.2), progress = T, parallel = TRUE)
# metablism on observed and detided
apaobseco <- ecometab(apaobs, tz = tz, lat = lat, long = long, DO_var = 'DO_obs', gasex = 'Wanninkhof')
apadtdeco <- ecometab(apadtd, tz = tz, lat = lat, long = long, DO_var = 'DO_nrm', gasex = 'Wanninkhof')
# save the files
save(apaobs, file = here('data/apaobs.RData'))
save(apadtd, file = here('data/apadtd.RData'))
save(apaobseco, file = here('data/apaobseco.RData'))
save(apadtdeco, file = here('data/apadtdeco.RData'))Evaluate anomalies and PAR
Load the files saved from above.
load(file = here('data/apaobs.RData'))
load(file = here('data/apadtd.RData'))
load(file = here('data/apaobseco.RData'))
load(file = here('data/apadtdeco.RData'))Evaluate the anomalous estimates using the meteval function. The new output is a list with elements named cmp for an assessment of the complete time series and mos for an assessment by month.
metevalobs <- meteval(apaobseco)
metevaldtd <- meteval(apadtdeco)
metevalobs## $cmp
## meanPg sdPg anomPg meanRt sdRt anomRt
## 1 87.8302 121.3034 15.76074 -99.92574 126.9163 13.29219
##
## $mos
## month DOcor Pgcor Rtcor meanPg sdPg anomPg
## 1 01 0.087982234 -0.097237542 -0.25595627 36.40752 55.98628 22.977346
## 2 02 -0.040758116 -0.117197547 -0.20832784 50.00434 77.96171 19.191919
## 3 03 -0.048923067 0.061525532 -0.03071626 40.91978 88.43726 23.913043
## 4 04 -0.057263722 -0.004342643 -0.07195823 48.92897 111.90752 23.145401
## 5 05 -0.070241828 -0.245350213 -0.27459306 115.84581 131.21467 14.697406
## 6 06 -0.004416874 -0.329797411 -0.27108494 156.05554 156.62486 13.008130
## 7 07 -0.111209715 -0.381278973 -0.29977834 158.40810 165.46500 13.612565
## 8 08 -0.041274781 -0.490888594 -0.32815730 157.17640 140.57877 9.433962
## 9 09 0.018298544 -0.264254168 -0.25344807 105.67122 123.78441 10.294118
## 10 10 -0.017284551 -0.040145118 -0.05997625 68.47499 70.13867 11.475410
## 11 11 -0.216393690 -0.167999145 -0.27282125 46.19542 44.08784 10.526316
## 12 12 -0.222358393 -0.216505224 -0.35617318 37.68683 56.65281 20.661157
## meanRt sdRt anomRt
## 1 -40.68854 57.07977 19.417476
## 2 -56.03307 94.87794 17.845118
## 3 -52.91373 92.17032 18.322981
## 4 -68.97112 119.11145 18.100890
## 5 -129.69097 131.69865 13.544669
## 6 -176.17531 162.67481 11.111111
## 7 -172.89116 171.74280 13.350785
## 8 -172.59493 145.13455 7.547170
## 9 -122.06751 128.91199 6.127451
## 10 -75.82820 64.68978 7.103825
## 11 -53.57855 51.17175 11.695906
## 12 -42.76015 58.95822 19.008264metevaldtd## $cmp
## meanPg sdPg anomPg meanRt sdRt anomRt
## 1 81.90767 69.50446 8.54498 -93.59255 70.26454 5.95775
##
## $mos
## month DOcor Pgcor Rtcor meanPg sdPg anomPg
## 1 01 0.23106432 -0.03953780 0.04277763 40.92664 52.00938 14.239482
## 2 02 -0.05119115 0.07012787 -0.06206281 59.21723 53.00934 10.437710
## 3 03 0.03739361 0.16038978 0.06442153 52.23732 68.94571 18.322981
## 4 04 0.19707126 0.19650589 0.19229113 73.57002 60.85701 9.198813
## 5 05 0.10362273 0.18016679 0.16929434 111.66845 73.85764 3.170029
## 6 06 0.31593515 -0.01210406 0.02935626 109.12110 77.10458 8.130081
## 7 07 0.19581442 -0.03040586 -0.01234911 109.67155 87.77617 11.518325
## 8 08 0.15093891 -0.06607275 -0.08576549 111.02449 77.63989 8.355795
## 9 09 0.15344058 0.01478369 -0.02785023 104.47858 70.32448 4.656863
## 10 10 0.08788572 0.12008843 0.03105974 80.37817 43.67748 1.366120
## 11 11 -0.15128437 0.01300963 -0.15168801 58.25643 39.95394 7.017544
## 12 12 -0.10302883 0.04820132 -0.01532509 52.78524 38.68320 8.539945
## meanRt sdRt anomRt
## 1 -43.90274 45.00728 10.0323625
## 2 -64.07535 56.15968 8.4175084
## 3 -66.36783 65.62194 11.1801242
## 4 -89.48839 58.93149 5.3412463
## 5 -122.56419 67.83005 2.0172911
## 6 -127.18300 80.72389 6.2330623
## 7 -125.49364 86.70345 9.4240838
## 8 -129.64935 76.40604 5.3908356
## 9 -120.89757 71.13922 2.6960784
## 10 -88.13501 41.01235 0.5464481
## 11 -65.17666 39.79606 4.9707602
## 12 -57.32292 38.87510 6.8870523Compile the complete and monthly results for comparison.
# total anomalies
anomcmp <- list(
'Observed' = metevalobs$cmp,
'Detided' = metevaldtd$cmp
) %>%
enframe %>%
unnest('value') %>%
select(matches('name|anom')) %>%
gather('var', 'val', -name) %>%
mutate(
var = gsub('anom', '', var),
name = factor(name, levels = c('Observed', 'Detided'))
)
# monthly anomalies
anommos <- list(
'Observed' = metevalobs$mos,
'Detided' = metevaldtd$mos
) %>%
enframe %>%
unnest('value') %>%
select(matches('name|anom|month')) %>%
gather('var', 'val', -name, -month)%>%
mutate(
var = gsub('anom', '', var),
name = factor(name, levels = c('Observed', 'Detided'))
)
anomcmp## # A tibble: 4 x 3
## name var val
## <fct> <chr> <dbl>
## 1 Observed Pg 15.8
## 2 Detided Pg 8.54
## 3 Observed Rt 13.3
## 4 Detided Rt 5.96anommos## # A tibble: 48 x 4
## name month var val
## <fct> <chr> <chr> <dbl>
## 1 Observed 01 Pg 23.0
## 2 Observed 02 Pg 19.2
## 3 Observed 03 Pg 23.9
## 4 Observed 04 Pg 23.1
## 5 Observed 05 Pg 14.7
## 6 Observed 06 Pg 13.0
## 7 Observed 07 Pg 13.6
## 8 Observed 08 Pg 9.43
## 9 Observed 09 Pg 10.3
## 10 Observed 10 Pg 11.5
## # ... with 38 more rowsMake a plot to compare the complete and monthly results by anomalies.
ggplot(anomcmp, aes(x = var, y = val, fill = name)) +
geom_bar(stat = 'identity', position = 'dodge', alpha = 0.7) +
thm +
theme(axis.title.x = element_blank()) +
labs(
y = '% anomalies',
subtitle = 'Complete time series'
)
ggplot(anommos, aes(x = month, y = val, fill = name)) +
geom_bar(stat = 'identity', position = 'dodge', alpha = 0.7) +
facet_wrap(~var, ncol = 2) +
thm +
theme(axis.title.x = element_blank()) +
labs(
y = '% anomalies',
subtitle = 'By month'
)
For fun, let’s plot the change in anomalies from observed to detided by month.
toplo <- anommos %>%
spread(name, val) %>%
mutate(
`Delta anomalies` = Detided - Observed
)
ggplot(toplo, aes(x = month, y = `Delta anomalies`)) +
geom_bar(stat = 'identity', position = 'dodge', alpha = 0.7) +
facet_wrap(~var, ncol = 2) +
thm +
theme(axis.title.x = element_blank()) +
labs(
subtitle = 'By month, anomaly change from observed to detided'
)
Compile and summarize observed/detided metabolism results by month.
apaecosum <- list(
Observed = apaobseco,
Detided = apadtdeco
) %>%
enframe %>%
unnest('value') %>%
mutate(
mo = month(Date),
name = factor(name, levels = c('Observed', 'Detided'))
) %>%
select(matches('name|mo|^Pg$|^Rt$')) %>%
gather('var', 'val', -name, -mo) %>%
group_by(name, mo, var) %>%
summarise(
ecomean = mean(val, na.rm = T),
ecolo = t.test(val)$conf.int[1],
ecohi = t.test(val)$conf.int[2],
.groups = 'drop'
)
apaecosum## # A tibble: 48 x 6
## name mo var ecomean ecolo ecohi
## <fct> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 Observed 1 Pg 36.4 30.1 42.7
## 2 Observed 1 Rt -40.7 -47.1 -34.3
## 3 Observed 2 Pg 50.0 41.1 58.9
## 4 Observed 2 Rt -56.0 -66.9 -45.2
## 5 Observed 3 Pg 40.9 31.2 50.6
## 6 Observed 3 Rt -52.9 -63.0 -42.8
## 7 Observed 4 Pg 48.9 36.9 60.9
## 8 Observed 4 Rt -69.0 -81.7 -56.2
## 9 Observed 5 Pg 116. 102. 130.
## 10 Observed 5 Rt -130. -144. -116.
## # ... with 38 more rowsMake a plot of the summaries.
ggplot(apaecosum, aes(x = factor(mo), y = ecomean, color = var)) +
geom_errorbar(aes(ymin = ecolo, ymax = ecohi)) +
geom_point() +
facet_wrap(~name) +
thm +
labs(
y = 'mmol/m2/d',
x = NULL,
subtitle = 'Metabolism by month (+/- 95% CI)'
)
According to our “rules of thumb”, detiding is assumed to be effective if the mean estimates are similar between observed and detided results, whereas the variance is expected to decrease. Long-term averages (e.g., at the month or longer scale) should not change dramatically with detiding because tidal advection acts primarily at shorter time scales. As a result, noise (or variance) of the estimates will be reduced, but overall rates (means) should be similar between results. Changes in averages may suggest over-smoothing has occurred with detiding.
The plots show mean and standard deviation of production and respiration between observed and detided results. The estimates are grouped by season (as quarters) with 1:1 shown as dashed lines and regressions shown as blue lines.
# monthly mean and sd metab
meansdmos <- list(
'Observed' = metevalobs$mos,
'Detided' = metevaldtd$mos
) %>%
enframe %>%
unnest('value') %>%
select(matches('name|sd|mean|month')) %>%
gather('var', 'val', -name, -month)%>%
mutate(
var = gsub('Pg$', ' Pg', var),
var = gsub('Rt$', ' Rt', var),
name = factor(name, levels = c('Observed', 'Detided')),
val = case_when(
var == c('mean Rt') ~ abs(val),
T ~ val
),
qrt = case_when(
month %in% c('01', '02', '03') ~ 'JFM',
month %in% c('04', '05', '06') ~ 'AMJ',
month %in% c('07', '08', '09') ~ 'JAS',
month %in% c('10', '11', '12') ~ 'OND',
),
qrt = factor(qrt, levels = c('JFM', 'AMJ', 'JAS', 'OND'))
) %>%
spread(name, val) %>%
separate(var, c('stat', 'var'))
toplo1 <- meansdmos %>%
filter(stat == 'mean')
toplo2 <- meansdmos %>%
filter(stat == 'sd')
rng1 <- toplo1 %>%
pull(Observed, Detided) %>%
range
rng2 <- toplo2 %>%
pull(Observed, Detided) %>%
range
p1 <- ggplot(toplo1, aes(x = Observed, y = Detided)) +
geom_abline(intercept = 0, slope = 1, linetype = 'dashed') +
geom_point(size = 4, alpha = 0.5, aes(color = qrt)) +
geom_smooth(method = 'lm', se = F) +
facet_wrap(~var) +
coord_cartesian(xlim = rng1, ylim = rng1) +
labs(
subtitle = 'Mean monthly metabolism (mmol/m2/d)',
color = 'Quarter'
)
p2 <- ggplot(toplo2, aes(x = Observed, y = Detided)) +
geom_abline(intercept = 0, slope = 1, linetype = 'dashed') +
geom_point(size = 4, alpha = 0.5, aes(color = qrt)) +
geom_smooth(method = 'lm', se = F) +
facet_wrap(~var) +
coord_cartesian(xlim = rng2, ylim = rng2) +
labs(
subtitle = 'St. dev. monthly metabolism (mmol/m2/d)',
color = 'Quarter'
)
p1 + p2 + plot_layout(ncol = 1, guides = 'collect') &
thm +
theme(
panel.border = element_rect(fill = NA),
panel.grid.major.x = element_line(),
legend.position = 'right',
legend.title = element_text()
) 
The summer months are pretty suspect for the detided results. Looking at the results from the evalcor function suggests that this is a bad time for detiding. But then again, so are the winter months…
toplo <- apaobs %>%
mutate(
dy = yday(DateTimeStamp),
yr = year(DateTimeStamp)
)
ggplot(toplo, aes(x = dy, y = cordat)) +
geom_line() +
geom_hline(yintercept = 0, col = 'lightblue') +
geom_hline(yintercept = 0.25, col = 'tomato1') +
geom_hline(yintercept = -0.25, col = 'tomato1') +
thm +
theme(
panel.grid.major.y = element_blank()
) +
labs(
x = 'Day of year',
y = 'Sun angle/tidal height correlation'
)
We can dig deeper to evaluate associations of production with PAR, both to evaluate differences in detided and observed DO and to understand the potential role of daily PAR intensity with anomalies. First we summarize PAR on the daily scale to combine with the metabolism estimates. For the weighted regression results, there’s a column called solar_period that indicated whether a time step was during the night (sunset) or day (sunrise). We filter only the sunrise time period, then take the average within each day.
apapar <- apadtd %>%
filter(solar_period == 'sunrise') %>%
mutate(
Date = as.Date(DateTimeStamp)
) %>%
group_by(Date) %>%
summarise(
`PAR mean` = mean(PAR, na.rm = T) * day_hrs * 4, # convert PAR to mmol/m2/d (hourly obs is mmol/m2/0.25hr, so x 4 x day hours)
.groups = 'drop'
)
apapar## # A tibble: 58,636 x 2
## Date `PAR mean`
## <date> <dbl>
## 1 2006-01-01 8141.
## 2 2006-01-01 8141.
## 3 2006-01-01 8141.
## 4 2006-01-01 8141.
## 5 2006-01-01 8141.
## 6 2006-01-01 8141.
## 7 2006-01-01 8141.
## 8 2006-01-01 8141.
## 9 2006-01-01 8141.
## 10 2006-01-01 8141.
## # ... with 58,626 more rowsA plot of average daily PAR over time.
ggplot(apapar, aes(x = Date, y = `PAR mean`)) +
geom_line() +
thm +
labs(
x = NULL,
y = 'PAR (mmol/m2/d)'
)
Now we combine with the metabolism observed and detided metabolism estimates.
apaeco <- list(
Observed = apaobseco,
Detided = apadtdeco
) %>%
enframe %>%
unnest('value') %>%
left_join(apapar, by = 'Date') %>%
mutate(
name = factor(name, levels = c('Observed', 'Detided'))
)
apaeco## # A tibble: 117,314 x 8
## name Date Pg Rt NEM Pg_vol Rt_vol `PAR mean`
## <fct> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Observed 2005-12-31 NA NA NA NA NA NA
## 2 Observed 2006-01-01 118. -168. -49.2 54.6 -77.2 8141.
## 3 Observed 2006-01-01 118. -168. -49.2 54.6 -77.2 8141.
## 4 Observed 2006-01-01 118. -168. -49.2 54.6 -77.2 8141.
## 5 Observed 2006-01-01 118. -168. -49.2 54.6 -77.2 8141.
## 6 Observed 2006-01-01 118. -168. -49.2 54.6 -77.2 8141.
## 7 Observed 2006-01-01 118. -168. -49.2 54.6 -77.2 8141.
## 8 Observed 2006-01-01 118. -168. -49.2 54.6 -77.2 8141.
## 9 Observed 2006-01-01 118. -168. -49.2 54.6 -77.2 8141.
## 10 Observed 2006-01-01 118. -168. -49.2 54.6 -77.2 8141.
## # ... with 117,304 more rowsPlotting Pg vs PAR mean may be telling.
ggplot(apaeco, aes(x = `PAR mean`, y = Pg)) +
geom_point(size = 0.5, alpha = 0.6) +
facet_wrap(~ name, ncol = 2) +
thm +
geom_smooth(method = 'lm') +
labs(
x = 'PAR (mmol/m2/d)',
y = 'Pg (mmol/m2/d)'
)
Another way to look at it is binning the PAR values into categories to see if there’s an association of anomalous values with different ranges of PAR. Below we bin PAR into thirds by the different summary metrics.
toplo <- apaeco %>%
select(name, Date, Pg, `PAR mean`) %>%
filter(!is.na(`PAR mean`)) %>%
mutate(
`PAR level` = cut(`PAR mean`, breaks = c(-Inf, quantile(`PAR mean`, c(0.33, 0.66), na.rm = T), Inf), labels = c('lo', 'md', 'hi'))
) %>%
group_by(name, `PAR level`) %>%
summarise(
`% anomalies` = 100 * sum(Pg < 0, na.rm = T) / n(),
.groups = 'drop'
)
toplo## # A tibble: 6 x 3
## name `PAR level` `% anomalies`
## <fct> <fct> <dbl>
## 1 Observed lo 19.6
## 2 Observed md 11.8
## 3 Observed hi 11.8
## 4 Detided lo 11.3
## 5 Detided md 6.24
## 6 Detided hi 6.28And a plot of the results.
ggplot(toplo, aes(x = `PAR level`, y = `% anomalies`, fill = name)) +
geom_bar(stat = 'identity', position = 'dodge', alpha = 0.7) +
thm
And lastly, repeating the analysis by season because of seasonal variation in PAR and metabolic rates.
toplo <- apaeco %>%
select(name, Date, Pg, `PAR mean`) %>%
filter(!is.na(`PAR mean`)) %>%
mutate(
qrt = quarter(Date),
qrt = factor(qrt, levels = c('1', '2', '3', '4'), labels = c('JFM', 'AMJ', 'JAS', 'OND'))
) %>%
mutate(
`PAR level` = cut(`PAR mean`, breaks = c(-Inf, quantile(`PAR mean`, c(0.33, 0.66), na.rm = T), Inf), labels = c('lo', 'md', 'hi'))
) %>%
group_by(name, qrt, `PAR level`) %>%
summarise(
`% anomalies` = 100 * sum(Pg < 0, na.rm = T) / n(),
.groups = 'drop'
)
ggplot(toplo, aes(x = `PAR level`, y = `% anomalies`, fill = name)) +
geom_bar(stat = 'identity', position = 'dodge', alpha = 0.7) +
facet_wrap(~ qrt) +
thm
Evaluate sun angle and tidal correlations
# site metadata
locs <- SWMPr::stat_locs %>%
filter(station_code == 'apacp')
lat <- locs$latitude
long <- locs$longitude
# get sun elevation angle
apasun <- apaobs %>%
mutate(
sunang = sunAngle(DateTimeStamp, longitude = long, latitude = lat)$altitude,
Date = as.Date(DateTimeStamp)
)
# plot
dts1 <- ymd_hms(c('2009-02-01 00:00:00','2009-03-01 00:00:00'))
dts2 <- as.Date(c('2009-08-01 00:00:00','2009-09-01 00:00:00'))
p <- ggplot(apasun, aes(x = DateTimeStamp, y = sunang)) +
geom_line() +
geom_hline(yintercept = 0, color = 'lightblue') +
theme_bw() +
labs(
x = NULL,
y = 'Sun elevation (degrees)'
)
p + facet_zoom(x = DateTimeStamp >= dts1[1] & DateTimeStamp <= dts1[2], zoom.size = 1)
p + facet_zoom(x = DateTimeStamp >= dts2[1] & DateTimeStamp <= dts2[2], zoom.size = 1)
Maybe we should foor the sun elevation at zero?
Now let’s look at the tidal harmonics using the oce R package and the tidem function to predict major tidal harmonics.
datsl <- as.sealevel(elevation = apaobs$Tide, time = apaobs$DateTimeStamp)
plot(datsl)
We want to predict the major diurnal and semidiurnal components, plus the solar annual constituent.
# tidal components to estimate
constituents <- c('M2', 'S2', 'K1', 'O1', 'SA')
# loop through tidal components, predict each with tidem
preds <- sapply(constituents, function(x){
mod <- tidem(t = datsl, constituent = x)
pred <- predict(mod)
pred - mean(pred)
})Here are the major diurnal/semidiurnal components for the same month slices from the sun angle plot.
# prep for plot
toplo <- preds %>%
data.frame %>%
mutate(DateTimeStamp = apaobs$DateTimeStamp) %>%
gather('component', 'estimate', -DateTimeStamp) %>%
mutate(component = factor(component, level = c('Estimated', constituents)))
toplo1 <- toplo %>%
filter(DateTimeStamp >= dts1[1] & DateTimeStamp <= dts1[2])
toplo2 <- toplo %>%
filter(DateTimeStamp >= dts2[1] & DateTimeStamp <= dts2[2])
ggplot(toplo1, aes(x = DateTimeStamp, y = estimate, group = component)) +
geom_line() +
geom_hline(yintercept = 0, color = 'lightblue') +
facet_wrap(~component, ncol = 1) +
thm +
labs(
x = NULL,
y = 'Height (m)'
)
ggplot(toplo2, aes(x = DateTimeStamp, y = estimate, group = component)) +
geom_line() +
geom_hline(yintercept = 0, color = 'lightblue') +
facet_wrap(~component, ncol = 1) +
thm +
labs(
x = NULL,
y = 'Height (m)'
)
Looking only at the solar annual component.
toplo3 <- toplo %>%
filter(component %in% 'SA')
ggplot(toplo3, aes(x = DateTimeStamp, y = estimate, group = component)) +
geom_line() +
geom_hline(yintercept = 0, color = 'lightblue') +
facet_wrap(~component, ncol = 1) +
thm +
labs(
x = NULL,
y = 'Height (m)'
)
Half-baked conclusions
- Detiding removes more anamolous estimates in the cooler months.
- Mid-summer reduction Pg is suspect, likely related to over-smoothing from weighted regression
- More anamolous estimates seem to occur at lower PAR values, detiding occurs across all PAR values
- Need to better understand the
evalcorfunction - do we floor the sun elevation angles at zero? Is it related to the SA component?