Geospatial mapping using ggplot faceting in R?

Question

I have a Watershed that has 33 subbasins. I wasn't able to come up with a code to draw a reproducible shapefile for the catchment so i am attaching my Shapefile. I have four models that generate Evapotranspiration (ET) data for the years 2005-2008. I would like to compare the four models products of each year using ggplot faceting functionality. I tried a few things (see my sample code) but couldn't succeeded. I would appreciate a way forward.

library(sf)
library(tidyverse)

shape <- read_sf(dsn = ".", layer = "Watershed")


ETM1 = data.frame(Subbasin = 1:33, Yr2005 = runif(33, 500,700), Yr2006 = runif(33, 600,750), Yr2007 = runif(33, 450,750),
                  Yr2008 = runif(33, 550,800), Model = rep("M1", 33))

ETM2 = data.frame(Subbasin = 1:33, Yr2005 = runif(33, 600,750), Yr2006 = runif(33, 550,750), Yr2007 = runif(33, 600,750),
                  Yr2008 = runif(33, 700,800), Model = rep("M2", 33))

ETM3 = data.frame(Subbasin = 1:33, Yr2005 = runif(33, 500,750), Yr2006 = runif(33, 650,750), Yr2007 = runif(33, 700,750),
                  Yr2008 = runif(33, 500,800), Model = rep("M3", 33))

ETM4 = data.frame(Subbasin = 1:33, Yr2005 = runif(33, 400,750), Yr2006 = runif(33, 450,750), Yr2007 = runif(33, 300,750),
                  Yr2008 = runif(33, 400,800), Model = rep("M4", 33))

ETData = rbind(ETM1,ETM2,ETM3,ETM4)
Combine = gather(ETData, key = "Variable", value = "Value", -c("Model","Subbasin"))

SpData = merge(shape, Combine, by='Subbasin')

ggplot() + 
  geom_polygon(SpData, aes(x = Lat, y = Long_, fill = Value))+
  facet_wrap(~Model, nrow = 4, ncol = 4)

here is a picture of my shapefile that i draw using plot(shape$geometry)

Here is an approximate figure that i would like to construct, though it has only 3 rows.

I made a figure (hand drawn) that reflects my ultimate goal- each oval shape (sort of) represent my Catchment shapefile with division as subbasins. I apologies for my bad drawing.

Answer 1

You were quite close to the solution. Your final data contains everything to get your plot.

Basically, you can use facet_grid to separate both "Year" and "Model" variables and geom_sf that will use the column "geometry" of your dataset to plot the shape (more informations here: https://ggplot2.tidyverse.org/reference/ggsf.html). Passing the value as a filling into the aes will make the rest.

Here, I used scale_fill_gradient to set the color and few functions to get the plot as close as your desired plot:

library(sf)
library(tidyverse)

ggplot() +
  geom_sf(data = SpData, aes(fill = Value))+
  facet_grid(Variable~Model)+
  scale_fill_gradient(name =  "Evapotranspiration (ET)", low = "green", high = "red",
                      limits = c(300,900), 
                      breaks = c(300, 500, 700, 900))+
  theme_void()+
  theme(legend.position = "bottom")