r - How to subtract 3 rasters from 1 raster with different extents and resolutions

Question

I have 4 rasters with different resolutions and extents. Can anyone help me how to subtract 3 rasters (a,b,c) from raster (d1) to get the new output raster named as "e"

like 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

`

标签： rrstudio

Answer 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