如何栅格化区域变量的矢量？

Question

我用过很多次相同的光栅化过程，效果相当好：

raster <- rasterize(vect(shapefile.shp), base_grid, "my_variable")

其中 raster 是栅格化的 shapefile，shapefile.shp 是原始矢量，base_grid 是栅格骨架，“我的变量”是要考虑的变量。对于与多边形面积无关的变量，这种方法是令人满意的，因为它使用均值计算来重新排列数据（例如：人口、生产、产量、温度、降水量）。

但是，现在我正在尝试将矢量多边形转换为具有可变收获面积的栅格多边形，严格来说，这不是多边形的面积，但可以认为它与总多边形面积成正比。上述方法在考虑收获区域时会产生夸大的值（相应矢量的 3-4 倍），可能是因为多边形通常大于网格单元。因此，一个有 100 个单位的大多边形被分成 10 个网格单元，每个单元将给出 100 个单元，而我希望它们每个有 10 个单元（因为它是一个区域）。

我想我的方法是这样工作的：“在每个网格单元中，根据它们在网格单元中的存在成比例地对所有多边形值进行加权平均”

但我正在寻找的是：“对于每个网格单元中的每个多边形部分，计算网格单元内多边形的分数（wrt 到总多边形面积）并对网格单元内的所有值求和（因为它是一个面积单位)"。

非常感谢任何帮助。

更新：

矢量数据视图。栅格其实是多方面的，因为我有很多年：

Simple feature collection with 9382 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XYZ
Bounding box:  xmin: -67.38379 ymin: -41.03791 xmax: -53.63737 ymax: -21.99877
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84
First 10 features:
       ADM2_REF anio my_variable                       geometry
2  Tres Arroyos 1978                     180 MULTIPOLYGON Z (((-60.16947...
3  Tres Arroyos 1979                       0 MULTIPOLYGON Z (((-60.16947...
4  Tres Arroyos 1988                    1000 MULTIPOLYGON Z (((-60.16947...
5  Tres Arroyos 1989                    1000 MULTIPOLYGON Z (((-60.16947...
6  Tres Arroyos 1990                    3000 MULTIPOLYGON Z (((-60.16947...
7  Tres Arroyos 1991                    1500 MULTIPOLYGON Z (((-60.16947...
8  Tres Arroyos 1992                    2800 MULTIPOLYGON Z (((-60.16947...
9  Tres Arroyos 1993                    2800 MULTIPOLYGON Z (((-60.16947...
10 Tres Arroyos 1994                    2500 MULTIPOLYGON Z (((-60.16947...
11 Tres Arroyos 1995                    1250 MULTIPOLYGON Z (((-60.16947...

将上面的数据帧转换为栅格的步骤是：

baserast <- rast(nrows=nrows, ncol=nrows,
                 extent= extent,
                 crs="+proj=longlat +datum=WGS84",
                 vals=NA)

rasters <- rast(lapply(1978:2019, 
                       function(x)
                         rasterize(vect(shp.soy.yld.arg %>% 
                                          filter(anio==x)), baserast, "my variable")))

Link 一年的数据 .gpkg（对于所有年份来说都太大了）。

Answer 1

如果我很好地理解你的问题（一个可重现的例子将不胜感激），你希望光栅化多边形中的所有像素总和为收获的值（代码中的“my_variable”）。

在这里，我创建了一个玩具示例来向您展示我的推理：

首先加载库
使用总面积和收获面积示例创建玩具数据
计算每个像素被多边形覆盖的比例

将每个覆盖率除以多边形的总面积，然后乘以采伐面积

    library(sf)
    library(raster)
    library(exactextractr)

    rast <- raster::raster(matrix(rep(1,100), ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
    pol  <- sf::st_sfc(sf::st_polygon(list(cbind(c(0.5,4,7,0.5),c(1,0,4,1)))))
    pol  <- st_sf(data.frame(area = st_area(pol),harvest=0.7, geom=pol))
    raster::plot(rast)
    raster::plot(pol,add=T)

    cov_frac <- exactextractr::coverage_fraction(rast, pol)[[1]]
    raster::plot(cov_frac)
    raster::plot(pol,add=T)
    result <- cov_frac/st_area(pol)*pol$harvest
    sum(values(result))

如您所见，光栅化多边形中所有像素的总和等于采伐面积

Answer 2

在这种情况下，您应该栅格化密度而不是面积

示例数据：

library(terra)
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
# all areas have 100 ha of the crop
v$crop_area <- 100
r <- rast(v, res=.01)

计算密度

e <- expanse(v, unit="ha") # or "km" to avoid small numbers
v$density <- v$crop_area / e

栅格化

x <- rasterize(v, r, "density")

返回区域

ra <- cellSize(r, unit="ha")         
area <- x * ra

检查数字是否合理（大面积/小小区误差应该最小）。每个多边形的预期值为 100。

extract(area, v, sum, na.rm=TRUE, exact=TRUE) |> round()
#      ID density
# [1,]  1      99
# [2,]  2     101
# [3,]  3      99
# [4,]  4      95
# [5,]  5      99
# [6,]  6      98
# [7,]  7      96
# [8,]  8      99
# [9,]  9      98
#[10,] 10      98
#[11,] 11     101
#[12,] 12     100

下面我将展示如何在多年的循环中执行此操作

示例数据：

library(terra)
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
# all areas have 100 ha of the crop
v$crop_area <- 100
v$year <- rep(c(1990,1991, 1992), each=4)
r <- rast(v, res=.01)

解决方案

ra <- cellSize(r, unit="ha")         
e <- expanse(v, unit="ha") 
v$density <- v$crop_area / e

years <- unique(v$year)
out <- list()
for (i in 1:length(years)) {
   vv <- v[v$year == years[i], ]
   x <- rasterize(vv, r, "density")
   out[[i]] <- x * ra
}
out <- rast(out)
names(out) <- years

out
#class       : SpatRaster 
#dimensions  : 73, 78, 3  (nrow, ncol, nlyr)
#resolution  : 0.01, 0.01  (x, y)
#extent      : 5.74414, 6.52414, 49.44781, 50.17781  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#sources     : memory  
#              memory  
#              memory  
#names       :      1990,      1991,      1992 
#min values  : 0.2544571, 0.3028134, 0.3200223 
#max values  : 1.0492719, 0.6249076, 0.4335355

如何栅格化区域变量的矢量？

How to rasterize a vector for area variable?

r

raster

r-raster

rasterizing

terra