如何用 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()修复。
备注:内容仅供参考。
