Subdivide tiles in ggplot2 geom_tile into quarters

I am trying to create a more elegant correlation plot for seasonal data using ggplot2 and geom_tile(). I want to split each tile in the plot into four quadrants (one quadrant for each seasonal value).

I already have code to plot each season's correlations separately:

library(dplyr)
library(ggplot2)

# simplified data for 1 season
dat.season1 <- data.frame(Var1 = c("species2", "species3", "species4", "species5",
                                            "species3", "species4", "species5",
                                            "species4", "species5", "species5"),
                          Var2 = c("species1", "species1", "species1", "species1",
                                            "species2", "species2", "species2",
                                            "species3", "species3", "species4"),
                          value = runif(10, -1, 1))

#plot Season 1
my.colors = colorRampPalette(c("#00002d", '#001a6d', '#99cce2','#ffffff', "#FED18A",'#fda416', "#ff6600"),
                             space = "rgb")

colorLevels<-9
cols_to_use= my.colors(colorLevels)

ggplot(data = dat.season1, mapping = aes_string(x = "Var2",y = "Var1", fill = "value")) +
  geom_tile(color = "gray60") +
  labs(x = NULL, y = NULL) +
  guides(size = "none")+
  scale_x_discrete(expand = c(0, 0), position = 'top') +
  scale_y_discrete(expand = c(0, 0), limits = rev) +
  scale_fill_gradientn(colors = cols_to_use,
                       space = "Lab",
                       guide = guide_colorbar(frame.colour = "gray60", frame.linewidth = 2))

This gives me a plot that looks like this: Example single season correlation plot

What I want is a plot that looks like this: example_subdivided_seasonal_plot

However, I have searched stackoverflow and googled, but can't figure out a way to subdivide each tile into quadrants based on a factor column.

Here is a sample of the combined seasonal dataset:

dat.season2 <- data.frame(Var1 = c("species2", "species3", "species4", "species5",
                                   "species3", "species4", "species5",
                                   "species4", "species5", "species5"),
                          Var2 = c("species1", "species1", "species1", "species1",
                                   "species2", "species2", "species2",
                                   "species3", "species3", "species4"),
                          value = runif(10, -1, 1))

dat.season3 <- data.frame(Var1 = c("species2", "species3", "species4", "species5",
                                   "species3", "species4", "species5",
                                   "species4", "species5", "species5"),
                          Var2 = c("species1", "species1", "species1", "species1",
                                   "species2", "species2", "species2",
                                   "species3", "species3", "species4"),
                          value = runif(10, -1, 1))

dat.season4 <- data.frame(Var1 = c("species2", "species3", "species4", "species5",
                                   "species3", "species4", "species5",
                                   "species4", "species5", "species5"),
                          Var2 = c("species1", "species1", "species1", "species1",
                                   "species2", "species2", "species2",
                                   "species3", "species3", "species4"),
                          value = runif(10, -1, 1))

all.dat <- dplyr::bind_rows(list(season1 = dat.season1, season2 = dat.season2,
                                  season3 = dat.season3, season4 = dat.season4),
                                  .id = "Season")

I found a couple threads that get close (here and here), but neither of them splits each tile into quadrants.

The second seems the most promising, giving a plot that looks like this.

But I can't figure out how to modify the code in the answer to give quadrants instead of vertical or horizontal bars.

One option would be to use facetting, i.e. facet by Var1 and Var2, map your seasons on "x" and "y" and use some styling to get rid of the axes and süacing between the panels:

library(ggplot2)
library(dplyr)

set.seed(123)

all.dat1 <- all.dat |>
  mutate(
    x = if_else(Season %in% paste0("season", c(1, 3)), 1, 2),
    y = if_else(Season %in% paste0("season", c(1, 2)), 1, 2),
    across(c(x, y), factor)
  )

season_labels <- all.dat1 |> 
  distinct(Season, x, y) |> 
  mutate(Var1 = "species2", Var2 = "species1")

ggplot(data = all.dat1, aes(x = x, y = y)) +
  geom_tile(aes(fill = value), color = "gray60") +
  geom_text(data = season_labels, aes(label = Season), color = "white") +
  facet_grid(Var1 ~ Var2, switch = "y") +
  labs(x = NULL, y = NULL) +
  guides(size = "none") +
  scale_x_discrete(expand = c(0, 0), position = "top") +
  scale_y_discrete(expand = c(0, 0), limits = rev) +
  scale_fill_gradientn(
    colors = cols_to_use,
    space = "Lab",
    guide = guide_colorbar(frame.colour = "gray60", frame.linewidth = 2)
  ) +
  theme(
    axis.text = element_blank(), 
    axis.ticks = element_blank(), 
    axis.ticks.length = unit(0, "pt"),
    panel.spacing = unit(0, "pt"), 
    strip.background = element_rect(color = "black")
  )