如何在R中计算随机网格单元内BiomeRealm属性的面积（km²）

嘿，我来帮你搞定这13个随机网格单元的生态区面积统计！先梳理下现有代码里可以优化的地方，再一步步实现你要的功能：

第一步：修正网格单元格大小的定义 #

你原来用英尺转米的方式有点绕，因为咱们已经把坐标系转成了UTM（单位是米），所以直接把单元格大小设为c(2000, 2000)就代表2km×2km的网格，不用额外转英尺啦。

第二步：把随机抽取的点转换成网格多边形 #

你现在得到的s是随机点，但咱们需要的是这些点对应的完整网格单元多边形，这样才能和WWF生态区做空间交集计算面积。可以用SpatialGrid的属性来生成每个点对应的网格多边形。

第三步：计算每个网格内各生态区的面积 #

对每个网格多边形，用intersect函数和WWF生态区图层做交集，然后提取BiomeRealm属性并计算面积（注意把平方米转换成平方公里，除以1e6）。

完整代码实现 #

library(raster)
library(rgdal)
library(sp) # 补充sp包，处理空间对象更方便

# 加载行政边界并转换为UTM坐标系
shp <- getData(country = "HND", level = 0)
shp <- spTransform(shp, CRSobj = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

# 创建2km×2km的网格（直接用米为单位）
cs <- c(2000, 2000) # 2000米 = 2公里
grdpts <- makegrid(shp, cellsize = cs)
spgrd <- SpatialPoints(grdpts, proj4string = CRS(proj4string(shp)))
spgrdWithin <- SpatialPixels(spgrd[shp,])
spgrdWithin <- as(spgrdWithin, "SpatialGrid")

# 选取13个随机网格单元的中心点
s <- spsample(spgrdWithin, 13, "random")

# 把随机点转换为对应的网格多边形
# 获取网格的单元格大小和原点信息
grid_params <- spgrdWithin@grid
cell_width <- grid_params@cellsize[1]
cell_height <- grid_params@cellsize[2]
origin_x <- grid_params@origin[1]
origin_y <- grid_params@origin[2]

# 为每个随机点生成网格多边形
grid_polygons <- lapply(1:length(s), function(i) {
  # 计算当前点对应的网格边界
  x <- coordinates(s)[i,1]
  y <- coordinates(s)[i,2]
  xmin <- x - cell_width/2
  xmax <- x + cell_width/2
  ymin <- y - cell_height/2
  ymax <- y + cell_height/2
  # 创建多边形
  Polygon(cbind(c(xmin, xmax, xmax, xmin, xmin), c(ymin, ymin, ymax, ymax, ymin)))
})
grid_polygons <- SpatialPolygons(lapply(1:length(grid_polygons), function(i) {
  Polygons(list(grid_polygons[[i]]), ID = as.character(i))
}), proj4string = CRS(proj4string(shp)))

# 读取并转换WWF生态区坐标系
HNDWWF <- readOGR('.','HONDURAS_WWF')
HNDWWFTRANS <- spTransform(HNDWWF, CRSobj = proj4string(shp))

# 计算每个网格内各BiomeRealm的面积并输出
for (i in 1:length(grid_polygons)) {
  cat(paste0("网格单元", i, ": "))
  # 交集计算
  intersect_poly <- intersect(grid_polygons[i,], HNDWWFTRANS)
  if (!is.null(intersect_poly)) {
    # 计算面积并转换为平方公里
    intersect_poly$area_km2 <- gArea(intersect_poly, byid = TRUE)/1e6
    # 按BiomeRealm汇总面积
    area_summary <- tapply(intersect_poly$area_km2, intersect_poly$BiomeRealm, sum)
    # 格式化输出
    output_str <- paste(names(area_summary), "- 面积 -", round(area_summary, 2), "km²", collapse = " ")
    cat(output_str, "\n")
  } else {
    cat("无匹配的生态区\n")
  }
}

代码说明 #

• 我们先把随机点转换成了完整的网格多边形，这样才能准确计算和生态区的重叠面积
• 用gArea计算交集面积，除以1e6把平方米转换成平方公里
• 用tapply按BiomeRealm汇总每个网格内的生态区面积
• 最后按你想要的格式输出每个网格的统计结果

运行这段代码后，你会得到类似这样的输出：

网格单元1: NT14 - 面积 - 3.2 km² NT7 - 面积 - 0.8 km²
网格单元2: NT7 - 面积 - 4.0 km²
...

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