need to plot multiple spatial maps side-by-side in a single window using sf plot function in R

Question

I like to plot all my 5 spatial maps of drought indices in a single window using sf plot function in R. I am able to plot all maps seperately but unable to plot altogether. All indexes are showing drought indices. I really appreciate your help!

###Agroindices
indices_det<-readRDS("D:/Ranjeet/IAMV6/output/agroindices/rice_grid_agindex_kharif_det.RDS")
View(indices_det)
cellID_yield_2000<-filter(indices_det, yearRef==2000)
view(cellID_yield_2000)
names(cellID_yield_2000)

#sum of drought indices for same GDMID using dplyr package

 drought_indices<-cellID_yield_2000%>%group_by(GDMID)%>%
 transmute(r_count_lt_p05_harvest_end=sum(r_count_lt_p05_harvest.end), 
        r_count_lt_p05_harvest = sum(r_count_lt_p05_harvest),
        r_count_lt_p05_main = sum(r_count_lt_p05_main),
        r_count_lt_p05_plant.start = sum(r_count_lt_p05_plant.start),
        r_count_lt_p05_plant = sum(r_count_lt_p05_plant))%>%unique()%>%
        arrange(GDMID, r_count_lt_p05_harvest_end, r_count_lt_p05_harvest,
                r_count_lt_p05_main, r_count_lt_p05_plant.start, r_count_lt_p05_plant)
view(drought_indices)


#India shape file for districts
library(sf)
shp_dist <- read_sf("D:/Ranjeet/IAMV6/input/shapefile/District_sim.shp")

#merging using merge function
merge_drought_indices <- merge(shp_dist, drought_indices)
view(merge_drought_indices)

#plot in a single window
library(sp)
par(mfrow=c(2,2))
 plot(merge_drought_indices["r_count_lt_p05_harvest_end"], 
 pal = colorRampPalette(c("red", "white", "blue","green")), 
 main = "End harvest index variation in India during drought year 2000",
 key.pos = 4, axes = TRUE, key.width = lcm(1.3), key.length = 1.0)
 
 plot(merge_drought_indices["r_count_lt_p05_harvest"], 
 pal = colorRampPalette(c("red", "white", "blue","green")), 
 main = "Harvest variation index in India during drought year 2000",
 key.pos = 4, axes = TRUE, key.width = lcm(1.3), key.length = 1.0)

 plot(merge_drought_indices["r_count_lt_p05_main"], 
 pal = colorRampPalette(c("red", "white", "blue","green")), 
 main = "Main growth stage index variation in India during drought year 2000",
 key.pos = 4, axes = TRUE, key.width = lcm(1.3), key.length = 1.0)

 plot(merge_drought_indices["r_count_lt_p05_plant.start"], 
 pal = colorRampPalette(c("red", "white", "blue","green")), 
 main = "Plant start index variation in India during drought year 2000",
 key.pos = 4, axes = TRUE, key.width = lcm(1.3), key.length = 1.0)

 plot(merge_drought_indices["r_count_lt_p05_plant"], 
 pal = colorRampPalette(c("red", "white", "blue","green")), 
 main = "main plant stage index variation in India during drought year 2000",
 key.pos = 4, axes = TRUE, key.width = lcm(1.3), key.length = 1.0)

Answer 1

One solution using ggplot (for lack of example data just the concept):

library(ggplot2)
library(dplyr)

merge_drought_indices %>%
  ## stack columns with drought indices rather than having them
  ## side by side:
  pivot_longer(cols = starts_with("r_count"),
               names_to = "index_type",
               values_to = "index_value"
               ) %>%
  ggplot() +
  geom_sf(aes(fill = index_value)) +
  ## since drought indices of various types are now stacked
  ## you can conveniently plot an array of panels-per-type:
  facet_wrap(~ index_type)
               
## save plot in desired format (deducted from file extension): 
ggsave("map_array.png", ...)

To adapt to your specific data, please read up on pivot_longer of {dplyr} and facet_wrap of {ggplot2}.