sosoc / croc

Ocean colour with R
3 stars 1 forks source link

move out compute jobs #8

Open mdsumner opened 7 years ago

mdsumner commented 7 years ago

This can be its own package,


#' MODISA chlorophyll-a in the southern hemisphere. 
#' 
#' Compute chlorophyll data in the L3 bins, using both Johnson and Oceancolor algorithms. 
#' 
#' Use `raadtools::ocfiles(product = "MODISA", varname = "RRS", type = "L3b")` to get all daily MODISA L3 bin RRS files, and
#' explore the `date` column to generate a set of `days` for this function. 
#' @param days dates to calculate
#' @param maxbin a maximum bin to use, for specialist use only
#' @return tibble of bins and chlorophyll-a values
#' @export
#' @importFrom tibble tibble
#' @importFrom progress progress_bar
#' @importFrom dplyr mutate filter select transmute
#' @importFrom rlang .data
chla_compute <- function(days,  maxbin = NULL) {
  ## data frame of all L3 RRS files for MODISA
  files <- ocfiles(time.resolution = "daily", product = "MODISA", varname = "RRS", type = "L3b", bz2.rm = TRUE, ext = "nc") 
  days <- as.Date(as.POSIXct(days, tz = "GMT"))

  files <- files[as.Date(files$date) %in% days, ]

  ## initialize the bin logic for MODISA
  init <- initbin(NUMROWS = 4320)
  ## counts up from the south
  if (is.null(maxbin)) maxbin <- init$totbin/2
  ## accumulate sum and count for mean
  count_rj <- cumul_rj <- numeric(maxbin)
  count_nasa <- cumul_nasa <- numeric(maxbin)

  ## track progress
  pb <- progress_bar$new(total = nrow(files))
  for (ifile in seq(nrow(files))) {
    d <- read_L3_file(files$fullname[ifile], compound_vars = c("Rrs_443", "Rrs_488", "Rrs_555", "Rrs_547")) %>% 
      filter(.data$bin_num <= maxbin) %>% 
      mutate(chla_johnson = chla(.data, sensor = "MODISA", algo = "johnson")) %>% 
      mutate(chla_nasa = chla(.data, sensor = "MODISA", algo = "oceancolor")) %>% 
      dplyr::select(.data$bin_num, .data$chla_johnson, .data$chla_nasa) %>% 
      filter(.data$chla_johnson > 0)

    ## we accumulate outside of the data frame to ensure speed
    ## the bin_num indexes the vector structure directly since the south pole is bin 1
    count_rj[d$bin_num] <- count_rj[d$bin_num] + 1
    cumul_rj[d$bin_num] <- cumul_rj[d$bin_num] + d$chla_johnson

    count_nasa[d$bin_num] <- count_nasa[d$bin_num] + 1
    cumul_nasa[d$bin_num] <- cumul_nasa[d$bin_num] + d$chla_nasa
    ## track progress
    pb$tick()
  }

  tibble(bin_num = seq_len(maxbin), count_nasa = count_nasa, cumul_nasa = cumul_nasa, 
         count_rj = count_rj, cumul_rj = cumul_rj) %>% filter(.data$count_rj > 0)
}
mdsumner commented 7 years ago

Probably regridder too (replace with Jessica method).

# Regrid ocean colour. 
## 
## Regrid NASA ocean colour RRS values to standard Mapped image with the Johnson chl-a algorithm. 
#### @param tgrid target grid
## @param file L3 bin file name with raw RRS wavelengths
## @param agg aggregate values?
## @param NROWS number of L3 rows
## @param fun function to aggregate by if agg is \code{TRUE}
##
## @importFrom rgdal project
## @importFrom raster setValues
## @export
# regridder <- function(tgrid, file, agg = FALSE, fun = mean, NROWS) {
#   .Deprecated()
#   xy <- lonlatFromCell(tgrid, spatial = FALSE)
#   
#   stop("need to use read_l3_file somehow")   #MDS 2017-05-25
#   #NROWS <- readL3(file, vname = "NUMROWS", bins = FALSE)$NUMROWS
#   binmap <- lonlat2bin(xy[,1], xy[,2], NROWS)
#   function(file) {
#     sens <- if(x$NUMROWS == 4320) "MODISA" else "SeaWiFS"
#     if(x$NUMROWS != NROWS)  stop("file doesn't match this regridder\nfile NUMROWS: ", x$NUMROWS, "\nfunction NUMROWS:", NROWS)
#     
#     if (agg) {
#       vals <-  chla(x, sensor = sens, algo = "johnson")
#       
#       binxy <- do.call(cbind, bin2lonlat(x$bin_num, x$NUMROWS))
#       if (!is.na(projection(tgrid)) & !isLonLat(tgrid)) binxy <- project(binxy, projection(tgrid))
#       gcell <- extract(tgrid, binxy, cellnumbers = TRUE)[,"cells"]
#       tvals <- tapply(vals, gcell, fun)
#       
#       r <- tgrid
#       r[unique(gcell)] <- tvals
#       r[!r>0] <- NA
#     } else {
#       r <- setValues(tgrid, chla(x, sensor = sens, algo = "johnson")[match(binmap, x$bin_num)])
#     }
#     r
#   }
# }