ggplot spatial data and ggsave unique id's in a for loop

I have a dataset with watersheds across 3 provinces. Each watershed has a corresponding shapefile, and score from 0 to 1 for ecosystem disruption under 3 rcp (climate projection models). Here is the head of my dataframe:

ws            rcp  ecodisrup                       geometry
ANNAPOLIS_sp rcp26 0.09090909 MULTIPOLYGON (((-64.86239 4...
BARRINGTON_C rcp26 0.00000000 MULTIPOLYGON (((-65.4722 43...
CHETICAMP_RI rcp26 0.09090909 MULTIPOLYGON (((-60.59559 4...
CLAM_HRB_ST. rcp26 0.27272727 MULTIPOLYGON (((-61.38909 4...
COUNTRY_HARB rcp26 0.09090909 MULTIPOLYGON (((-61.96444 4...
EAST_INDIAN_ rcp26 0.09090909 MULTIPOLYGON (((-63.94016 4...

The watersheds and their shapefiles repeat for each rcp (26, 45, 85), eg:

ws            rcp  ecodisrup                       geometry
ANNAPOLIS_sp rcp26 0.09090909 MULTIPOLYGON (((-64.86239 4...
...          rcp26 0.00000000 MULTIPOLYGON (((-65.4722 43...
46th ws      rcp26 0.09090909 MULTIPOLYGON (((-60.59559 4...
ANNAPOLIS_sp rcp45 0.27272727 MULTIPOLYGON (((-64.86239 4...
...          rcp45 0.09090909 MULTIPOLYGON (((-65.4722 43...
46th ws      rcp45 0.09090909 MULTIPOLYGON (((-60.59559 4...
ANNAPOLIS_sp rcp85 0.09090909 MULTIPOLYGON (((-64.86239 4...
...          rcp85 0.00000000 MULTIPOLYGON (((-65.4722 43...
46th ws      rcp85 0.09090909 MULTIPOLYGON (((-60.59559 4...

I'd like to map the watersheds, coloured by ecodisrup, under each of the three rcp scenarios by writing a for loop. The problem I run into is including a ggsave function in the for loop, and also the loop I'm writing is just producing the same map three times.

Later on in my analysis I'll need the same type of for loop for mapping, but I'll add in a species column (produce a map coloured by ecodisrup for each species across all watersheds, under 3 rcp scenarios).

Here is the code I've tried before, expected an output of three maps (one for each rcp). Instead I got three maps, but they all had the exact same colouring (not sure if they just coloured by the last one? rcp85?).

# make empty list to fill in with plots
ecodis_rcp_plots = list()

# define the different rcps
ecodis_rcp = unique(ecodis$rcp)

# begin for loop
for (rcp in ecodis_rcp){
  ecodis_rcp_plots[[rcp]] = ggplot(data = ecodis, aes(geometry = geometry)) + 
    geom_sf(aes(fill = ecodisrup)) +
    scale_fill_viridis_c(option = "viridis", limits = c(0, 1))
  ggsave(paste("C:/Users/myname/Desktop/EcoDis", rcp, ".png"),ecodis_rcp_plots[[rcp]],width=8,height=8,units="in",dpi=300)
}

Any insight would be greatly appreciated, thank you!

Solution

One option is to split the dataset first and then cycle through the list of resulting sf objects. Following is based on a nc dataset from sf package for a more generic example:

library(ggplot2)
library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

# prepare some example data
nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
nc$sample_group <- paste("Group", rep(1:5, each = 20))
nc$area <- scales::rescale(nc$AREA)
nc[,c("NAME", "sample_group", "area")]
#> Simple feature collection with 100 features and 3 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27
#> First 10 features:
#>           NAME sample_group       area                       geometry
#> 1         Ashe      Group 1 0.36180905 MULTIPOLYGON (((-81.47276 3...
#> 2    Alleghany      Group 1 0.09547739 MULTIPOLYGON (((-81.23989 3...
#> 3        Surry      Group 1 0.50753769 MULTIPOLYGON (((-80.45634 3...
#> 4    Currituck      Group 1 0.14070352 MULTIPOLYGON (((-76.00897 3...
#> 5  Northampton      Group 1 0.55778894 MULTIPOLYGON (((-77.21767 3...
#> 6     Hertford      Group 1 0.27638191 MULTIPOLYGON (((-76.74506 3...
#> 7       Camden      Group 1 0.10050251 MULTIPOLYGON (((-76.00897 3...
#> 8        Gates      Group 1 0.24623116 MULTIPOLYGON (((-76.56251 3...
#> 9       Warren      Group 1 0.38190955 MULTIPOLYGON (((-78.30876 3...
#> 10      Stokes      Group 1 0.41206030 MULTIPOLYGON (((-80.02567 3...

# split by some factor, resulting list is named according to factor levels:
nc_split <- split(nc, ~sample_group)

# create named list of nc_split names for lapply(), this helps us to get 
# a named list from lapply() while providing access to sample_group values
nc_groups <- setNames(names(nc_split), names(nc_split))
nc_groups
#>   Group 1   Group 2   Group 3   Group 4   Group 5 
#> "Group 1" "Group 2" "Group 3" "Group 4" "Group 5"

# cycle though the list of sf objects and generate list of plots
plots <- lapply(nc_groups, \(grp_name){
  ggplot(nc_split[[grp_name]], aes(fill = area)) +
    geom_sf() +
    scale_fill_viridis_c(limits = c(0, 1)) +
    ggtitle(grp_name) +
    theme_void() +
    theme(legend.position = "none")
})
# resulting list of plots, named:
str(plots, max.level = 1)
#> List of 5
#>  $ Group 1:List of 9
#>   ..- attr(*, "class")= chr [1:2] "gg" "ggplot"
#>  $ Group 2:List of 9
#>   ..- attr(*, "class")= chr [1:2] "gg" "ggplot"
#>  $ Group 3:List of 9
#>   ..- attr(*, "class")= chr [1:2] "gg" "ggplot"
#>  $ Group 4:List of 9
#>   ..- attr(*, "class")= chr [1:2] "gg" "ggplot"
#>  $ Group 5:List of 9
#>   ..- attr(*, "class")= chr [1:2] "gg" "ggplot"

# visualise all plots:
patchwork::wrap_plots(plots, nrow = 3)


# save all
for (plot_name in names(plots)) {
  ggsave(paste0("nc_plot ", plot_name, ".png"),plots[[plot_name]])
}
#> Saving 7 x 5 in image
#> Saving 7 x 5 in image
#> Saving 7 x 5 in image
#> Saving 7 x 5 in image
#> Saving 7 x 5 in image

# check resulting files
fs::dir_info(glob = "nc_plot*")[,1:3]
#> # A tibble: 5 × 3
#>   path                type         size
#>   <fs::path>          <fct> <fs::bytes>
#> 1 nc_plot Group 1.png file        84.1K
#> 2 nc_plot Group 2.png file        96.7K
#> 3 nc_plot Group 3.png file        96.3K
#> 4 nc_plot Group 4.png file        88.9K
#> 5 nc_plot Group 5.png file        93.9K

^{Created on 2023-06-16 with reprex v2.0.2}