R: Raster Mosaic from List of Rasters

Mosaicking rasters from the same date in a list of rasters

You can do it inside a loop. Using a reproducible example:

set.seed(123)

rlist <- list() # store all rasters

rlist[[1]] <- raster(nrows=90, ncols=180, xmn=-180, xmx=0, ymn=-90, ymx=90,vals=rnorm(16200))
rlist[[2]] <- raster(nrows=90, ncols=180, xmn=0, xmx=180, ymn=-90, ymx=90,vals=rnorm(16200))
rlist[[3]] <- raster(nrows=90, ncols=180, xmn=-180, xmx=0, ymn=-90, ymx=90,vals=rnorm(16200))
rlist[[4]] <- raster(nrows=90, ncols=180, xmn=0, xmx=180, ymn=-90, ymx=90,vals=rnorm(16200))

dates_2 <- c("20160908","20160918","20160918","20160925") # you dates

dates_2u <- unique(dates_2) # create a vector with unique dates

rlist2 <- list() # a second list to store processed rasters

for(i in seq_along(dates_2u)){
  idx <- which(dates_2 %in% dates_2u[i]) #some index
  if(length(idx) > 1){
    rlisttemp <- rlist[idx] # to create a temporal list
    rlisttemp$fun <- mean # set function to make mosaic
    rlist2[[i]] <- do.call(mosaic,rlisttemp)
  }else{
    rlist2[[i]] <- rlist[[idx]] # for non-mosaicking images
  }
}

And check results:

plot(rlist2[[1]]);plot(rlist2[[2]]);plot(rlist2[[3]])

In you case, you can load all raster inside a list (to maintain example's code) with something like:

rlist <- lapply(list.files('/your/path/to/rasters',
                           pattern = '.tif$',full.names = T,
                           recursive = T),FUN=raster)

Or:

rlist <- lapply(lst_B5, FUN=raster)

Work with large raster mosaics in R without merging them to a single file (like lidR catalog)

You should look into the terra package which provides exactly the functionality you're looking for through virtual raster tiles (VRTs). We can use them to treat a collection of raster files on disk as a single raster file while taking advantage of the API to do a majority of the same tasks as you can do through the raster package.

First, let's create a sample of 4 rasters using the example straight from the ?terra::vrt() documentation.

library(terra)

r <- rast(ncols=100, nrows=100)
values(r) <- 1:ncell(r)
x <- rast(ncols=2, nrows=2)
filename <- paste0(tempfile(), "_.tif")
ff <- makeTiles(r, x, filename)
ff
#> [1] "/var/folders/b7/_6hwb39d43l71kpy59b_clhr0000gn/T//RtmpACJYNv/filedf6b65d4fca4_1.tif"
#> [2] "/var/folders/b7/_6hwb39d43l71kpy59b_clhr0000gn/T//RtmpACJYNv/filedf6b65d4fca4_2.tif"
#> [3] "/var/folders/b7/_6hwb39d43l71kpy59b_clhr0000gn/T//RtmpACJYNv/filedf6b65d4fca4_3.tif"
#> [4] "/var/folders/b7/_6hwb39d43l71kpy59b_clhr0000gn/T//RtmpACJYNv/filedf6b65d4fca4_4.tif"

Now, we'll read them in as a VRT, again, straight from the same example. This allows

vrtfile <- paste0(tempfile(), ".vrt")
v <- vrt(ff, vrtfile)
head(readLines(vrtfile))
#> [1] "<VRTDataset rasterXSize=\"100\" rasterYSize=\"100\">"                                                                                                                                                                                                                                                                                                             
#> [2] "  <SRS dataAxisToSRSAxisMapping=\"2,1\">GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]</SRS>"
#> [3] "  <GeoTransform> -1.8000000000000000e+02,  3.6000000000000001e+00,  0.0000000000000000e+00,  9.0000000000000000e+01,  0.0000000000000000e+00, -1.8000000000000000e+00</GeoTransform>"                                                                                                                                                                             
#> [4] "  <VRTRasterBand dataType=\"Float32\" band=\"1\">"                                                                                                                                                                                                                                                                                                                
#> [5] "    <NoDataValue>nan</NoDataValue>"                                                                                                                                                                                                                                                                                                                               
#> [6] "    <ColorInterp>Gray</ColorInterp>"
v
#> class       : SpatRaster 
#> dimensions  : 100, 100, 1  (nrow, ncol, nlyr)
#> resolution  : 3.6, 1.8  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : filedf6b216a737.vrt 
#> name        : filedf6b216a737 
#> min value   :               1 
#> max value   :           10000

Now, we can just create a simple example polygon to restrict our region of interest from -180 to 180 to -90 to 90 longitude.

library(sf)

pl <- list(rbind(c(-90,-90), c(-90,90), c(90,90), c(90,-90), c(-90,-90)))
roi <- st_sfc(st_polygon(pl), crs = "EPSG:4326")

crop(v, roi)
#> class       : SpatRaster 
#> dimensions  : 100, 50, 1  (nrow, ncol, nlyr)
#> resolution  : 3.6, 1.8  (x, y)
#> extent      : -90, 90, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : memory 
#> name        : filedf6b216a737 
#> min value   :              26 
#> max value   :            9975

Reduce memory usage for mosaic on large list of rasters

Rather than loading all rasters into memory at once (in the mosaic() call), can you process them one at a time? That way, you have your mosaic that updates each time you bring one more raster into memory, but then you can get rid of the new raster and just keep the continuously updating mosaic raster.

Assuming that your rlist object is a list of rasters, I'm thinking of something like:

Pseudocode

1. Initialize an updating_raster object as the first raster in the list
2. Loop through each raster in the list in turn, starting from the 2nd raster
3. Read the ith raster into memory called next_raster
4. Update the updating_raster object by overwriting it with the mosaic of itself and the next raster using a weighted mean

R code

Testing with the code in the mosaic() help file example...

First generate some rasters and use the standard mosaic method.

library(raster)

r <- raster(ncol=100, nrow=100)
r1 <- crop(r, extent(-10, 11, -10, 11))
r2 <- crop(r, extent(0, 20, 0, 20))
r3 <- crop(r, extent(9, 30, 9, 30))

r1[] <- 1:ncell(r1)
r2[] <- 1:ncell(r2)
r3[] <- 1:ncell(r3)

m1 <- mosaic(r1, r2, r3, fun=mean)

Put the rasters in a list so they are in a similar format as I think you have.

rlist <- list(r1, r2, r3)

Because of the NA handling of the weighted.mean() function, I opted to create the same effect by breaking down the summation and the division into distinct steps...

First initialize the summation raster:

updating_sum_raster <- rlist[[1]]

Then initialize the "counter" raster. This will represent the number of rasters that went into mosaicking at each pixel. It starts as a 1 in all cells that aren't NA. It should properly handle NAs such that it only will increment for a given pixel if a non-NA value was added to the updating sum.

updating_counter_raster <- updating_sum_raster
updating_counter_raster[!is.na(updating_counter_raster)] <- 1

Here's the loop that doesn't require all rasters to be in memory at once. The counter raster for the raster being added to the mosaic has a value of 1 only in the cells that aren't NA. The counter is updated by summing the current counter raster and the updating counter raster. The total sum is updated by summing the current raster values and the updating raster values.

for (i in 2:length(rlist)) {

  next_sum_raster <- rlist[[i]]
  next_counter_raster <- next_sum_raster
  next_counter_raster[!is.na(next_counter_raster)] <- 1

  updating_sum_raster <- mosaic(x = updating_sum_raster, y = next_sum_raster, fun = sum)
  updating_counter_raster <- mosaic(updating_counter_raster, next_counter_raster, fun = sum)

}

m2 <- updating_sum_raster / updating_counter_raster

The values here seem to match the use of the mosaic() function

identical(values(m1), values(m2))
> TRUE

But the rasters themselves aren't identical:

identical(m1, m2)
> FALSE

Not totally sure why, but maybe this gets you closer?

Perhaps compareRaster() is a better way to check:

compareRaster(m1, m2)
> TRUE

Hooray!

Here's a plot!

plot(m1)
text(m1, digits = 2)
plot(m2)
text(m2, digits = 2)

A bit more digging in the weeds...

From the mosaic.R file:

It looks like the mosaic() function initializes a matrix called v to populate with the values from all the cells in all the rasters in the list. The number of rows in matrix v is the number of cells in the output raster (based on the full mosaicked extent and resolution), and the number of columns is the number of rasters to be mosaicked (11,000) in your case. Maybe you're running into the limits of matrix creation in R?

With a 1000 x 1000 raster (1e6 pixels), the v matrix of NAs takes up 41 GB. How big do you expect your final mosaicked raster to be?

r <- raster(ncol=1e3, nrow=1e3)
x <- 11000
v <- matrix(NA, nrow=ncell(r), ncol=x)
format(object.size(v), units = "GB")
[1] "41 Gb"

Error in mosaic of rasters from different extent using terra package in R?

The error was because the rasters were not aligned, and due to a bug. I now get

library(terra)
#terra version 1.2.1
crs <- "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m"
r1 <- rast(nrow=1667, ncol=1667, ext=c(-1800000, -1749990, -600010, -550000), crs=crs)
r2 <- rast(nrow=1667, ncol=1667, ext=c(-1800000, -1749990, -550010, -5e+05), crs=crs)
values(r1) <- 1:ncell(r1)
values(r2) <- 1:ncell(r2)
m <- mosaic(r1, r2)
#Warning message:
#[mosaic] rasters did not align and were resampled

m
#class       : SpatRaster 
#dimensions  : 3334, 1667, 1  (nrow, ncol, nlyr)
#resolution  : 30, 30  (x, y)
#extent      : -1800000, -1749990, -600010, -499990  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#name        :   lyr.1 
#min value   :       1 
#max value   : 2778889 

Also

@Elia's suggestion is a good work-around:

r1 <- writeRaster(r1, "test1.tif", overwrite=TRUE)
r2 <- writeRaster(r2, "test2.tif", overwrite=TRUE)
v <- vrt(c("test1.tif", "test2.tif"), "test.vrt", overwrite=TRUE)

v
#class       : SpatRaster 
#dimensions  : 3334, 1667, 1  (nrow, ncol, nlyr)
#resolution  : 30, 30  (x, y)
#extent      : -1800000, -1749990, -600020, -5e+05  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source      : test.vrt 
#name        :    test 
#min value   :       1 
#max value   : 2778889 

Why is Mosaic rasters not working in some files?

The raster stacks need to have the same number of layers

nlayers(test_1)
# [1] 4
nlayers(test_2)
# [1] 3

If we add another dummy layer filled with NA, then the mosaic works. I assumed that the missing layer is the fourth, but you may need to figure out which one is missing in your specific case.

test_2a = stack(test_2, init(test_2[[1]], fun=function(x) rep(NA,x)))
joint <- mosaic(test_1, test_2a, fun=mean, filename = "joint.tif", overwrite=TRUE)

Or, we can use package terra instead of raster

library(terra)
test_1 <- rast(files[1])
test_2 <- rast(files[2])
joint <- mosaic(test_1, test_2)

How to extract rasters with the same extent from a list of rasters in R

You can use sapply:

same_as_r1 <- sapply(list_raster, function(x) extent(x) == extent(list_raster$raster1))

And subset your list with it:

group1 <- list_raster[same_as_r1]
group2 <- list_raster[!same_as_r1]

---

Related Topics