Maps of covariates

knitr::opts_chunk$set(cache = TRUE)
library(eucs)
library(raster)
library(sp)

rasters <- read.csv("../../data-cache/rasters.csv", stringsAsFactors = FALSE)

grampians <- subset(map_ibra_southeast, IBRA_SUB_N == "Greater Grampians")
grampians_json <- tempfile(fileext = ".json")
rgdal::writeOGR(grampians, grampians_json, "grampians", driver = "GeoJSON")

se_aus <- subset(
  map_ibra_southeast,
  IBRA_SUB_N %in% unique(
    subset(occupancy_southeast, ibra_subregion_coverage > .5)$ibra_subregion
  )
)
se_aus@data$ID <- 1
se_aus <- aggregate(se_aus, by = "ID")
se_aus_json <- tempfile(fileext = ".json")
rgdal::writeOGR(se_aus, se_aus_json, "se_aus", driver = "GeoJSON")

plot_r <- function(x) {

  src_file <- paste0("../../", rasters$File[[x]])

  info <- attr(rgdal::GDALinfo(src_file, returnStats = FALSE), "df")

  gdtype <- as.character(info$GDType[[1]])

  r_crs <- if (tools::file_ext(src_file) == "vrt") {
    gdalUtils::gdalsrsinfo(src_file, as.CRS = TRUE)
  } else {
    crs(raster::raster(src_file))
  }

  r_crs <- if (is.na(r_crs)) "+init=epsg:4326" else as.character(r_crs)

  tmp <- tempfile()

  srcnodata <- if (gdtype == "Byte") 255 else c(-3.40282346638528e+38, -9999)
  if (x %in% 70:71)  srcnodata <- info$NoDataValue[[1]]

  gdalUtils::gdalwarp(
    src_file, tmp, r_crs, "+init=epsg:4326", ot = "Float32",
    srcnodata = srcnodata, dstnodata = -3.40282346638528e+38,
    cutline = grampians_json, crop_to_cutline = TRUE, dstalpha = TRUE,
    overwrite = TRUE
  )

  r_gr <- raster(tmp)

  tmp2 <- tempfile()

  gdalUtils::gdalwarp(
    src_file, tmp2, r_crs, "+init=epsg:4326", ot = "Float32",
    srcnodata = srcnodata, dstnodata = -3.40282346638528e+38,
    cutline = se_aus_json, crop_to_cutline = TRUE, dstalpha = TRUE,
    overwrite = TRUE
  )

  r_se <- raster(tmp2)

  par(mfrow = 1:2, oma = c(3, 3, 1, 1), mar = c(1, 1, 1, 1))
  plot(
    grampians,
    font.main = 1,
    panel.first = grid()
  )
  plot(
    r_gr, col = grey.colors(256, 1, 0), add = TRUE,
    horizontal = TRUE, legend.mar = 0
  )
  plot(
    se_aus,
    font.main = 1,
    panel.first = grid()
  )
  plot(
    r_se, col = grey.colors(256, 1, 0), add = TRUE,
    horizontal = TRUE, legend.mar = 0
  )
  title(rasters$Attribute[[x]], outer = TRUE)
}

for (i in seq(nrow(rasters))) plot_r(i)