如何基于栅格单元内长度占比对矢量线进行栅格化加权求和？

我有一个矢量线数据集，想要栅格化时对每个栅格单元内所有线的目标字段数值求和。但部分线会跨越多个栅格单元，为避免重复计算，需要按线在各栅格单元内的长度占比进行加权求和——比如某条线字段值为5，若平分在栅格单元A和B中，则A、B的求和结果分别增加2.5。发现R的terra::rasterize()无法实现这个逻辑，求R或Python中的现成实现或自定义方法。

补充可复现代码（当前直接累加触碰单元的线值）：

library(terra)
library(tidyterra)

f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
l <- as.lines(v)

# 给线要素赋值
l$values <- 1
l <- l[, 'values']
l <- project(l, '+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs')

# 栅格化（直接累加）
dummy_raster <- rast(l, res = 20000)
l_rast <- rasterize(l, dummy_raster, field = 'values', fun = 'sum')

# 可视化
ggplot() +
    geom_spatraster(data = l_rast) +
    scale_fill_hypso_b() +
    geom_spatvector(data = l, aes(colour = as.factor(sample(seq_along(l)))), show.legend = F) +
    geom_spatvector_label(data = l, mapping = aes(label = values)) +
    NULL

R实现方法 #

核心思路是将线要素切割到对应栅格单元，计算每段线的长度占原线的比例后加权，再进行栅格化求和，步骤如下：

1. 将目标栅格转换为矢量多边形（每个栅格单元对应一个面）
2. 用terra::intersect()切割线要素，得到每个栅格单元内的线段
3. 计算原线总长度、当前线段长度，得出长度占比
4. 计算加权值（原字段值 × 长度占比）
5. 对加权值执行栅格化求和

代码示例：

library(terra)
library(tidyterra)

# 加载并预处理数据（同原示例）
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
l <- as.lines(v)
l$values <- 1
l <- l[, 'values']
l <- project(l, '+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs +type=crs')

# 创建目标栅格
dummy_raster <- rast(l, res = 20000)

# 栅格转面要素（每个单元为一个多边形）
raster_polys <- as.polygons(dummy_raster, dissolve = FALSE)

# 切割线要素，得到栅格内的线段
cut_lines <- intersect(l, raster_polys)

# 计算长度占比与加权值
cut_lines$total_length <- length(l)[cut_lines$ID]
cut_lines$segment_length <- length(cut_lines)
cut_lines$weighted_value <- cut_lines$values * (cut_lines$segment_length / cut_lines$total_length)

# 栅格化加权值并求和
weighted_rast <- rasterize(cut_lines, dummy_raster, field = 'weighted_value', fun = 'sum')

# 可视化结果
ggplot() +
  geom_spatraster(data = weighted_rast) +
  scale_fill_hypso_b() +
  geom_spatvector(data = cut_lines, aes(colour = as.factor(ID)), show.legend = F) +
  geom_spatvector_label(data = cut_lines, mapping = aes(label = round(weighted_value, 2))) +
  NULL

Python实现方法 #

用geopandas结合rasterio实现，逻辑与R一致：切割线段、计算长度占比、加权后栅格化。

代码示例：

import geopandas as gpd
import rasterio
from rasterio.features import rasterize
from rasterio.transform import from_bounds
import numpy as np

# 加载数据并预处理
lux_gdf = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
lux_gdf = lux_gdf[lux_gdf.name == 'Luxembourg']
lines_gdf = lux_gdf.boundary.to_frame('geometry')
lines_gdf['values'] = 1
# 投影到指定坐标系
crs = '+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'
lines_gdf = lines_gdf.to_crs(crs)

# 设置栅格参数
res = 20000
xmin, ymin, xmax, ymax = lines_gdf.total_bounds
width = int((xmax - xmin) / res)
height = int((ymax - ymin) / res)
transform = from_bounds(xmin, ymin, xmax, ymax, width, height)

# 生成栅格单元多边形
grid_cells = []
for x in range(width):
    for y in range(height):
        x_left = xmin + x * res
        x_right = x_left + res
        y_bottom = ymin + y * res
        y_top = y_bottom + res
        grid_cells.append({'geometry': gpd.GeoSeries.from_wkt(f'POLYGON(({x_left} {y_bottom}, {x_right} {y_bottom}, {x_right} {y_top}, {x_left} {y_top}, {x_left} {y_bottom}))').iloc[0]})
grid_gdf = gpd.GeoDataFrame(grid_cells, crs=crs)

# 切割线段
cut_lines = gpd.overlay(lines_gdf, grid_gdf, how='intersection')

# 计算长度占比与加权值
cut_lines['total_length'] = lines_gdf.geometry.length[cut_lines.index.get_level_values(0)].values
cut_lines['segment_length'] = cut_lines.geometry.length
cut_lines['weighted_value'] = cut_lines['values'] * (cut_lines['segment_length'] / cut_lines['total_length'])

# 栅格化加权值
shapes = [(row.geometry, row.weighted_value) for _, row in cut_lines.iterrows()]
weighted_raster = rasterize(
    shapes,
    out_shape=(height, width),
    transform=transform,
    fill=0,
    dtype=np.float64,
    merge_alg=rasterio.enums.MergeAlg.add
)

# 保存结果
with rasterio.open(
    'weighted_raster.tif',
    'w',
    driver='GTiff',
    height=height,
    width=width,
    count=1,
    dtype=np.float64,
    crs=crs,
    transform=transform,
) as dst:
    dst.write(weighted_raster, 1)

内容的提问来源于stack exchange，提问作者tlhenvironment