7 Composite all PPIs
With all the preprocessing done, we can now make composites of the PPIs to ‘solve’ the overlapping areas of the Den Helder and Herwijnen radars. While we’re at it, we will also render out these composites for VIR, to visualise the en-masse take-off of birds.
7.2 Generating the composite PPIs
We will loop over all the Herwijnen and Den Helder PPIs we have generated so far to create composite PPIs for all included parameters. We generate these composites at a 500m, 1000m and 2000m resolution to later test at which resolution the models perform best. Additionally we will save the composites as a bunch of .png files that we can then use separately or inclued in the animated GIF below.
source("R/comp_ppi.R")
hrw_ppis <- Sys.glob(file.path("data/processed/final-ppis", "*NL62*"))
dhl_ppis <- Sys.glob(file.path("data/processed/final-ppis", "*NL61*"))
generate_composites <- function(hrw_ppis, dhl_ppis, res, maxrange, output_dir_ppi, output_dir_vis) {
basemap <- NULL
for (i in seq_along(hrw_ppis)) {
ppi_hrw <- readRDS(hrw_ppis[i])
print(hrw_ppis[i])
ppi_dhl <- readRDS(dhl_ppis[i])
print(dhl_ppis[i])
# Set all columns to NA if further than maxrange from radar
ppi_hrw$data@data[ppi_hrw$data@data$dist_radar > maxrange, ] <- NA
ppi_dhl$data@data[ppi_dhl$data@data$dist_radar > maxrange, ] <- NA
params <- c("VIR", "VID", "R", "overlap", "eta_sum", "eta_sum_expected", "dist_radar", "class", "land",
"urban", "agricultural", "semiopen", "forests", "wetlands", "waterbodies", "dist_urban", "human_pop",
"wb_area_id", "wb_area_nr", "ptt_route", "wb_area_ha", "wb_total_biomass", "ptt_total_biomass", "total_biomass",
"wb_weighted_mean_weight", "ptt_weighted_mean_weight", "weighted_mean_weight",
"wb_total_crs", "ptt_total_crs", "total_crs",
"wb_weighted_mean_crs", "ptt_weighted_mean_crs", "weighted_mean_crs",
"wb_total_birds", "ptt_total_birds", "total_birds")
# All mean methods except for factors and urban area (set to max), because we want to strictly filter out fireworks
methods <- c("mean", "mean", "mean", "mean", "mean", "mean", "min", "min", "min",
"max", "mean", "mean", "mean", "mean", "mean", "mean", "mean",
"factor", "factor", "factor", "factor", "mean", "mean", "mean",
"mean", "mean", "mean",
"mean", "mean", "mean",
"mean", "mean", "mean",
"mean", "mean", "mean")
cppi <- comp_ppi(list(ppi_hrw, ppi_dhl), param = params, method = methods, res = c(res, res), coverage = "count")
# Set rain and background pixels to NA
cppi$data$VIR[cppi$data$class < 2] <- NA
# Add coordinates
coords_cppi <- raster::coordinates(cppi$data)
cppi$data$x <- coords_cppi[, 1]
cppi$data$y <- coords_cppi[, 2]
# Add pixel ID
cppi$data@data %>%
mutate(pixel = row_number()) -> cppi$data@data
# Solve different factors
solve_factors <- function(x) {
r <- x[1]
if (is.na(x[1]) && is.na(x[2])) { r <- NA }
if (is.na(x[1]) && !is.na(x[2])) { r <- x[2] }
if (!is.na(x[1]) && is.na(x[2])) { r <- x[1] }
if (!is.na(x[1]) && !is.na(x[2])) {
if (x[1] != x[2]) {
r <- NA
} else {
r <- x[1]
}
}
r
}
for (j in which(methods == "factor")) {
cppi$data@data[params[j]] <- apply(cppi$data@data[params[j]], 1, solve_factors)
}
saveRDS(cppi, file = paste(output_dir_ppi, res, "m/", strftime(ppi_hrw$datetime, format = "%Y%m%d%H%M"), ".RDS", sep = ""))
if (i == 1) {
basemap <- download_basemap(cppi, alpha = 0.3)
}
cppi$data$VIR <- log10(cppi$data$VIR)
cppi$data$VIR[is.na(cppi$data$VIR)] <- 0
bioRad::map(cppi, map = basemap, radar_size = 1, xlim = c(3.1, 6.8), ylim = c(51, 54), zlim = c(0, 4.5),
palette = viridis(256, option = "viridis", alpha = 0.6)) +
labs(title = "Fireworks NYE 2017-2018",
subtitle = paste(ppi_hrw$datetime, ' UTC', sep = ""))
ggsave(paste(output_dir_vis, res, "m/", strftime(ppi_hrw$datetime, format = "%Y%m%d%H%M"), ".png", sep = ""))
}
}
r <- lapply(c(500, 1000, 2000), function(x) { generate_composites(hrw_ppis, dhl_ppis, res = x, maxrange = 66000,
output_dir_ppi = "data/processed/composite-ppis/",
output_dir_vis = "data/plots/vir-ppis/")})We should also generate the composite PPIs for the baseline PPIs using the same procedure.
hrw_ppis <- Sys.glob(file.path("data/processed/final-ppis-baseline", "NLHRW*"))
dhl_ppis <- Sys.glob(file.path("data/processed/final-ppis-baseline", "NLDHL*"))
extract_dt <- function(filepath) {
file <- basename(filepath)
str_split(basename(file), "_")[[1]][3]
}
hrw_dts <- unlist(lapply(hrw_ppis, extract_dt))
dhl_dts <- unlist(lapply(dhl_ppis, extract_dt))
hrw_selection <- sort(hrw_ppis[hrw_dts %in% dhl_dts])
dhl_selection <- sort(dhl_ppis[dhl_dts %in% hrw_dts])
r <- lapply(c(500, 1000, 2000), function(x) { generate_composites(
hrw_selection, dhl_selection, res = x, maxrange = 66000,
output_dir_ppi = "data/processed/composite-ppis-baseline/",
output_dir_vis = "data/plots/vir-ppis-baseline/")})7.3 Visualising the composites
In order to run this comparatively simple conversion to a GIF of the PPIs, it may be necessary to increase the memory available to ImageMagick by changing the memory policy in /etc/ImageMagick-6/policy.xml. Additionally, we will reduce the resolution of the animated GIF and we’ll stick to the 1000m resolution PPIs.
width <- 800
composites <- list.files(path = "data/plots/vir-ppis/1000m/", pattern = "*.png", full.names = TRUE)
process_image <- function(img_path) {
image_read(img_path) %>%
image_resize(geometry_size_pixels(width = width)) %>%
image_write(path = paste(tools::file_path_sans_ext(img_path), ".jpeg", sep = ""), format = "jpeg")
}
for (file in list.files(path = "data/plots/vir-ppis/", pattern = "*.png", full.names = TRUE)) {
process_image(file)
}
list.files(path = "data/plots/vir-ppis/", pattern = "*.jpeg", full.names = T) %>%
purrr::map(image_read) %>%
image_join() %>%
image_animate(fps = 2) %>%
image_write("data/plots/vir-ppis/NYE-2017-2018.gif")