PAYATHON 2026

如何用 geom_sf() 国家边界筛选或裁剪 geom_tile() 网格?

数据阿明

可以先将每个 tile 转成矩形多边形,再用 sf::st_intersects() 判断它与国家区域是否相交,并根据结果筛选原始数据。这样可以自动删除完全位于境外的网格,保留跨越国界的完整网格。

如果希望网格的显示范围严格止于国界,可以用 sf::st_intersection() 求交集,再用 geom_sf() 绘制裁剪后的多边形。

为什么叠加国界线不会裁剪网格?

geom_sf(data = country, fill = NA) 只会叠加国家轮廓,不会对其他图层执行空间筛选。调整坐标范围也只能限制显示窗口,无法沿不规则国界删除网格。

在 geom_tile() 中,long 和 lat 表示网格中心。要判断整个网格是否完全位于境外,需要使用网格的矩形范围。只检查中心点,会漏掉中心在境外、但部分面积落在境内的网格。

自动删除完全位于境外的网格

下面的代码适用于问题中的规则经纬度网格。先根据相邻中心的间距确定宽高,再为每一行构建对应的矩形。

library(giscoR)
library(ggplot2)
library(tidyr)
library(sf)

country <- gisco_get_countries(country = "ITA")
country <- st_transform(country, 4326)

df <- expand_grid(
  long = seq(6, 19, length.out = 10),
  lat = seq(35, 48, length.out = 10)
)

set.seed(42)
df$z <- sample(100)

# 当前示例中,网格宽高等于相邻中心的间距
# 每个方向至少需要两个不同坐标,并且间距应保持一致
 dx <- diff(sort(unique(df$long)))[1]
 dy <- diff(sort(unique(df$lat)))[1]

tile_geometry <- st_sfc(
  lapply(seq_len(nrow(df)), function(i) {
    x <- df$long[i]
    y <- df$lat[i]

    ring <- rbind(
      c(x - dx / 2, y - dy / 2),
      c(x + dx / 2, y - dy / 2),
      c(x + dx / 2, y + dy / 2),
      c(x - dx / 2, y + dy / 2),
      c(x - dx / 2, y - dy / 2)
    )

    st_polygon(list(ring))
  }),
  crs = 4326
)

tiles <- st_sf(df, geometry = tile_geometry)

# 对意大利示例,转到 EPSG:3857 后执行平面空间关系判断
# 此处不进行面积或距离计算
 tiles_m <- st_transform(tiles, 3857)
 country_m <- st_transform(country, 3857)

keep <- lengths(st_intersects(tiles_m, country_m)) > 0

ggplot() +
  geom_tile(
    data = df[keep, ],
    aes(x = long, y = lat, fill = z),
    width = dx,
    height = dy
  ) +
  geom_sf(data = country, color = "red", fill = NA) +
  coord_sf(
    crs = st_crs(4326),
    default_crs = st_crs(4326)
  )

expand_grid() 来自 tidyr,因此完整示例需要加载该包。

st_intersects() 默认返回一个列表,每个元素记录对应网格与哪些国家要素相交。lengths(...) > 0 表示至少存在一个相交对象,结果可以直接用于筛选 df。

“相交”包括接触边界,因此只与国界相切的网格也会保留。跨越国界的网格,其境外部分仍会显示,符合只删除完全位于境外的网格这一要求。

让网格严格沿国界裁剪

如果不希望显示跨界网格的境外部分,可以继续使用上面创建的 tiles_m 和 country_m:

# 合并国家几何,避免按多个国家要素分别处理
country_area <- st_union(country_m)

# 保留每个网格落在国家区域内的部分
clipped <- st_intersection(tiles_m, country_area)

ggplot() +
  geom_sf(
    data = clipped,
    aes(fill = z),
    color = NA
  ) +
  geom_sf(
    data = country_m,
    color = "red",
    fill = NA
  ) +
  coord_sf(crs = st_crs(4326))

裁剪后,边缘网格可能变成不规则多边形,所以这里使用 geom_sf()。原来的 z 会作为网格属性保留下来,仍可用于填色;求交集本身不会重新计算这些数值。

使用时注意哪些条件?

  • 上面的矩形构造方式依赖规则网格。如果实际绘图指定了不同的 width、height,构造多边形时也要使用相同尺寸。
  • 经度作为 x、纬度作为 y,并按 WGS 84 经纬度赋予 EPSG:4326。如果输入坐标来自其他坐标系,需要填写真实 CRS,再进行转换。
  • EPSG:3857 在这里用于意大利示例的相交和裁剪操作,不应用这段代码的结果直接进行精确面积统计。换成其他区域时,应根据区域和计算目标选择投影。
  • 如果空间运算报告几何无效,可以先用 st_is_valid() 检查,再按需要用 st_make_valid() 修复。

备注:内容仅供参考。