Create a gradient mosaic from two raster images in R

Question

I have two different raster layers with different masked areas. They overlap by a few pixels, which is desired. So i want to merge both layers with the overlapping values "smoothly" transitioning. However, with the terra package best i can get is fun = "mean". Ive constructed an example, which shows what i did:

library(terra)

r <- rast(ncols = 36, nrows = 18)
r[450] <- 1
r <- 1 * buffer(r, width = 5000000) 
r[r == 0] <- NA
plot(r)

r

b <- rast(ncols = 36, nrows = 18)
b[345] <- 1
b <- 1 * buffer(b, width = 5000000)
b[b == 1] <- 5
b[b == 0] <- NA
plot(b)

b

m <- mosaic(r, b, fun = "mean")
plot(m)

m

The problem now is, that i want the transition between r and b to be smooth and not just mean. So there should be a transition from e.g. 2, 3, 4 from r to b in the overlapping pixels. Can anyone help me here? I haven't found any similar problem online, to be honest, i don't even know what to search for in this case.

Thanks in regards.

Answer 1

library(cgwtools)
library(terra)

rb_intersect <- intersect(r, b)
rb_intersect
class       : SpatRaster 
dimensions  : 18, 36, 1  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
source(s)   : memory
name        : layer 
min value   : FALSE 
max value   :  TRUE 

which(rb_intersect[r] == TRUE)
 [1]  2  3  4  5  8  9 10 11 12 16 17 18 19 20 21 26 27 28 29 30 31 32 38 39 40
[26] 41 42 43 50 51 52 53 54

These series look 'rowish', two rle-like views, cgwtools::seqle gives rle of sequences

rle(which(rb_intersect[r] == TRUE))
Run Length Encoding
  lengths: int [1:33] 1 1 1 1 1 1 1 1 1 1 ...
  values : int [1:33] 2 3 4 5 8 9 10 11 12 16 ...

cgwtools::seqle(which(rb_intersect[r] == TRUE))
Run Length Encoding
  lengths: int [1:6] 4 5 6 7 6 5
  values : int [1:6] 2 8 16 26 38 50

The lengths are of interest as they suggest the allocation of 2,3,4 transitioning to 5, but is the logic that it should happen up or across? Given that 50:55 is both downmost and closest to 1 values, perhaps these are your 2(s). How 2,3,4,5 are allocated, up or across, more 2s than 3s, might be indicated by your data.

Consider that this might be a completely foolish approach, though one I'd likely employ myself, and others may provide better guidance. But, continuing on into the fog

rb_mtx = as.matrix(rb_intersect, wide = TRUE) # TRUE/FALSE/NA
rb_mtx_01 = +rb_mtx # convert T/F/NA to 1/0/NA
rb_01_idx = which(rb_mtx_01 == TRUE, arr.ind = TRUE)
rb_01_idx_ord = rb_01_idx[order(rb_01_idx[,1]), ]
head(rb_01_idx_ord)
     row col
[1,]   9  17
[2,]   9  18
[3,]   9  19
[4,]   9  20
[5,]  10  17
[6,]  10  18

Now have something to reason about relative a gradient, and looked at in context of plotting r, b, and rb_intersect we see defensible decisions need making where rb_mtx_01[9, 17] <- could be 3, but [9, 18:20] <- should be 5. Make your decisions for these TRUE(s). Make another index for the FALSE(s) assigning 1 or 5. rb_gradient <- rast(rb_mtx_01), and you've got a gradient that hopefully reflects reality. With your FALSE index

seqle(rb_01_idx_F_ord[which(rb_01_idx_F_ord == 9), 2])
Run Length Encoding
  lengths: int [1:2] 1 5
  values : int [1:2] 16 21 # col 16 = 1, col 21:(21+5-1) = 5

rb_mtx_01[9, 16] <- 1
rb_mtx_01[9, 21:25] <- 5

# can use `length` to test which rows are going to require this hand work
length(seqle(rb_01_idx_F_ord[which(rb_01_idx_F_ord == 9), 2])$lengths)
[1] 2
length(seqle(rb_01_idx_F_ord[which(rb_01_idx_F_ord == 6), 2])$lengths)
[1] 1

As mentioned, plenty of hand work, even in this small example. It is useful, none the less, when producing counter-factual data where (returning to elevation), the satellite data says the river goes here, the river got 'moved' 500 years ago and you're trying to model the river 800 years ago. Hand work.