Merge pull request #393 from barthoekstra/composite-ppi

Composite PPIs for multiple parameters at once

R/composite_ppi.R

@@ -5,13 +5,18 @@
 #' radars.
 #'
 #' @inheritParams integrate_to_ppi
-#' @param x A list of \code{ppi} objects.
-#' @param param Scan parameter to composite.
-#' @param method string. Compositing method, one of "mean", "min", "max" or "idw"
+#' @param x A list of `ppi` objects.
+#' @param param The scan parameter name(s) to composite. An atomic vector of character strings
+#' can be provided to composite multiple scan parameters at once. To composite all available
+#' scan parameters use 'all' (default).
+#' @param method string. Compositing method, one of "mean", "min", "max" or "idw". Provide a list of methods
+#' names of length(param) to apply different methods to each of the parameters.
 #' @param idw_max_distance numeric. Maximum distance from the radar to consider in
 #' inverse distance weighting. Measurements beyond this distance will have a
 #' weighting factor of zero.
-#' @param idp numeric. inverse distance weighting power
+#' @param idp numeric. inverse distance weighting power.
+#' @param coverage logical. When TRUE adds an additional "coverage" parameter to the \code{ppi} with the
+#' number of PPIs covering a single composite \code{ppi} pixel.
 #'
 #' @return A \code{\link[=summary.ppi]{ppi}}.
 #'
@@ -38,7 +43,10 @@
 #' }
 #'
 #' The coordinates system of the returned \code{ppi} is a WGS84
-#' (lat,lon) datum.
+#' (lat, lon) datum, unless a different \code{crs} is provided. If only
+#' \code{res} is provided, but no \code{crs} is set, \code{res} is in
+#' meter units and the origin of the composite \code{ppi} is set to the
+#' mean (lat, lon) location.
 #'
 #' This function is a prototype and under active development
 #'
@@ -63,7 +71,7 @@
 #' # plot the calculated max product on the basemap
 #' map(my_composite, bm)
 #' }
-composite_ppi <- function(x, param = "DBZH", nx = 100, ny = 100, xlim, ylim, res, crs, raster = NA, method = "max", idp = 2, idw_max_distance=NA) {
+composite_ppi <- function(x, param = "all", nx = 100, ny = 100, xlim, ylim, res, crs, raster = NA, method = "max", idp = 2, idw_max_distance = NA, coverage = FALSE) {
   if (FALSE %in% sapply(x, is.ppi)) {
     stop("'composite' expects objects of class ppi only")
   }
@@ -80,19 +88,25 @@ composite_ppi <- function(x, param = "DBZH", nx = 100, ny = 100, xlim, ylim, res
   if (!missing(res)) {
     assert_that(is.numeric(res))
     assert_that(length(res) <= 2)
+    t_res <- res
   } else {
-    res <- NA
+    t_res <- NULL
   }
-
   # check crs argument as in raster::raster()
   if (!missing(crs)) {
-    crs <- CRS(as.character(raster::projection(crs)))
-  }
-  else {
-    crs <- CRS("+proj=longlat +datum=WGS84")
+    t_crs <- CRS(as.character(raster::projection(crs)))
+  } else {
+    t_crs <- NULL
   }
+  if (!all(method %in% c("max", "min", "mean", "idw"))) stop("'method' should be one or multiple of 'max', 'mean', 'min' or 'idw'")
+  if (length(method) != length(param) & length(method) != 1) stop("'method' should be of length 1 or length(param)")
+  if (!is.logical(coverage)) stop("'coverage' should be a logical")
 
+  if (length(param) == 1 && param == "all") {
+    param <- names(x[[1]]$data)
+  }
   ppis <- lapply(x, `[.ppi`, i = param)
+
   lons <- sapply(ppis, function(x) x$geo$bbox["lon", ])
   lats <- sapply(ppis, function(x) x$geo$bbox["lat", ])
   if(!missing(xlim)) lons <- xlim
@@ -106,53 +120,92 @@ composite_ppi <- function(x, param = "DBZH", nx = 100, ny = 100, xlim, ylim, res
   )
 
   if (!are_equal(raster, NA)) {
-    r <- raster(raster)
+    r <- raster::raster(raster)
   } else {
-    if (missing(res) | is.na(res)) {
-      r <- raster(ncols = nx, nrows = ny, ext = raster::extent(c(min(lons),max(lons),min(lats),max(lats))), crs = crs)
-    }
-    else {
-      r <- raster(ncols = nx, nrows = ny, ext = raster::extent(c(min(lons),max(lons),min(lats),max(lats))), crs = crs, res = res)
+    d_crs <- CRS("+proj=longlat +datum=WGS84")
+    if (!is.null(t_res) && !is.null(t_crs)) {
+      r <- raster(ext = raster::extent(c(min(lons), max(lons), min(lats), max(lats))), crs = t_crs, resolution = t_res)
+    } else if (!is.null(t_crs) && is.null(t_res)) {
+      r <- raster(ncols = nx, nrows = ny, ext = raster::extent(c(min(lons), max(lons), min(lats), max(lats))), crs = t_crs)
+    } else if (is.null(t_crs) && !is.null(t_res)) {
+      r <- raster(ext = raster::extent(c(min(lons), max(lons), min(lats), max(lats))), crs = d_crs)
+      t_crs <- CRS(paste0("+proj=aeqd +units=m +ellps=WGS84 +lat_0=", mean(lats), " +lon_0=", mean(lons)))
+      r <- raster::projectExtent(r, t_crs)
+      raster::res(r) <- t_res
+    } else {
+      r <- raster(ncols = nx, nrows = ny, ext = raster::extent(c(min(lons), max(lons), min(lats), max(lats))), crs = d_crs)
     }
   }
 
   # initialize all values of the grid to NA
-  raster::values(r) <- NA
-  spGrid = as(r,'SpatialGridDataFrame')
+  suppressWarnings(r <- raster::setValues(r, NA))
+  spGrid = as(r, 'SpatialGridDataFrame')
   names(spGrid@data) <- names(ppis[[1]]$data)[1]
 
+  if (coverage) {
+    ppis <- lapply(ppis, function(x) {x$data$coverage <- 1; return(x)})
+    param <- c(param, "coverage")
+  }
+
   # merge
-  projs <- suppressWarnings(sapply(
-    ppis,
+  projs <- sapply(ppis,
     function(x) {
       over(
-        spTransform(
-          spGrid,
-          CRS(proj4string(x$data))
+        suppressWarnings(
+          spTransform(
+            spGrid,
+            CRS(proj4string(x$data))
+          )
         ),
-        x$data
+        x$data[param]
       )
     }
-  ))
+  )
 
-  if(method == "max") spGrid@data[, 1] <- do.call(function(...) pmax(..., na.rm = TRUE), projs)
-  if(method == "min") spGrid@data[, 1] <- do.call(function(...) pmin(..., na.rm = TRUE), projs)
-  if(method == "mean") as.data.frame(projs) %>% rowMeans(na.rm=T) -> spGrid@data[, 1]
-  if(method == "idw"){
-    brick_data = raster::brick(raster::brick(spGrid),nl=length(projs))
-    brick_weights = brick_data
-    #weights<-raster::pointDistance(as.matrix(data.frame(x=lons.radar,y=lats.radar)), coordinates(raster(spGrid)),lonlat=T)
-    for(i in 1:length(projs)){
-      brick_data <- raster::setValues(brick_data, projs[[i]], layer=i)
-      latlon.radar <- unique(data.frame(lat=c(lats.radar), lon=c(lons.radar)))
-      weights<-raster::pointDistance(as.matrix(data.frame(x=latlon.radar$lon,y=latlon.radar$lat))[i,], coordinates(raster(spGrid)),lonlat=T)
-      if(!is.na(idw_max_distance)) weights[weights>idw_max_distance]=NA
-      weights = 1/(weights^idp)
+  for (p in param) {
+    if (p == "coverage") next()
+    if (length(param) > 1) {
+      merged <- projs[p, ]
+    } else {
+      merged <- projs
+    }
 
-      brick_weights <- raster::setValues(brick_weights, weights, layer=i)
+    if (length(method) > 1) {
+      param_method <- method[match(p, param)]
+    } else {
+      param_method <- method
     }
-    spGrid <- as(raster::weighted.mean(brick_data, brick_weights, na.rm=T),"SpatialGridDataFrame")
-    names(spGrid@data) <- names(ppis[[1]]$data)[1]
+
+    if(param_method == "max") spGrid@data[, p] <- do.call(function(...) pmax(..., na.rm = TRUE), merged)
+    if(param_method == "min") spGrid@data[, p] <- do.call(function(...) pmin(..., na.rm = TRUE), merged)
+    if(param_method == "mean") as.data.frame(merged) %>% rowMeans(na.rm=TRUE) -> spGrid@data[, p]
+    if(param_method == "idw"){
+      brick_data <- suppressWarnings(raster::brick(raster::brick(spGrid), nl = length(merged)))
+      brick_weights <- brick_data
+      #weights<-raster::pointDistance(as.matrix(data.frame(x=lons.radar,y=lats.radar)), coordinates(raster(spGrid)),lonlat=T)
+      for(i in 1:length(merged)){
+        brick_data <- raster::setValues(brick_data, merged[[i]], layer=i)
+        latlon.radar <- unique(data.frame(lat = c(lats.radar), lon = c(lons.radar)))
+        if (is.null(t_res)) {
+          weights <- suppressWarnings(raster::pointDistance(as.matrix(data.frame(x = latlon.radar$lon, y = latlon.radar$lat))[i, ], coordinates(raster(spGrid)), lonlat = TRUE))
+        } else {
+          d <- data.frame(lon = latlon.radar$lon, lat = latlon.radar$lat)
+          coordinates(d) <- c("lon", "lat")
+          proj4string(d) <- d_crs
+          proj.radar <- as.data.frame(spTransform(d, t_crs))
+          weights <- suppressWarnings(raster::pointDistance(as.matrix(data.frame(x = proj.radar$lon, y = proj.radar$lat))[i, ], coordinates(raster(spGrid)), lonlat = FALSE))
+        }
+        if(!is.na(idw_max_distance)) weights[weights > idw_max_distance] <- NA
+        weights <- 1 / (weights ^ idp)
+        brick_weights <- raster::setValues(brick_weights, weights, layer=i)
+      }
+      spGrid@data[, p] <- as.vector(raster::weighted.mean(brick_data, brick_weights, na.rm = TRUE))
+    }
+  }
+
+  if (coverage) {
+    cov <- !is.na(do.call("cbind", projs["coverage", ]))
+    spGrid@data$coverage <- rowSums(cov)
   }
 
   ppi.out <- list(data = spGrid, geo = list(