R语言ggplot绘制大量GeoJSON多边形时遇栈溢出问题求解决方案

问题背景 #

需要在R中绘制GeoJSON文件里的大量多边形（约1000个及以上），原代码通过循环逐个添加geom_sf图层，绘制超过250个多边形时出现栈溢出错误：

Error: C stack usage 7970756 is too close to the limit

原代码如下：

library ("ggplot2")
library ("sf")
library ("rnaturalearth")
library ("rnaturalearthdata")
library ("ggspatial")
library ("geojsonR")

europe <- ne_countries (continent = "europe", scale = "medium", returnclass = "sf")
class (europe)

p <- ggplot (data = europe)
p <- p + geom_sf ()

js <- FROM_GeoJson (url_file_string = "file.geojson")

i <- 0
for (f in js$features) {
    if (!is.null (f$properties$InterferenceLevel)) {
        if (i < 250) {
            overlay_polygon <- f$geometry$coordinates
            overlay_polygon <- sf::st_polygon (list(f$geometry$coordinates))
            poly <- sf::st_as_sfc (list (overlay_polygon), crs = 4326)

            p <- p + geom_sf (data = poly, fill=alpha ("red", 0.2), color=alpha ("red", 0.2))
        }
        i <- i + 1
    }
}

p <- p + coord_sf(xlim = c(-20, 42), ylim = c(35, 72))
plot (p)

核心问题原因 #

原代码通过循环逐个向ggplot对象添加独立的geom_sf图层，每一个多边形对应一个图层。当多边形数量超过几百个时，累积的图层会耗尽栈内存，触发报错。

解决方案 #

方法1：sf+ggplot优化（推荐） #

直接将GeoJSON文件读取为sf对象，一次性绘制所有多边形，避免循环添加图层：

library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)

# 读取欧洲底图
europe <- ne_countries(continent = "europe", scale = "medium", returnclass = "sf")

# 直接读取GeoJSON为sf对象，自动整合所有多边形
overlay_sf <- st_read("file.geojson")
# 过滤出包含InterferenceLevel属性的多边形
overlay_sf <- overlay_sf[!is.na(overlay_sf$InterferenceLevel), ]

# 一次性绘制所有多边形，仅需一个图层
p <- ggplot(data = europe) +
  geom_sf() +
  geom_sf(data = overlay_sf, fill = alpha("red", 0.2), color = alpha("red", 0.2)) +
  coord_sf(xlim = c(-20, 42), ylim = c(35, 72))

plot(p)

方法2：sf原生base plot（轻量高效） #

无需依赖ggplot，直接用sf的原生绘图功能，适合处理大量空间数据：

library(sf)
library(rnaturalearth)
library(rnaturalearthdata)

# 读取底图和多边形数据
europe <- ne_countries(continent = "europe", scale = "medium", returnclass = "sf")
overlay_sf <- st_read("file.geojson")
overlay_sf <- overlay_sf[!is.na(overlay_sf$InterferenceLevel), ]

# 绘制底图
plot(st_geometry(europe), xlim = c(-20, 42), ylim = c(35, 72))
# 添加透明多边形，add=TRUE表示叠加在已有图上
plot(st_geometry(overlay_sf), 
     col = adjustcolor("red", alpha.f = 0.2), 
     border = adjustcolor("red", alpha.f = 0.2), 
     add = TRUE)

方法3：tmap包（专业地图工具） #

tmap是专为空间数据可视化设计的包，处理大量多边形效率更高：

library(tmap)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)

# 读取数据
europe <- ne_countries(continent = "europe", scale = "medium", returnclass = "sf")
overlay_sf <- st_read("file.geojson")
overlay_sf <- overlay_sf[!is.na(overlay_sf$InterferenceLevel), ]

# 绘制地图
tm_shape(europe) +
  tm_borders() +
  tm_shape(overlay_sf) +
  tm_fill(col = "red", alpha = 0.2) +
  tm_borders(col = "red", alpha = 0.2) +
  tm_xlim(-20, 42) +
  tm_ylim(35, 72)

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