如何从 1 个具有不同范围和分辨率的栅格中减去 3 个栅格

我有 4 个具有不同分辨率和范围的栅格。谁能帮助我如何从栅格（d1）中减去3个栅格（a，b，c）以获得名为“e”的新输出栅格

就像 e= d1-a-b-c。 `

d1
class      : RasterLayer 
dimensions : 180, 360, 64800  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : -1.110223e-16, 360, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : Liquid_Water_Equivalent_Thickness 
values     : -249.2061, 806.3248  (min, max)

> a
class      : RasterLayer 
dimensions : 39, 46, 1794  (nrow, ncol, ncell)
resolution : 0.826087, 1.153846  (x, y)
extent     : 48, 86, 6, 51  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 0, 0.4987984  (min, max)

> b
class      : RasterLayer 
dimensions : 39, 46, 1794  (nrow, ncol, ncell)
resolution : 0.826087, 1.153846  (x, y)
extent     : 48, 86, 6, 51  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 0, 555.5283  (min, max)

> c
class      : RasterLayer 
band       : 1  (of  4  bands)
dimensions : 46, 39, 1794  (nrow, ncol, ncell)
resolution : 1, 1  (x, y)
extent     : 47.5, 86.5, 5.5, 51.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
source     : D:/GLDAS Data/2002_GLDAS/GLDAS_NOAH10_M.A200201.001.grb.SUB.nc4 
names      : Average.layer.1.soil.moisture 
z-value    : 2 
zvar       : SoilMoist 
level      : 1

`

为了实现这一点，您需要确保每个栅格的范围和分辨率相同，下面是我的代码，它生成示例栅格并通过它们进行处理以使它们具有相同的范围和分辨率。这是带注释的代码：

#### ----- Making Sample rasters ----- ####
e1 <- extent(-1.110223e-16, 360, -90, 90)
d1 <- raster(e1)
res(d1) <- 1
d1[] <- runif(64800, min = 0, max = 1)
#plot(d1)

e2 <- extent(48, 86, 6, 51)
a <- raster(e2)
res(a) <- c(0.826087, 1.153846)
a[] <- rnorm(1794, 5, 1)
#plot(a)

e3 <- extent(48, 86, 6, 51)
b <- raster(e3)
res(b) <- c(0.826087, 1.153846)
b[] <- rnorm(1794, 12, 1)

e3 <- extent(47.5, 86.5, 5.5, 51.5)
c <- raster(e3)
res(c) <- 1
c[] <- rnorm(1794, 12, 1)

crs(d1) = crs("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
crs(a) = crs("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
crs(b) = crs("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
crs(c) = crs("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
#### ----- Making Sample rasters ----- ####


##### Cropping all rasters to the lowest extent #####
ex <- intersect(intersect(intersect(extent(a), extent(b)), extent(c)), extent(d1))

d1_crop = crop(d1, ex)
a_crop = crop(a, ex)
b_crop = crop(b, ex)
c_crop = crop(c, ex)
##### Cropping all rasters to the lowest extent #####


##### Reprojecting the rasters to make them the same resolution, making them the same as d1's resolution #####
a_res = projectRaster(a_crop, d1_crop)
b_res = projectRaster(b_crop, d1_crop)
c_res = projectRaster(c_crop, d1_crop)
##### Reprojecting the rasters to make them the same resolution, making them the same as d1's resolution #####

##### Doing your calculations after the resolution and extents are the same.
e = d1_crop-a_res-b_res-c_res