a = raster("~/a_rstr.tif")

b = raster("~/b_rstr.tif")

# Stack covariates
a_b_stack = stack(a,b)
names(a_b_stack) <- c("a", "b")

# Use the function on two rasters at the same time. We define a third raster (a_b_nb) that records the positions of values from 1 to ncell. focal will thus extract the positions for which we can get values of the two rasters (as they now have the same extent, origin and resolution).

a_b_nb = raster(a_b_stack, 1)
values(a_b_nb  = 1:ncell(a_b_stack)

matrix_a_b= values(a_b_stack) # stack as raster


# Define a function

focal_cor = focal(
  x = a_b_nb,
  w = matrix(1, 5, 5),
  fun = function(x, y = matrix_a_b){ # [MW]
    cor(y[x, 1], y[x, 2], # [MW]
        use = "na.or.complete")
  },
  filename = file.path("~/focal_cor.tif"),
  overwrite = TRUE
)

focal_cor = raster("~focal_cor.tif")
plot(focal_cor)

# Get data for ggplot
focal_cor_dat = SDMSelect::gplot_data(focal_cor, maxpixels = 50000) %>%
  mutate(variable = "Correlation") %>%
  filter(!is.na(value))

# Plot
ggplot(focal_cor_dat) +
  geom_tile(aes(x, y, fill = value)) +
  scale_fill_gradient2("Corr",
    low = "#d7191c",
    mid = "#ffffbf",
    high = "#1a9641",
    midpoint = 0) +
  # geom_polygon(data = Europe_border,
  #              aes(long, lat, group = group),
  #              fill = "white",
  #              colour = "grey20", size = 0.1) +
  coord_quickmap(
    xlim = range(focal_cor_dat$x),
    ylim = range(focal_cor_dat$y)
  ) + xlab("") + ylab("")