Skip to main content

Fireworks

11 Figure 1: Data overview

Most of Figure 1 in the paper is made using QGIS and Adobe Illustrator, but here we export VIR and total_biomass estimates. Furthermore, we use the land use map that is stored in data/processed/landuse/landuse_hrw_reclassified.tif, so that doesn’t need exporting anymore.

11.1 Processing environment

11.2 Workflow

We load the radar PPIs.

hrw <- readRDS("data/processed/final-ppis/RAD_NL62_VOL_NA_201712312305_ODIM.RDS")
dhl <- readRDS("data/processed/final-ppis/RAD_NL61_VOL_NA_201712312305_ODIM.RDS")

Composite the PPIs and log-transform VIR and total_biomass for improved visualisation.

cppi <- comp_ppi(list(hrw, dhl), param = c("VIR", "total_crs", "dist_urban"), res = 500, 
                 method = c("mean", "mean", "min"))
cppi$data$VIR <- log10(cppi$data$VIR)
cppi$data$VIR[is.infinite(cppi$data$VIR)] <- 0
cppi$data$VIR[cppi$data$VIR == 0] <- NA
cppi$data$total_crs[cppi$data$total_crs == 0] <- NA
cppi$data$total_crs <- log10(cppi$data$total_crs)

Make plots to confirm the composite and transformation works as intended.

plot(cppi, param = "VIR", zlim = c(0, 10))
plot(cppi, param = "total_crs", zlim = c(0, 6))
plot(cppi, param = "dist_urban", zlim = c(0, 10000))

And now we write out raster files containing these parameters, which we can then stitch together in QGIS and Adobe Illustrator.

raster::writeRaster(as(cppi$data["VIR"], "RasterLayer"), filename = "data/plots/paper/fig1_VIR.tif", overwrite = TRUE)
raster::writeRaster(as(cppi$data["total_crs"], "RasterLayer"), filename = "data/plots/paper/fig1_total_crs.tif", overwrite = TRUE)
raster::writeRaster(as(cppi$data["dist_urban"], "RasterLayer"), filename = "data/plots/paper/fig1_dist_urban.tif", overwrite = TRUE)

To visualize the disturbance event as a movie, let’s generate raster files for each of the timestamps that we can then animate in AI/AE.

hrw_ppis <- Sys.glob(file.path("data/processed/final-ppis", "*NL62*"))
dhl_ppis <- Sys.glob(file.path("data/processed/final-ppis", "*NL61*"))

create_raster <- function(hrw_fp, dhl_fp) {
  hrw <- readRDS(hrw_fp)
  dhl <- readRDS(dhl_fp)
  cppi <- comp_ppi(list(hrw, dhl), param = c("VIR"), res = 500, 
                 method = c("mean"))
  
  cppi$data$VIR <- log10(cppi$data$VIR)
  cppi$data$VIR[is.infinite(cppi$data$VIR)] <- 0
  cppi$data$VIR[cppi$data$VIR == 0] <- NA
  
  filename <- paste("data/plots/vir-ppis-raster/", strftime(hrw$datetime, format = "%Y%m%d%H%M"), ".tif", sep = "")
  
  raster::writeRaster(as(cppi$data["VIR"], "RasterLayer"), filename = filename, overwrite = TRUE)
}

mapply(create_raster, hrw_ppis, dhl_ppis)