瑞士地形SHP文件绘图报错：无法分配超大内存向量求助

错误原因 #

你遇到的内存溢出问题，核心是用错了绘图方法：geom_raster要求输入是规则网格的栅格数据，但dhm25_p.shp是离散的高程点数据集。直接将点数据转成data.frame后用geom_raster，R会尝试为所有可能的x/y组合生成栅格单元，这必然导致内存爆炸（错误提示的6539542.5 Gb就是明证）。

解决方案 #

1. 直接可视化离散点（快速简单） #

如果仅需要展示高程点分布，用geom_point替代geom_raster，无需转换为Spatial对象：

library(sf)
library(ggplot2)

# 读取SHP点数据
aux <- st_read(dsn="dhm25_p.shp")

# 查看数据列名，找到存储高程的字段（比如Z、H或elevation）
names(aux)

ggplot(aux) +
  geom_point(aes(color = Z), size = 0.1) + # 替换Z为实际的高程列名
  scale_color_viridis_c(option = "plasma") +
  theme_minimal()

2. 插值为栅格后绘制连续地形 #

如果需要生成平滑的地形表面，用插值工具将离散点转换为规则栅格：

library(sf)
library(terra)
library(ggplot2)

# 读取点数据
aux <- st_read(dsn="dhm25_p.shp")
# 转换为terra支持的矢量格式
v <- vect(aux)
# 创建对应DHM25精度的空栅格模板（25米分辨率）
r <- rast(v, resolution = 25)
# 反距离加权插值生成高程栅格
r_elev <- interpolate(r, v, method = "idw")

# 转换为ggplot可用的数据格式
r_df <- as.data.frame(r_elev, xy = TRUE)

ggplot(r_df) +
  geom_raster(aes(x = x, y = y, fill = Z)) + # Z为插值后的高程列名
  scale_fill_viridis_c(option = "plasma") +
  coord_sf(crs = st_crs(aux)) +
  theme_minimal()

3. 内存优化技巧 #

• 裁剪数据范围：如果不需要全瑞士地形，先裁剪到目标区域再处理：

  # 定义裁剪范围（示例坐标，替换为你需要的区域）
  bbox <- st_bbox(c(xmin=600000, ymin=200000, xmax=650000, ymax=250000), crs=st_crs(aux))
  aux_cropped <- st_crop(aux, bbox)
  

• 降低栅格分辨率：内存不足时，可将resolution设为更大的值（比如50米）

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