PAYATHON 2026

如何在 R 中将 WGS84 空间数据转换为等积投影并完整显示研究区域？

数据小陈

可以让研究区域完整显示，同时保留栅格单元。不过，更换投影后无法同时保持原来的形状、距离和面积。如果主要目标是“等积，而且布局接近第一张经纬度地图”，可以先尝试圆柱等积投影 +proj=cea：将水深数据保留为 SpatRaster，用 terra::project() 重投影后再绘制。

原代码中还有两个前提需要纠正：coord_sf(crs = 4326) 使用的是 WGS84 地理坐标系，不是 Mercator；Aitoff 也不是等积投影。容易与它混淆的 Hammer 投影，也称 Hammer-Aitoff，才是等积投影。Aitoff 官方说明、Hammer 官方说明

先统一 CRS 定义

描述中的 ESRI:54043 与代码中的 53043 不一致。sf 接收纯数字 CRS 时会按 EPSG 编号解释，因此不能用裸数字代替 ESRI 编号。如果确实需要描述中的 Aitoff CRS，应明确写出命名空间：

aitoff_crs <- sf::st_crs("ESRI:54043")

# 仅在确实需要 Aitoff 时使用；它不满足等积要求
bathy_aitoff <- sf::st_as_sf(
  bathy,
  coords = c("V1", "V2"),
  crs = 4326
) |>
  sf::st_transform(aitoff_crs)

st_transform() 会转换坐标。st_set_crs() 或直接修改 CRS 属性只会重新标记坐标，不能完成重投影。sf CRS 文档

选择接近原图布局的等积投影

研究区域横跨 180° 经度，很难找到一种投影，让整个区域都保持接近原来的形状。圆柱等积投影的经纬网仍是矩形布局，可以先尝试这种投影。它保面积，但形状和纵横比例仍会变化。

下面以中央经线 70°、标准纬线 -40° 作为示例起点，其中 70° 是范围中点。这些参数并非经过优化得出的唯一选择。lat_ts 控制真比例纬线，lon_0 控制中央经线。PROJ 圆柱等积投影文档

target_crs <- sf::st_crs(
  "+proj=cea +lat_ts=-40 +lon_0=70 +datum=WGS84 +units=m +no_defs +type=crs"
)

如果希望改用扇形布局，可以评估 Albers Equal Area；如果绘制整个南极周边，可以考虑以南极为中心的等积方位投影。这两种投影都不会保留第一张图的矩形经纬网外观。

水深数据保留为栅格后再重投影

st_as_sf(..., coords = c("V1", "V2")) 会把每行坐标转换成 POINT。转换后不再有单元边界信息，geom_sf() 只能画点。

以下代码假定 bathy 是原示例中 marmap::as.xyz() 的结果，V1、V2 是规则网格的经纬度中心坐标，V3 是水深；同时假定 World 和 fronts 都带有正确的源 CRS。

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

# 从完整的规则 XYZ 数据恢复栅格
r <- terra::rast(
  as.data.frame(bathy[, c("V1", "V2", "V3")]),
  type = "xyz",
  crs = "EPSG:4326"
)
names(r) <- "depth"

# 与原图一致，只显示负值水深
r[r >= 0] <- NA

# 20000 表示示例输出分辨率为 20 km，需按源数据和制图需求调整
r_proj <- terra::project(
  r,
  target_crs$wkt,
  method = "bilinear",
  res = 20000
)

plot_data <- as.data.frame(r_proj, xy = TRUE, na.rm = TRUE)
e <- as.vector(terra::ext(r_proj))

p <- ggplot() +
  geom_raster(
    data = plot_data,
    aes(x = x, y = y, fill = depth)
  ) +
  geom_sf(data = World) +
  geom_sf(
    data = fronts,
    aes(color = factor(name)),
    linetype = 2
  ) +
  coord_sf(
    crs = target_crs,
    default_crs = target_crs,
    xlim = e[1:2],
    ylim = e[3:4],
    expand = FALSE
  ) +
  labs(
    fill = "Bathymetry (m)",
    color = "ACC fronts (Park & Durand, 2019)"
  ) +
  theme_bw() +
  theme(axis.title = element_blank())

p

这里先将水深重采样到目标投影下的规则网格，再用 geom_raster() 绘制。连续水深可以使用 bilinear，分类栅格通常应使用 near。重投影会估算新的单元值，输出分辨率也需要明确选择。更细的输出网格不意味着增加了真实数据细节。terra 重投影文档

如果需要叠加多个栅格或进行逐单元计算，应共用一个输出模板，使 CRS、范围、分辨率和网格原点一致。例如，可以将 terra::rast(r_proj) 作为后续 terra::project() 的第二个参数。

如何避免研究区域显示不全

非线性投影会把经纬度矩形的边界变成曲线。如果只转换几个角点，得到的坐标范围可能漏掉曲线上的极值。地图看起来偏向一侧，也可能与投影中央经线的位置有关，不能仅凭外观认定转换错误。

对于仅包含 sf 图层的地图，可以继续用经纬度指定范围，让 coord_sf() 计算能容纳整个经纬度框的投影范围：

ggplot() +
  geom_sf(data = World) +
  geom_sf(data = fronts, aes(color = factor(name))) +
  coord_sf(
    crs = target_crs,
    default_crs = sf::st_crs(4326),
    xlim = c(-20, 160),
    ylim = c(-60, -20),
    lims_method = "box",
    expand = FALSE
  )

lims_method = "box" 会选择能容纳整个指定框的显示范围，可能增加空白。经纬度框投影后的形状未必是矩形，因此这个显示范围不等同于严格按研究区边界裁剪。coord_sf() 文档

两段示例中的 default_crs 设置不同：栅格示例中的普通 x、y 已是投影坐标，因此使用 target_crs；上面的矢量示例用经纬度设置范围，因此使用 4326。不能把已投影的 geom_raster() 坐标再按经纬度解释。纬度范围也建议统一写成从小到大的 c(-60, -20)。

必须使用 sf 单元面时怎样转换

如果确实需要将每个栅格单元变成独立面，可以从重投影前的栅格生成多边形：

cells <- terra::as.polygons(
  r,
  aggregate = FALSE,
  values = TRUE,
  na.rm = TRUE
)

cells_sf <- sf::st_as_sf(cells)

# 用此图层替换前面地图中的水深 geom_raster() 图层
cell_layer <- geom_sf(
  data = cells_sf,
  aes(fill = depth),
  color = NA
)

aggregate = FALSE 防止相同数值的相邻单元被合并；fill 用于填充单元，color = NA 用于去掉边线。terra 单元面转换文档

这种方式保留源网格的单元划分，绘图时由 coord_sf(crs = target_crs) 转换。如果从 r_proj 转面，得到的则是重采样后的目标网格单元。将大量单元转换为面会明显增加内存和绘图开销，单纯显示水深时通常保留 SpatRaster 更合适。

备注：内容仅供参考。