无法从单个 y 坐标创建光栅几何

Question

我使用如下所示的多边形裁剪了一个栅格。

我想调整 2 个栅格像元的值以说明它们对多边形的面积贡献。例如如果单元格 1 有 90% 的面积在多边形内，我想将它的值乘以 0.9，另一个单元格也一样。

我采用的方法是使用

terra::extract

和

xy = T

和

weights = T

来提取单元格值及其权重

library(terra)
library(dplyr)
  
crop_rast_df <- terra::extract(crop_rast, my_shp , weights = T, touches = T, xy = T)
crop_rast_df
  
ID GDP_PPP_3       x      y   weight
1  79067856 -85.875 32.125 0.874914
1 344945408 -85.625 32.125 0.932744
  
# adjust the value
crop_rast_df <- crop_rast_df %>% 
                  dplyr::mutate(eff_gdp = GDP_PPP_3 * weight) %>% # account for cell contribution
                  dplyr::select(x, y, eff_gdp)
crop_rast_df
x      y   eff_gdp
-85.875 32.125  38269172
-85.625 32.125 177990394
  

但是，当我想再次将其转回光栅时，出现以下错误：

crop_rast_adj <- terra::rast(crop_rast_df, type = 'xyz', crs = crs(temp_gdp))
Error: [rast] cannot create a raster geometry from a single y coordinate

我该如何纠正？

Answer 1

错误重现

library(terra)
library(dplyr)

crop_rast <- rast(matrix(1:2, nrow = 1, ncol = 2))
names(crop_rast) <- "GDP_PPP_3"
my_shp <- vect("POLYGON ((.5 0, .5 1, 1.5 0, .5 0))")

plot(crop_rast)
lines(my_shp)

crop_rast_df <- terra::extract(crop_rast, my_shp , weights = T, touches = T, xy = T)
crop_rast_df

  ID GDP_PPP_3   x   y weight
1  1         1 0.5 0.5  0.375
2  1         2 1.5 0.5  0.125

crop_rast_df <- crop_rast_df %>% 
  dplyr::mutate(eff_gdp = GDP_PPP_3 * weight) %>% # account for cell contribution
  dplyr::select(x, y, eff_gdp)
crop_rast_df

    x   y eff_gdp
1 0.5 0.5   0.375
2 1.5 0.5   0.250

crop_rast_adj <- terra::rast(crop_rast_df, type = 'xyz', crs = crs(temp_gdp))
Error: [rast] cannot create a raster geometry from a single y coordinate

建议的解决方案

我认为没有办法

terra::rast

从具有单个 y（或 x）坐标的数据集创建栅格。
因此，我宁愿建议复制您的初始栅格并使用调整后的栅格更改其值，如下所示：

crop_rast_adj <- crop_rast
values(crop_rast_adj) <- crop_rast_df$eff_gdp
crop_rast_adj

class       : SpatRaster 
dimensions  : 1, 2, 1  (nrow, ncol, nlyr)
resolution  : 1, 1  (x, y)
extent      : 0, 2, 0, 1  (xmin, xmax, ymin, ymax)
coord. ref. :  
source(s)   : memory
name        : GDP_PPP_3 
min value   :     0.250 
max value   :     0.375 

plot(crop_rast_adj)
lines(my_shp)

Answer 2

L--的示例数据

library(terra)

r <- rast(matrix(1:2, nrow = 1, ncol = 2), crs="+proj=utm +zone=1")
v <- vect("POLYGON ((.5 0, .5 1, 1.5 0, .5 0))")

方法 1

e <- extract(r, v, exact=TRUE, cell=TRUE)
# if there were multiple polygons, you would need to sum the fractions
e <- aggregate(e[,"fraction", drop=FALSE], e[,"cell", drop=FALSE], sum)

x <- rast(r)
x[e$cell] <- e$fraction
x * r
#class       : SpatRaster 
#dimensions  : 1, 2, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 0, 2, 0, 1  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=utm +zone=1 +datum=WGS84 +units=m +no_defs 
#source(s)   : memory
#name        :     lyr.1 
#min value   : 0.2499978 
#max value   : 0.3750118 

方法 2

a <- rasterizeGeom(v, r, "area")
b <- cellSize(a, transform=FALSE) # use TRUE with real data
d <- a/b
d * r
#class       : SpatRaster 
#dimensions  : 1, 2, 1  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 0, 2, 0, 1  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=utm +zone=1 +datum=WGS84 +units=m +no_defs 
#source(s)   : memory
#name        :  area 
#min value   : 0.250 
#max value   : 0.375 

如果不是所有单元格都至少部分被多边形覆盖，这些方法也有效。