In jwbannister/owensReports:

airsci_loc <- Sys.getenv("R_AIRSCI")
suppressMessages(devtools::load_all(airsci_loc))
load_all("~/code/owensReports")
library(pander)
library(tidyverse)
library(lubridate)
panderOptions('knitr.auto.asis', FALSE)
panderOptions('keep.line.breaks', TRUE)
panderOptions('table.alignment.rownames', 'left')

flux_df <- load_site_data(area, start_date, end_date)
flux_df[flux_df$sand_flux<=0, ]$sand_flux <- 0
area_polys <- pull_onlake_polygons() %>% 
    rename(id2=dca_name) %>%
    mutate(id1=id2, id3=id2)
area_labels <- pull_onlake_labels() %>% 
    mutate(id2=label, id1=label, id3=label)
if (area=='twb2'){
    area_polys <- group_twb2_areas(area_polys)
    area_labels <- group_twb2_areas(area_labels)
}
csc_locs <- load_sites(area, area_polys, end_date) 
csc_locs <- csc_locs %>% arrange(csc) %>% 
    left_join(distinct(select(area_polys, objectid, id1, id2, id3)), by="objectid")
csc_collections <- load_collections(area, start_date, end_date)
# here is where sites can be intentionally excluded from the report
if (TRUE){
    excluded_sites = paste(c(1920, 1921, seq(1938, 1940), seq(1942, 1947),
                             1960, 1963, 1964, 1321))
    flux_df <- flux_df %>% filter(!(csc %in% excluded_sites))
    csc_locs <- csc_locs %>% filter(!(csc %in% excluded_sites))
    csc_collections <- csc_collections %>% filter(!(deployment %in% excluded_sites))
}
# this is a one-off adjustment to check fluxes for sites that weren't accessible 
# prior to May 2019
if (FALSE){
    print("Running false flooded branch...")
    library(pbapply)
    source('~/code/owensReports/scripts/twb2_may_adjust.R')
    new_start <- as.Date('2019-01-01')
    flux_df <- load_site_data(area, new_start, end_date) %>%
        filter(!invalid)
    flux_df[flux_df$sand_flux<=0, ]$sand_flux <- 0
    csc_collections <- load_collections(area, new_start, end_date)

    site_list <- c('1600', '1601', '1607', '1643', '1644', '1645')
    new_collections <- csc_adjust(csc_collections, site_list)
    csc_collections <- csc_collections %>%
        filter(!(csc_data_id %in% new_collections$csc_data_id)) %>%
        filter(!(deployment %in% site_list & dwp_mass == -999))
    csc_collections <- rbind(csc_collections, new_collections)

    a <- unlist(pbapply(flux_df, 1, apply_func, col='dwp_mass'))
    flux_df$new_mass <- as.numeric(a)
    b <- unlist(pbapply(flux_df, 1, apply_func, col='sumpc_total'))
    flux_df$new_sumpc_total <- sapply(b, function(x) x[1])
    geom_factor = 1.2
    flux_df$sand_flux <- (flux_df$new_mass/geom_factor) * 
        (flux_df$sumpc/flux_df$new_sumpc_total)
    flux_df <- filter(flux_df, as.Date(datetime - 1) >= start_date)
    csc_collections <- filter(csc_collections, 
                              as.Date(collection_datetime - 1) >= start_date)
}

source("~/code/owensReports/scripts/sandmass.R")

# report_index is set after daily sand flux calculations
source("~/code/owensReports/scripts/sandflux.R")

# don't run twb2 paired TEOM analysis if too early in the month (interim 
# reports). TEOM data may not be in database and will throw error.
#paired_flag <- ifelse(Sys.Date() < start_date %m+% weeks(2), FALSE, TRUE)
paired_flag <- TRUE
if (area=='twb2' & paired_flag) source("~/code/owensReports/scripts/twb2_pairs.R")

if (area=='twb2') source("~/code/owensReports/scripts/twb2_surface.R")

lidar_flag <- TRUE
if (area=='twb2' & lidar_flag) source("~/code/owensReports/scripts/twb2_lidar.R")

```r
# One-off adjustments specific to report
# CHECK EACH MONTH AND CHANGE IF NECESSARY
url <- "ftp://proc2.airsci.com//changes/"
ftp_login = Sys.getenv("OWENS_REPORTS_FTP")
ftp_output <- RCurl::getURL(url, userpwd = ftp_login, 
                               ftp.use.epsv=FALSE, dirlistonly = TRUE) 
change_list <- strsplit(ftp_output, "\n")[[1]]
change_file <- paste0(area, month(start_date), year(start_date), ".R")
if (change_file %in% change_list){
    file_url <- paste0(url, change_file)
    fl <- RCurl::getBinaryURL(file_url, userpwd = ftp_login)
    writeBin(fl, paste0(tempdir(), change_file))
    source(paste0(tempdir(), change_file))
}

img_files <- vector(mode="character", length=length(report_index))
names(img_files) <- report_index
pair_img_files <- img_files
surface_img_files <- img_files
for (i in report_index){
    img_files[[i]] <- paste0(tempfile(), ".png")
    png(filename=img_files[[i]], width=8, height=4, units="in", res=300)
    gridExtra::grid.arrange(flux_grobs[[i]], ncol=2)
    dev.off()
    if (area=='twb2' & i!='East' & paired_flag){
        # twb2 report also shows TEOM pair maps
        pair_img_files[[i]] <- paste0(tempfile(), ".png")
        png(filename=pair_img_files[[i]], width=8, height=4, units="in", res=300)
        gridExtra::grid.arrange(pair_grobs[[i]], ncol=2)
        dev.off()
    }
    if (area=='twb2'){
        # twb2 report also shows surface data trend plots
        surface_img_files[[i]] <- paste0(tempfile(), ".png")
        png(filename=surface_img_files[[i]], width=8, height=8, units="in", res=300)
        gridExtra::grid.arrange(surface_grobs[[i]]$rh_plot, 
                                surface_grobs[[i]]$rsrh_plot, 
                                surface_grobs[[i]]$clods_plot, ncol=2)
        dev.off()
    }
}

for (i in report_index){
    # build page titles
    if (area=='twb2'){
        ttle <- paste0(" \n# ", i, " TwB2 \n")
    } else{ 
        ttle <- paste0(" \n# ", i, " \n")
    }
    report_header(start_date, end_date, report_date, area)
    cat(ttle)
    # site maps on top of page
    cat(" \n## Monitoring Site Results \n")
    cat(paste0("![](", img_files[[i]], ")"))
        # tables for all other reports
        cat(" \n\n## Daily Sand Flux \n")
        if (area %in% c("channel", "t1a1")){
            area_index <- c('channel'='Channel Areas are', 't1a1'='T1A-1 is')
            cat(paste0("\n\nNote that the ", area_index[area], " not subject ", 
                       "to a specific sand flux threshold. However, this report ",
                       "includes dates when sand flux was observed in excess of ",
                       "1.0 g/cm<sup>2</sup>/day; the most conservative sand flux ",
                       "threshold applied to the approved Best Approved ",
                       "Control Measures.\n"))
        }
        if (area=='twb2'){
            cat(paste0(" \nThe sand flux threshold for maintenance is 0.5 ",
                       "g/cm<sup>2</sup>/day. The reflood threshold is 1.0 ",
                       "g/cm<sup>2</sup>/day (Rule 433, paragraph h.i.).\n"))
        }
        if (area %in% c("brine", "dwm")){
            cat(paste0(" \nErosion threshold to trigger BACM Shallow Flooding ", 
                       "is sand flux of 5.0 grams per square centimeter per day ", 
                       "(Rule 433, paragraph h.*i*.) \n"))
        } 
        if (area=='twb2'){
            # twb2 is broken up into regions
            tbl_flux <- daily_flux %>% filter(id3==i, sand.flux>=0.5) %>% 
                select(id2, csc, date, sand.flux) %>% arrange(desc(sand.flux))
            colnames(tbl_flux) <- c("DCA", "CSC Site", "Date", 
                                    "Sand Flux\\\n(g/cm<sup>2</sup>/day)") 
            tbl_mass <- csc_mass %>% filter(id3==i) %>% 
                select(csc, sand.mass, comment) %>% arrange(csc)
            if (nrow(tbl_mass)>0){
                for (row_num in 1:nrow(tbl_mass)){
                    if(!is.na(tbl_mass$comment[row_num])){
                        tbl_mass$sand.mass[row_num] <- tbl_mass$comment[row_num]
                    }
                }
            }
            tbl_mass <- select(tbl_mass, -comment)
            colnames(tbl_mass) <- c("CSC Site", "Sand Mass (g)") 
            row.names(tbl_mass) <- NULL
        } else{
            # all other reports are broken up by DCA
            tbl_flux <- daily_flux %>% filter(id2==i, sand.flux>=1) %>% 
                select(csc, date, sand.flux) %>% arrange(date, sand.flux) 
            tbl_flux <- tbl_flux[!duplicated(tbl_flux), ]
            colnames(tbl_flux) <- c("CSC Site", "Date", 
                                    "Sand Flux\\\n(g/cm<sup>2</sup>/day)") 
            row.names(tbl_flux) <- NULL
            tbl_mass <- csc_mass %>% filter(id2==i) %>% 
                select(csc, sand.mass, comment) %>% arrange(csc)
            if (nrow(tbl_mass)>0){
                for (row_num in 1:nrow(tbl_mass)){
                    if(!is.na(tbl_mass$comment[row_num])){
                        tbl_mass$sand.mass[row_num] <- tbl_mass$comment[row_num]
                    }
                }
            }
            tbl_mass <- select(tbl_mass, -comment)
            colnames(tbl_mass) <- c("CSC Site", "Sand Mass (g)") 
            row.names(tbl_mass) <- NULL
        }
        split_size <- 1
        if (nrow(tbl_flux)==0){
            if (area=='twb2'){
                cat(" \nNo days with measured sand flux >= 0.5 g/cm<sup>2</sup>/day. \n")
            } else{
                cat(" \nNo days with measured sand flux >= 1.0 g/cm<sup>2</sup>/day. \n")
            }
        } else{
            if (area=='twb2'){
                cat(" \nDays with recorded sand flux >= 0.5 g/cm<sup>2</sup>/day: \n")
            } else{
                cat(" \nDays with recorded sand flux >= 1.0 g/cm<sup>2</sup>/day: \n")
            }
            if (nrow(tbl_flux)<=15){
                tbl_flux_list <- vector(mode='list', length=1)
                tbl_flux_list[[1]] <- tbl_flux
                pandoc.table(as.data.frame(tbl_flux_list[[1]]), split.table=Inf)
            } else{
                # if flux table is too large, break into sections
                split_size <- 2 + (nrow(tbl_flux[-(1:15), ]) %/% 35)
                split_factors <- rep(1, 15) 
                for (k in 2:split_size){
                    split_factors <- c(split_factors, rep(k, 35))
                }
                tbl_flux_list <- split(tbl_flux, split_factors)
                for (j in 1:split_size){
                    if (j > 1){
                        report_header(start_date, end_date, report_date, area)
                        cat(" \n\n## Daily Sand Flux - Continued\n")
                    }
                    row.names(tbl_flux_list[[j]]) <- 1:nrow(tbl_flux_list[[j]])
                    pandoc.table(as.data.frame(tbl_flux_list[[j]]), 
                                 split.table=Inf)
                    if (j < split_size){
                        cat("<p style=\"page-break-after:always;\"></p> \n")
                    }
                }
            }
        }
    # check for comment file and see how long
    url <- "ftp://proc2.airsci.com//comments/"
    ftp_login = Sys.getenv("OWENS_REPORTS_FTP")
    ftp_output <- RCurl::getURL(url, userpwd = ftp_login, 
                                   ftp.use.epsv=FALSE, dirlistonly = TRUE) 
    comment_list <- strsplit(ftp_output, "\n")[[1]]
    comment_file <- gsub(" ", "_", paste0(area, month(start_date), year(start_date), i, ".txt"))
    if (comment_file %in% comment_list){
        file_url <- paste0(url, comment_file)
        fl <- RCurl::getURL(file_url, userpwd = ftp_login)
        comments <- strsplit(fl, "\r\n")[[1]]
        comment_length <- 1 + length(comments)
    } else{
        comment_length <- 1
    }

    row_limit <- if_else(split_size==1, 12, 31)
    flux_length <- ifelse(exists("tbl_flux_list"), 
                          nrow(tbl_flux_list[[split_size]]), 0)
    if (nrow(tbl_mass)>0){
        if (flux_length + nrow(tbl_mass) + comment_length < row_limit){
            cat(" \n\n## Monthly Sand Mass \n")
            pandoc.table(as.data.frame(tbl_mass))
        } else if (flux_length + nrow(tbl_mass) < row_limit){
            cat(" \n\n## Monthly Sand Mass \n")
            pandoc.table(as.data.frame(tbl_mass))
            cat("<p style=\"page-break-after:always;\"></p> \n")
            report_header(start_date, end_date, report_date, area)
            cat(ttle)
        } else{
            cat("<p style=\"page-break-after:always;\"></p> \n")
            report_header(start_date, end_date, report_date, area)
            cat(ttle)
            cat(" \n\n## Monthly Sand Mass \n")
            pandoc.table(as.data.frame(tbl_mass))
        }
    }
    if (comment_file %in% comment_list){
        if (nrow(tbl_mass) + comment_length > 30){
            cat("<p style=\"page-break-after:always;\"></p> \n")
        }
        cat(" \n\n## Comments \n")
        for (n in 1:length(comments)){
            cat(" \n", comments[n], 
                " \n", sep="")
        }
    }
    cat("<p style=\"page-break-after:always;\"></p> \n")

    if (area=='twb2' & i!="East" & paired_flag){
        report_header(start_date, end_date, report_date, area)
        cat(ttle)
        cat(" \n## Paired TEOM PM10 Results\n")
        cat(paste0("![](", pair_img_files[[i]], ")"))
        cat(" \n\n\n\n")
        tbl_pair <- daily_summary %>% filter(id3==i, pm10.delta>1) %>% 
            arrange(desc(pm10.delta)) %>% select(-id3, -status)
        colnames(tbl_pair) <- c("Date", 
                                "Upwind Cross-Tillage\\\nAvg. PM10 (micrograms/m<sup>3</sup>)", 
                                "Downwind Cross-Tillage\\\nAvg. PM10 (micrograms/m<sup>3</sup>)", 
                                "Increase in Avg.\\\n PM10 (micrograms/m<sup>3</sup>)", 
                                "Windspeed at Max Hourly\\\nPM10 Increase (m/s)")
        if (nrow(tbl_pair)>20){
            cat("<p style=\"page-break-after:always;\"></p> \n")
            report_header(start_date, end_date, report_date, area)
            cat(ttle)
            cat(" \n## Paired TEOM PM10 Results\n")
        }
        if (nrow(tbl_pair)==0){
            cat(" \n\n###No significant PM10 increases across area recorded. \n")
        } else{
        pandoc.table(as.data.frame(tbl_pair), split.table=Inf)
        }

        # Comments for Mar 2020 report only
        if (i=='South'){
        cat(" \n After review of PM10 and meteorological data, it is determined that the positive PM10 concentration difference between upwind and downwind TEOMs on 3/2/2020 is due to one hour (19:00). Sand fluxes in the South TwB2 areas were zero at all Sensit sites for the 19:00 hour in question, and were below the maintenance threshold for the entire day. Furthermore, the March ground survey data showed that ridge geometry in the South TwB2 areas meets required measurements. In aggregate, this information is indicative that the PM10 measured at the downwind (T2-1) monitor is due to sources outside the TwB2 areas. This event is noted for informational purposes but does not qualify as an exceedance of the TwB2 requirements per the Stipulated Judgement.\n")
        }

        cat("<p style=\"page-break-after:always;\"></p> \n")
    }
    if (area=='twb2'& exists('surface_df')){
        report_header(start_date, end_date, report_date, area)
        cat(ttle)
        cat(" \n## Surface Survey Data\n")
        coll_period <- surface_df[[i]]$coll_period[!is.na(surface_df[[i]]$coll_period)][1]
        tbl_surface <- surface_df[[i]] %>% 
            filter(!(is.na(rs_avg) & is.na(rh_avg) & 
                     is.na(rs_rh) & is.na(clods))) %>%
            select(id2, site, rs_avg, rh_avg, rs_rh, clods)
        missing_surface <- surface_df[[i]] %>% 
            filter((is.na(rs_avg) & is.na(rh_avg) & 
                    is.na(rs_rh) & is.na(clods)))
        colnames(tbl_surface) <- c("DCA", "Site", "Row Spacing (cm)", 
                                   "Row Height (cm)", "RS/RH", 
                                   "Clod Coverage (%)")
        cat("<h6> \n")
        pandoc.table(as.data.frame(tbl_surface), split.table=Inf)
        cat("</h6> \n")
        cat(" \n\n## Comments \n")
        cat(paste0(" \n Surface survey data was collected ", coll_period, ".\n"))
        cat(" \n RS and RH results are reported as rounded integer numbers. RS/RH ratio reported above is calculated from the actual (decimal) values.\n")
        comment_file_surface <- 
            gsub(" ", "", paste0(area, month(start_date), year(start_date), i, 
                                 "_surface.txt"))
        if (comment_file_surface %in% comment_list){
            file_url <- paste0(url, comment_file_surface)
            fl <- RCurl::getURL(file_url, userpwd = ftp_login)
            comments_surface <- strsplit(fl, "\r\n")[[1]]
            for (k in 1:length(comments_surface)){
                cat(" \n", comments_surface[k], " \n", sep="")
            }
        }
        cat("<p style=\"page-break-after:always;\"></p> \n")
        report_header(start_date, end_date, report_date, area)
        cat(ttle)
        cat(" \n## Surface Survey Data Trend Plots\n")
        cat(paste0("![](", surface_img_files[[i]], ")"))
        cat(" \n\n## Comments \n")
        cat(" \n Temporal figures reflect the average conditions at all current survey sites by area.  Month to month differences can be attributed to a combination of actual differences, as well as differences in the number and locations of survey sites.  The latter two can vary on a monthly basis, due to maintenance (sites not used) or wet surface conditions (lack of safe access).\n")
        cat("<p style=\"page-break-after:always;\"></p> \n")
    }
    if (area=='twb2' & lidar_flag){
        lidar_images <- lidar_files[[i]][!is.na(lidar_files[[i]])]
        if (length(lidar_images)>0){
            report_header(start_date, end_date, report_date, area)
            cat(ttle)
            cat(" \n## LiDAR Roughness Maps\n")
            cat(" \n Maps developed from most recently processed data.\n")
            for (n in 1:length(lidar_images)){
               # knitr::include_graphics(lidar_images[n])
               cat(paste0("<p>![](", lidar_images[n], "){height=8in}</p>"))
               # if ((n %% 2 == 0) & (n!=length(lidar_images))){
                if (n<length(lidar_images)){
                    cat("<p style=\"page-break-after:always;\"></p> \n")
                    report_header(start_date, end_date, report_date, area)
                    cat(ttle)
                    cat(" \n## LiDAR Roughness Maps (Continued)\n")
                    cat(" \n Maps developed from most recently processed data.\n")
                }
            }
            cat("<p style=\"page-break-after:always;\"></p> \n")
        }
    }
}

jwbannister/owensReports documentation built on Aug. 8, 2020, 11 p.m.