Maps and plots for the eucalypt project

knitr::opts_chunk$set(cache = TRUE, dev = c("png", "pdf", "tiff"))
library(eucs)
library(ggplot2)
library(GGally)
library(here)
library(oz)
library(png)
library(raster)
library(rgeos)
library(sp)
library(tidyr)
library(xtable)

Plot Locations

ibra_subregions <- unique(
  c(
    subset(
      occupancy_southeast,
      occupancy_southeast$ibra_subregion_coverage > .5,
      ibra_subregion
    )[[1]],
    occupancy_grampians$ibra_subregion
  )
)

tmp <- tempfile()
png(tmp, 600, 500, bg = rgb(1, 1, 1, 1))
oz(FALSE)
plot(
  subset(map_ibra_southeast, IBRA_SUB_N %in% ibra_subregions),
  col = "grey",
  border = NA,
  add = TRUE
)
dev.off()

par(mar = c(3, 3, 1, 1))
plot(
  map_southeast, xlim = c(141, 151), ylim = -c(39, 34), border = "black",
  lwd = .5, panel.first = grid()
)
plot(
  subset(map_ibra_southeast, IBRA_SUB_N %in% ibra_subregions),
  col = grey.colors(length(ibra_subregions), .1, .8),
  border = NA,
  add = TRUE
)
points(
  occupancy_grampians[c("longitude_gda94", "latitude_gda94")], pch = 16,
  cex = .1, col = rgb(.1, .1, .1, .1)
)
points(
  occupancy_grampians[c("longitude_gda94", "latitude_gda94")], pch = 16,
  cex = .01, col = rgb(1, 1, 1)
)


text(142.4, -38, "Greater Grampians", pos = 3)
points(
  subset(
    occupancy_southeast,
    ibra_subregion_coverage > .5,
    c("longitude_gda94", "latitude_gda94")
  ),
  cex = .1,
  pch = 16,
  col = rgb(.1, .1, .1, .1)
)
points(
  subset(
    occupancy_southeast,
    ibra_subregion_coverage > .5,
    c("longitude_gda94", "latitude_gda94")
  ),
  cex = .01,
  pch = 16,
  col = rgb(1, 1, 1)
)

points(144.963056, -37.813611, pch = 16)
text(144.963056, -37.813611, label = "Melbourne", cex = .5, pos = 2)
points(151.209444, -33.865, pch = 16)
text(151.209444, -33.865, label = "Sydney", cex = .5, pos = 2)

axis(1, tick = FALSE)
axis(2, tick = FALSE, las = 1)
rasterImage(readPNG(tmp), 149, -39.2, 151, -37.9)

Correlations among sampled covariates

hist_data <- do.call(
  rbind,
  list(
    with(
      covariates_grampians,
      data.frame(
        Region = "Grampians",
        mlq = mean_moisture_index_of_low_qtr,
        twi = topographic_wetness_index_3sec,
        r1k = relief_1000m_radius,
        tn  = total_nitrogen
      )
    ),
    with(
      covariates_southeast,
      data.frame(
        Region = "South-east Australia",
        mlq = mean_moisture_index_of_low_qtr,
        twi = topographic_wetness_index_3sec,
        r1k = relief_1000m_radius,
        tn  = total_nitrogen
      )
    )
  )
)

hist_data$Region <- relevel(hist_data$Region, "South-east Australia")

par(mar = rep(.5, 4), oma = c(3, 4, 4, 3), mfrow = c(3, 3))
ramps <- list()
ramps[[1]] <- colorRampPalette(
  c(rgb(.7, .7, .7, 0), rgb(0, 0, 0, 1)), alpha = TRUE
)
ramps[[2]] <- colorRampPalette(
  c(rgb(.3, .3, .3, 0), rgb(1, 1, 1, 1)), alpha = TRUE
)
rgns <- c("South-east Australia" = FALSE, "Grampians" = TRUE)
vars <- list()
vars[[1]] <- c("mlq", "twi", "r1k")
vars[[2]] <- c("tn", "r1k", "twi")
labels <- c(
  mlq = "Mean Moisture Index\nLowest Quarter",
  twi = "Topographic\nWetness Index",
  r1k = "Relief\n(1000m Radius)",
  tn  = "Total Nitrogen"
)

for (i in seq(vars[[1]])) {
  for (j in seq(vars[[2]])) {
    if (sum(i, j) > 4) {
      plot.new()
      if (i == j) {
        legend(
          "center", c("Southeast","Grampians"),
          fill = c("black", "white"),
          bty = "n", cex = 2, xpd = NA
        )
      }
    } else {
      for (k in seq(rgns)) {
        smoothScatter(
          subset(
            hist_data, Region == names(rgns)[k], c(vars[[1]][j], vars[[2]][i])
          ),
          colramp = ramps[[k]], nrpoints = 0, add = rgns[k], xaxt = "n",
          yaxt = "n",
          panel.first = {
            rect(
              par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4],
              col = "lightgrey"
            )
            grid(col = "grey")
          },
          xlim = range(hist_data[[vars[[1]][j]]], na.rm = TRUE),
          ylim = range(hist_data[[vars[[2]][i]]], na.rm = TRUE)
        )
      }
      if (sum(i, j) == 4) {
        axis(1)
        axis(4, las = 1)
      }
      if (i == 1) mtext(labels[vars[[1]][j]], 3, 1)
      if (j == 1) mtext(labels[vars[[2]][i]], 2, 1)
    }
  }
}

Histograms of sampled covariates

hist_data$Region <- relevel(hist_data$Region, "Grampians")
hist_data <- gather(hist_data, type, value, -Region)
hist_data$type <- factor(hist_data$type, rev(levels(factor(hist_data$type))))
hist_data$Region <- factor(hist_data$Region, rev(levels(hist_data$Region)))

ggplot(na.omit(hist_data)) +
aes(value, fill = Region) +
geom_histogram(bins = 25) +
scale_fill_grey() +
facet_grid(
  ~type,
  labeller = as_labeller(labels),
  scales = "free_x"
) +
xlab("") +
ylab("Plots") +
theme_bw(base_size = 14) +
theme(
  legend.position = c(.975, .95), legend.justification = c("right", "top"),
  legend.title = element_blank()
)

rgns <- subset(map_ibra_southeast, IBRA_SUB_N %in% ibra_subregions)
km2 <- 1000000
aa <- CRS("+init=epsg:3577")

summary_table <- data.frame(
  Region = ibra_subregions,
  "Area (km$^2$)" = as.integer(round(as.vector(tapply(
    gArea(spTransform(rgns, aa), TRUE) / km2, rgns@data$IBRA_SUB_N, sum
  )[ibra_subregions]), -2)),
  "No. of plots" = as.vector(table(
    c(occupancy_grampians[, 2], occupancy_southeast[, 2])
  )[ibra_subregions]),
  "No. of taxa" = c(
    apply(
      with(
        subset(modeldata_southeast, occupancy),
        table(taxon, ibra_subregion) > 0
      ),
      2,
      sum
    ),
    `Greater Grampians` = length(unique(modeldata_grampians$taxon))
  )[ibra_subregions],
  check.names = FALSE,
  stringsAsFactors = FALSE
)

summary_table <- xtable(
  summary_table, caption = "IBRA subregions", label = "tab:summary"
)

dir.create(
  here("notebooks/eucalypt_project_maps_files/figure-html"),
  showWarnings = FALSE, recursive = TRUE
)

print(
  summary_table,
  file = here(
    "notebooks/eucalypt_project_maps_files/figure-html/summarytab.tex"
  ),
  include.rownames = FALSE,
  table.placement = "htbp!",
  caption.placement = "top",
  sanitize.text.function = identity
)