如何在R中生成空间对象面积的纬度分布密度图或直方图？

下面给你两种靠谱的实现方法，都能准确统计每个纬度带的面积，解决你之前点采样只覆盖边界的问题：

方法一：用纬度条带空间相交计算 #

直接创建覆盖目标区域的水平纬度条带，和你的SF对象做空间相交，就能得到每个条带内的实际区域，再计算面积就行，结果精准且覆盖所有内部区域。

代码示例：

library(sf)
library(dplyr)
library(ggplot2)
library(rnaturalearth)

# 获取巴西的SF对象
brazil <- ne_countries(continent = "South America", returnclass = "sf") %>%
  filter(name == "Brazil")

# 确定巴西的纬度范围，设置1度间隔的纬度断点
lat_range <- st_bbox(brazil)[c("ymin", "ymax")]
lat_steps <- seq(lat_range[1], lat_range[2], by = 1)

# 生成每个纬度区间的矩形条带（覆盖全球经度范围）
lat_bands <- lapply(seq_along(lat_steps)[-length(lat_steps)], function(i) {
  st_polygon(list(rbind(
    c(-180, lat_steps[i]),
    c(180, lat_steps[i]),
    c(180, lat_steps[i+1]),
    c(-180, lat_steps[i+1]),
    c(-180, lat_steps[i])
  ))) %>%
    st_sfc(crs = st_crs(brazil)) %>%
    st_sf(id = i, lat_mid = (lat_steps[i] + lat_steps[i+1])/2)
}) %>%
  bind_rows()

# 计算每个条带与巴西的交集面积（转成平方公里）
brazil_lat_area <- st_intersection(brazil, lat_bands) %>%
  mutate(area_sqkm = st_area(geometry)/1e6) %>%
  group_by(id, lat_mid) %>%
  summarise(total_area = sum(area_sqkm)) %>%
  st_drop_geometry()

# 绘制柱状图（纬度为Y轴，面积为X轴）
ggplot(brazil_lat_area, aes(x = total_area, y = lat_mid)) +
  geom_col(width = 0.8) +
  labs(x = "面积（平方公里）", y = "纬度（度）", title = "巴西各纬度带的面积分布") +
  theme_minimal()

方法二：栅格化统计 #

如果处理大区域追求效率，可以把SF对象转成栅格，统计每个纬度行的有效单元格数量，再乘以单元格面积得到对应纬度的总面积，分辨率越高结果越准。

代码示例：

library(raster)

# 定义栅格分辨率（这里用0.5度，调小数值能提升精度）
res <- 0.5
# 创建覆盖巴西的空栅格
brazil_raster <- raster(extent(brazil), res = res, crs = st_crs(brazil)$wkt)
# 将巴西的SF对象栅格化（有效区域标记为1）
brazil_rasterized <- rasterize(brazil, brazil_raster, field = 1)

# 统计每个纬度行的单元格数，计算对应面积
lat_area_raster <- tibble(
  lat = yFromRow(brazil_rasterized),
  cell_count = rowSums(getValues(brazil_rasterized), na.rm = TRUE)
) %>%
  mutate(
    # 计算单个单元格的平方公里面积（考虑纬度对经度距离的影响）
    cell_area = res*res*(111.32^2)*cos(lat*pi/180),
    total_area = cell_count * cell_area
  )

# 绘制折线密度图
ggplot(lat_area_raster, aes(x = total_area, y = lat)) +
  geom_line(size = 1) +
  labs(x = "面积（平方公里）", y = "纬度（度）", title = "巴西各纬度的面积分布（栅格法）") +
  theme_minimal()

补充说明 #

你之前用st_cast()转多点的方法只能提取多边形的边界点，完全没法覆盖内部区域，统计出来的结果肯定不准。上面两种方法都能完整覆盖所有区域，结果更可靠。

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