Join multiple rasters at specific locations, when the number of rasters is variable

Given "n" rasters, each exactly 100x100 px, I want to combine them all in a single plot such that there are exactly 4 images per row. For n > 4, a new row should be created and so on. Thus, the dimensions of the empty plot created will be fixed along the x-axis, but the y-axis will depend on the number of rasters. I used image_montage() function from "magick" package to generate montages as follows:

mag_montage <- list()
for(ii in 1:n){
filelist_crop <- list.files()[grep(".png",list.files())]  
mag_montagetemp <- image_read(filelist_crop)
mag_montage[[ii]] <- image_montage(mag_montagetemp)
}

But I can't control the specific position of each individual file within the file montage this way. Knowing the locations is very crucial, as I need to pick certain xy coordinates (using "locator") from the combined rasters for some downstream processing. Any help will be much appreciated. Thanks.

This way you can plot all the rasters inside the list rlist into one plot with 4 columns:

library(raster)

n <- 26

rlist <- lapply(1:n,function(x) raster(system.file("external/test.grd", package="raster")))

par(mfrow=c(ceiling(n/4),4))

for (ii in 1:length(rlist)){

  plot(rlist[[ii]])

## additional options for plot to omit legend and box
#bty="n", box=FALSE, axes=F, legend=F

}