# read in data
df_wq <- read_wq_data(monthly = FALSE)

df_wq$DissAmmonia <- replace_rl(df_wq, 'DissAmmonia')
df_wq$DissNitrateNitrite <- replace_rl(df_wq, 'DissNitrateNitrite')
df_wq$DissOrthophos <- replace_rl(df_wq, 'DissOrthophos')

df_wq <- clean_df(df_wq)

Nutrients Included Are: * ‘Chlorophyll’ (Chlorophyll) * ‘DissNitrateNitrite’ (Dissolved Nitrate/Nitrite) * ‘DissAmmonia’ (Dissolved Ammonia) * ‘Salinity’ (Salinity) * ‘Secchi’ (Secchi) * ‘Temperature’ (Temperature) * ‘TotPhos’ (Total Phosphorous)

Check Temporal Coverage

# check temporal coverage
plt <- check_temporal_coverage(df_wq)

plt

Core stations are:

Subset out Relevant Stations

station_list <- c('P8','MD10','MD10A','D8','D7','D6','D41','D4','D28A','D26','C3','C3A','C10','C10A')
df_wq <- df_wq[df_wq$Station %in% station_list,]

# check temporal coverage
plt <- check_temporal_coverage(df_wq)

plt

Check Spatial Coverage

map <- create_station_map(df_wq)
map

Check if C3/C3A, MD10/MD10A, and C10/C10A can be combined

analytes <- unique(df_wq$Analyte)
stations_list <- list(c('MD10','MD10A'),c('C10','C10A'),c('C3','C3A'))

for (stations in stations_list){
  cat('\n##', stations, '{.tabset .tabset-fade .tabset-pills}')
  for (analyte in analytes){
    cat('\n###', analyte, '\n')
    df_check <- df_wq %>% filter(Station %in% stations, Analyte == analyte)
    plt <- ggplot(df_check) +
      geom_line(aes(Date, Value, color = Station)) +
      ylab(analyte)
  
    plot(plt)
    cat('\n')
  }
}

MD10 MD10A

Temperature

Chlorophyll

DissOrthophos

Secchi

Salinity

DissNitrateNitrite

DissAmmonia

C10 C10A

Temperature

Chlorophyll

DissOrthophos

Secchi

Salinity

DissNitrateNitrite

DissAmmonia

C3 C3A

Temperature

Chlorophyll

DissOrthophos

Secchi

Salinity

DissNitrateNitrite

DissAmmonia

Looks good to me, so will combine the stations.

df_wq <- combine_wq_stations(df_wq)

# check temporal coverage
plt <- check_temporal_coverage(df_wq)
plt

Calculate annual indices (first pass).

df_wq <- df_wq %>%
region_assigner(analysis='annual')
## Loading required package: sf
## Warning: package 'sf' was built under R version 4.1.3
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

df_wq$Year <- lubridate::year(df_wq$Date)

df_wq$Value <- round(df_wq$Value, 3)


df_wq_regions <- df_wq %>%
    dplyr::group_by(Analyte, Year, Region) %>%
    dplyr::summarize(Value=mean(Value, na.rm=TRUE), .groups='drop')

df_wq_noregions <- df_wq %>%
    dplyr::group_by(Analyte, Year) %>%
    dplyr::summarize(Value=mean(Value, na.rm=TRUE), .groups='drop')


# ggplot(data=df_wq, aes(Year, Value, color = Analyte)) +
#   geom_point() +
#   facet_wrap( ~ Analyte, scales='free_y')
df_stations <- subset(df_wq, select = c(Station, Latitude, Longitude))
df_stations <- distinct(df_stations)

# subset out coords for merged stations
df_stations <- df_stations[!(df_stations$Latitude %in% c(38.04381, 37.67575, 38.34575)),]
df_wq_wide_regions <- pivot_wider(data = df_wq_regions, id_cols = c(Year,Region),
                                      names_from = Analyte, values_from = Value)

df_wq_wide_noregions <- pivot_wider(data = df_wq_noregions, id_cols = Year,
                                      names_from = Analyte, values_from = Value)

write_csv(df_wq_wide_regions, 'data/annual_averages/nutrient_data_yearly_regions.csv')
write_csv(df_wq_wide_noregions, 'data/annual_averages/nutrient_data_yearly_noregions.csv')
write_csv(df_stations, 'data/stations/stations_nutrients_annual.csv')