如何用R语言raster包crop函数提取栅格数据至矢量网格？

嗨，我看到你在尝试用R的raster包处理栅格和矢量数据时遇到了问题——其实你用错函数啦！crop()的作用并不是把栅格属性和矢量网格关联起来，而是裁剪栅格的空间范围，得到的仍然是一个栅格图层，这就是为什么你导出的CSV只有shapefile的属性~ 下面我来帮你一步步解决这个问题，得到每个5×5km网格对应的土地利用信息：

核心问题解释 #

你想要的是将每个矢量多边形（5×5km网格）内的栅格土地利用数据进行统计关联，而crop()只是缩小栅格的空间范围，完全不会把栅格数据和矢量的属性表结合。我们需要用extract()（或更高效的替代工具）来完成这个关联任务。

前提确认 #

首先检查你的数据坐标系：你的shapefile（EPSG:3006）和土地利用栅格的坐标系完全匹配，不需要做投影转换，这一步已经没问题啦！

解决方案（分两种方案，推荐第二种） #

方案1：用raster包的extract()函数直接统计 #

这个方法直接针对每个矢量多边形提取栅格值并统计，适合数据量中等的场景：

library(raster)
library(rgdal)
library(sp)

# 1. 读取数据（和你原代码一致）
sweref.def <- "+init=epsg:3006"
ekrut <- readOGR(dsn = "//storage-al.slu.se/student$/nilc0001/Desktop/Nina/Ekrut", layer = "ekrut_5x5_flat", p4s = sweref.def)
landuse <- raster("nmd2018bas_ogeneraliserad_v1_0.tif")

# 2. 整理土地利用类型映射表（把ID转换成对应的类型名称）
landuse_attr <- levels(landuse)[[1]]
# 过滤掉无类型名称的行（比如ID=0、1）
landuse_attr <- landuse_attr[!is.na(landuse_attr$Klass), ]
# 建立ID到类型名称的映射
id_to_klass <- setNames(landuse_attr$Klass, landuse_attr$ID)

# 3. 提取每个网格内的土地利用类型统计（用table统计各类型数量）
# 注意：fun参数让我们在提取时直接统计，避免加载大量栅格值到内存
landuse_stats <- extract(
  x = landuse, 
  y = ekrut, 
  fun = function(x) table(factor(x, levels = names(id_to_klass))), 
  na.rm = TRUE
)

# 4. 转换统计结果为数据框，并重命名列
landuse_stats_df <- as.data.frame(do.call(rbind, landuse_stats))
colnames(landuse_stats_df) <- id_to_klass

# 5. 合并原shapefile的属性和土地利用统计结果
result_df <- cbind(ekrut@data, landuse_stats_df)

# 6. 导出最终CSV
write.csv(result_df, "landuse_ekrut_stats.csv", row.names = FALSE)

方案2：用terra+fasterize高效处理大数据 #

你的土地利用栅格非常大（110亿个栅格单元），raster包可能会遇到内存瓶颈。推荐用更现代的terra包（raster的替代工具，性能更强）结合fasterize来处理：

library(terra)
library(fasterize)

# 1. 用terra读取数据（比raster更高效）
ekrut_terra <- vect("//storage-al.slu.se/student$/nilc0001/Desktop/Nina/Ekrut/ekrut_5x5_flat.shp")
landuse_terra <- rast("nmd2018bas_ogeneraliserad_v1_0.tif")

# 2. 整理土地利用类型映射表
landuse_attr <- cats(landuse_terra)[[1]]
landuse_attr <- landuse_attr[!is.na(landuse_attr$Klass), ]
id_to_klass <- setNames(landuse_attr$Klass, landuse_attr$ID)

# 3. 将矢量网格转成与栅格同分辨率的栅格图层（用唯一标识BK_flat作为值）
ekrut_raster <- fasterize(ekrut_terra, landuse_terra, field = "BK_flat")

# 4. 用zonal函数统计每个网格的土地利用类型
landuse_zonal <- zonal(
  x = landuse_terra, 
  z = ekrut_raster, 
  fun = function(x) table(factor(x, levels = names(id_to_klass))), 
  na.rm = TRUE
)

# 5. 合并原属性和统计结果
result_df <- merge(as.data.frame(ekrut_terra), landuse_zonal, by.x = "BK_flat", by.y = "zone")

# 6. 导出CSV
write.csv(result_df, "landuse_ekrut_zonal.csv", row.names = FALSE)

结果说明 #

导出的CSV会包含原shapefile的所有属性（AREA、PERIMETER、BK等），再加上每列对应一种土地利用类型的计数（即该网格内该类型的栅格数量，乘以100就是对应的面积，因为栅格分辨率是10m×10m=100m²）。如果需要占比，可以自己计算（比如某类型计数/总计数）。

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