I want to multiply Landsat 30m resolution image with an ERA5 raster of about 25km resolution. I can downscale ERA5 to 30m but it would increase processing. I am just wondering is it possible that whichever 30m cells overlap with 25km cell, these 30m cells get multiplied by single value from 25km cell.
library(terra)
### 1st raster
ETrF <- rast(ncols=4, nrows=4, xmin=73, xmax=75, ymin=31, ymax=33)
### 2nd raster
ET0 <- rast(ncols=2, nrows=2, xmin=73, xmax=75, ymin=31, ymax=33)
### Random values assigned to both rasters
values(ETrF) <- c(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8)
values(ET0) <- c(5,1,2,4)
### Now 1st cell of ET0 (value: 5) will overlap with 1st,2nd,5th and 6th cells (values: 0.1,0.2,0.5,0.6)
### of ETrF.Is it possible that 5 gets multiplied with 0.1,0.2,0.5 and 0.6 and final output
### keeps ETrFs resolution.
ETa <- ETrF*ET0.....?
A reason for not implementing that is that you need to chose how to resample the data; x * y cannot capture that.
In your case, the default bilinear interpolation is probably the most reasonable option when using resample (unless you need a more sophisticated methods that use elevation and other covariables).
But it does not give the same results as with the extract method. (you need to use resample( , method="near") for that). With the toy data, you can also use disagg.
If you want to use the extract route, that can be made much more efficient by avoiding the coercion to sf, at least for this toy example:
fastextr <- function(coarse, fine, method="simple") {
xy <- crds(fine)
value <- extract(coarse, xy, method=method) * values(fine)
rasterize(xy, fine, value)
}
microbenchmark::microbenchmark(
"res" = terra::resample(ET0, ETrF, "near", threads = TRUE) * ETrF,
"mutat" = extract_mutate_rasterize(ET0, ETrF),
"disagg" = terra::disagg(ET0, 2) * ETrF,
"fextr" = fastextr(ET0, ETrF),
times = 50
)
#Unit: milliseconds
# expr min lq mean median uq max neval
# res 6.2586 6.4362 6.948614 6.67490 6.8020 21.1767 50
# mutat 38.6304 39.4301 40.627638 40.08070 41.1208 58.3876 50
# disagg 4.5824 4.6490 5.158002 4.77420 4.9126 21.7821 50
# fextr 7.9329 8.1154 8.468164 8.44225 8.6150 11.1409 50
So your best option is to use bilinear resampling
bres <- terra::resample(ET0, ETrF, "bilinear", threads = TRUE) * ETrF
The same result can be obtained with bilinear extract
bextr <- fastextr(ET0, ETrF, "bilinear")
This extract approach is not much slower than the resample approach for this example. But I would not use the extract method because it may cause problems with memory allocation if you are using large raster data.