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

【问题标题】：How to subtract 3 rasters from 1 raster with different extents and resolutions如何从具有不同范围和分辨率的 1 个栅格中减去 3 个栅格
【发布时间】：2019-12-13 22:19:11
【问题描述】：

我有 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

【问题讨论】：

标签： r rstudio

【解决方案1】：

为了实现这一点，您需要确保每个栅格的范围和分辨率都相同，这是我下面的代码，它制作了示例栅格并通过它们使它们具有相同的范围和分辨率。这是带有 cmets 的代码：

#### ----- 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

【讨论】：