Week 12 practical: GIS

library(pacman)
p_load(tidyverse, sf, viridis, patchwork, RColorBrewer, rgdal, ggspatial, raster, GISTools, cowplot, here)

path <- glue::glue(here(), "/practicals/Data/Week12_Data/Week12_Data")

files <- list.files(path, "*", full.names = T)

agg.cov <-  read_csv(files[[1]], show_col_types = F)

head(agg.cov)

NA

p1 <- agg.cov %>%
  ggplot(aes(east, north, colour = SR)) +
  geom_point(shape = 15, size = 1) +
  coord_fixed() +
  scale_colour_viridis()

p1 <- p1 +
  theme(legend.title = element_blank(), 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.ticks = element_blank(),
        plot.background = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_blank(), 
        axis.text = element_blank(), 
        axis.title = element_blank()
        )
  
p1

p2 <-  agg.cov %>%
  ggplot(aes(east, north, colour = SR)) +
  geom_point(shape = 15, size = 1) +
  scale_colour_viridis()

p2

ireland <- st_read(files[[6]], quiet = T )

ireland_sf <- st_transform(st_as_sf(ireland), 29903)

p.future <- ireland_sf %>%
  ggplot() +
  geom_sf(fill = "white") +
  coord_sf()

p.future

load(files[[9]])
load(files[[10]])


change <- current.resistance %>%
  mutate(change = pred.resistance$pred.H45.50 - resist.fit) %>%
  dplyr::select(-resist.fit) %>%
  mutate(id = row_number()) %>%
  arrange(change)

o <- order(change$change)[c(1:26, 497:522)]


central <- change[-o, ]

extreme <- change[o,]

p.future <- p.future +
  theme_bw() +
  theme(axis.ticks = element_blank(),
        axis.text = element_blank(),
        panel.border = element_blank(),
        panel.grid = element_blank()) +
  labs(x = NULL, y = NULL)

p.future <- p.future +
  ggtitle("Change in resistance 2050")
  

pal <- brewer.pal(n = 11, name = 'PRGn')

p.future +
  geom_tile(data = central, aes(east, north, fill = change)) +
  geom_tile(data = extreme, aes(east, north, group = change)) +
  scale_fill_gradientn(limits = c(-0.0185, 0.0185), colours = pal) +
  theme(legend.title = element_text(size = 20), 
        legend.text = element_text(size = 20), 
        plot.title = element_text(size = 20), 
        plot.title.position = "plot") +
  annotation_north_arrow(location = 'tl', height = unit(1, "cm"), wirth = unit(1, "cm")) +
  annotation_scale(location = 'br', text_cex = 1.5)

plot_grid(p1, p2, nrow = 1, 
          labels = c('a)', 'b)'))

NA
NA

load(files[[10]])

head(stability)

NA
NA

df <- rasterFromXYZ(stability)

projection(df) <- CRS("+init=epsg:29903")

plot(df$evi.rec, box = FALSE,
     axes = FALSE, 
     main = "Recovery Time", 
     cex.main = 2,
     legend.width = 2, 
     axis.args = list(cex.axis = 2),
     col = viridis(15))
north.arrow(50000, 150000, len = 10000, lab ="N", col = "black")

scalebar(d = 100000,xy = c(200000, 50000), label = '100km')