Maps of model covariates

knitr::opts_chunk$set(cache = TRUE, dev = c("png", "pdf"))
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")

nms <- c(
  `62` = "Mean Moisture Index\nLowest Quarter",
  `24` = "Topographic Wetness\nIndex",
  `22` = "Relief 1000m\nRadius",
  `8`  = "Total Nitrogen"
)

zlims <- list(
  `62` = c(0.09091329, 1),
  `24` =  c(0, 19.50653),
  `22` = c(0, 1110.44),
  `8`  = c(0, 0.1795395)
)

par(mfrow = c(4, 2), oma = c(1, 1, 1, 1))
for (x in c(62, 24, 22, 8)) {
  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)

  zlim <- zlims[[as.character(x)]]
  cols <- grey.colors(256, 1, 0)

  par(mar = c(0, .5, 0, 0))
  plot(
    grampians,
    border = rgb(1, 1, 1, 0),
    panel.first = grid()
  )
  plot(
    r_gr, col=cols, add = TRUE,
    legend = FALSE, zlim = zlim
  )
  if (x == 8) axis(1, tick = FALSE, line = -1, hadj = 0)
  axis(2, las = 1, tick = FALSE, padj = 0, line = -2)
  text(141.6, -36.9, nms[[as.character(x)]], pos = 4)
  par(mar = c(0, .5, 0, 1))
  plot(
    se_aus,
    border = rgb(1, 1, 1, 0),
    panel.first = grid()
  )
  plot(
    r_se, col=cols, add = TRUE, zlim = zlim,
    legend.shrink = .87
  )
  axis(2, las = 1, tick = FALSE, padj = 0, line = -2)
  if (x == 8) axis(1, tick = FALSE, line = -1, hadj = 0)
}