WGS 84/EASE-Grid转Winkel-Tripel投影时几何图形无效问题求助

我编写了R代码，可在WGS 84 / NSIDC EASE-Grid 2.0 Global（EPSG:6933）投影下创建全球等面积六边形网格并绘制。但将网格重投影至Winkel-Tripel投影时，针对加拿大、俄罗斯等高纬度大国或全球范围会报错：

Evaluation error: IllegalArgumentException: Invalid number of points in LinearRing found 2 - must be 0 or >= 4.

尝试过多种底图、自动/手动生成的边界框，均无法解决；转换后部分几何图形失效，且st_make_valid()无法修复。

可运行示例（挪威区域） #

library(rnaturalearth)
library(sf)
library(dplyr)
library(ggplot2)
library(cowplot)

# 数据加载与转换
world = ne_countries(scale = 'medium', returnclass = 'sf') %>%
  st_transform(crs=4236)

map = world[world$admin == 'Norway', ]

# 生成边界框
bbox = st_bbox(map) %>%
  st_as_sfc() %>%
  st_sf()

# EPSG:6933下生成六边形网格
map_6933 = st_transform(map, 6933)
bbox_6933 = st_transform(bbox,6933)
bbox_grid_6933 <- st_make_grid(bbox_6933,
                               n = c(25,25),
                               what = 'polygons',
                               square = FALSE,
                               flat_topped = TRUE) %>%
  st_as_sf() %>%
  mutate(area = st_area(.)/c(1000*1000))

# 转换至Winkel-Tripel投影
wintri <- "+proj=wintri +datum=WGS84 +no_defs +over"
map_wintri = st_transform(map_6933, wintri)
bbox_wintri = st_transform(bbox_6933, wintri)
bbox_grid_wintri = st_transform(bbox_grid_6933,wintri)

# 绘图
p_6933 <- ggplot() + 
  geom_sf(data = bbox_grid_6933, aes(fill = units::drop_units(area)),alpha=0.75) +
  geom_sf(data = bbox_6933, fill = NA, color = 'white') +
  geom_sf(data = map_6933, fill = NA, color = 'white')

p_wintri <- ggplot() + 
  geom_sf(data = bbox_grid_wintri, aes(fill = units::drop_units(area)),alpha=0.75) +
  geom_sf(data = bbox_wintri, fill = NA, color = 'white') +
  geom_sf(data = map_wintri, fill = NA, color = 'white')

cowplot::plot_grid(p_6933, p_wintri, nrow = 2,labels = c('6933','wintri'))

示例1预期结果：挪威区域在EPSG:6933与Winkel-Tripel投影下的六边形网格对比图

报错示例 #

示例2：加拿大区域

将map = world[world$admin == 'Norway', ]改为map = world[world$admin == 'Canada', ]，运行报错：

Error in CPL_geos_is_empty(st_geometry(x)) :
Evaluation error: IllegalArgumentException: Invalid number of points in LinearRing found 2 - must be 0 or >= 4.

示例3：全球范围

将map = world[world$admin == 'Norway', ]改为map = world，运行报错：

Error in CPL_geos_op2(op, x, y) :
Evaluation error: IllegalArgumentException: point array must contain 0 or >1 elements.

根因 #

Winkel-Tripel是伪圆柱投影，在高纬度（尤其是接近南北极）区域的投影转换会导致几何严重变形：

• 部分六边形网格的顶点投影后重合，形成只有2个点的无效LinearRing（不符合GEOS要求的≥4个点的规则）
• 全球范围下，跨180°经线或覆盖两极的几何，投影后会出现拓扑断裂或空几何，触发point array must contain 0 or >1 elements错误

解决方案 #

针对高纬度国家或全球范围，需在投影前后增加几何校验与预处理步骤：

1. 预处理：过滤高纬度易失效网格（在EPSG:6933投影下）

在生成六边形网格后，先计算每个网格的中心坐标，过滤掉纬度超过±80°的网格（可根据需求调整阈值）：

# 在生成bbox_grid_6933后添加：
bbox_grid_6933 <- bbox_grid_6933 %>%
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(lat = st_coordinates(centroid)[,2]) %>%
  filter(abs(lat) < 80) %>% # 过滤高纬度网格
  select(-centroid, -lat)

2. 处理跨180°经线的几何

如果目标区域跨国际日期变更线（比如俄罗斯、全球），先在WGS84坐标系下用st_wrap_dateline拆分几何：

# 针对全球或跨180°的区域，在转换到6933前处理：
world <- ne_countries(scale = 'medium', returnclass = 'sf') %>%
  st_transform(crs=4326) %>% # 先回到WGS84
  st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) %>%
  st_transform(crs=4236)

3. 投影后清理无效几何

转换到Winkel-Tripel后，强制过滤无效几何：

# 转换网格后添加：
bbox_grid_wintri <- st_transform(bbox_grid_6933, wintri) %>%
  filter(st_is_valid(geometry)) # 过滤无效几何

4. 全球范围优化：拆分南北半球生成网格

全球范围直接生成网格容易触发投影问题，建议拆分南北半球分别生成再合并：

# 北半球bbox
bbox_north <- st_bbox(c(xmin=-180, ymin=0, xmax=180, ymax=80), crs=4326) %>%
  st_as_sfc() %>% st_sf() %>% st_transform(6933)
# 南半球bbox
bbox_south <- st_bbox(c(xmin=-180, ymin=-80, xmax=180, ymax=0), crs=4326) %>%
  st_as_sfc() %>% st_sf() %>% st_transform(6933)

# 分别生成网格
grid_north <- st_make_grid(bbox_north, n=c(50,25), square=FALSE, flat_topped=TRUE) %>% st_as_sf()
grid_south <- st_make_grid(bbox_south, n=c(50,25), square=FALSE, flat_topped=TRUE) %>% st_as_sf()

# 合并后投影
bbox_grid_6933 <- rbind(grid_north, grid_south) %>% mutate(area = st_area(.)/(1000*1000))

修改后的完整代码（以加拿大为例） #

library(rnaturalearth)
library(sf)
library(dplyr)
library(ggplot2)
library(cowplot)

# 数据预处理
world = ne_countries(scale = 'medium', returnclass = 'sf') %>%
  st_transform(crs=4326) %>%
  st_wrap_dateline(options = c("WRAPDATELINE=YES", "DATELINEOFFSET=180")) %>%
  st_transform(crs=4236)

map = world[world$admin == 'Canada', ]

# bounding box处理
bbox = st_bbox(map) %>%
  st_as_sfc() %>%
  st_sf()

# 转换到EASE-Grid 2.0
map_6933 = st_transform(map, 6933)
bbox_6933 = st_transform(bbox,6933)

# 生成六边形网格并过滤高纬度
bbox_grid_6933 <- st_make_grid(bbox_6933,
                               n = c(25,25),
                               what = 'polygons',
                               square = FALSE,
                               flat_topped = TRUE) %>%
  st_as_sf() %>%
  mutate(centroid = st_centroid(geometry)) %>%
  mutate(lat = st_coordinates(centroid)[,2]) %>%
  filter(abs(lat) < 80) %>%
  select(-centroid, -lat) %>%
  mutate(area = st_area(.)/c(1000*1000))

# 转换到Winkel-Tripel
wintri <- "+proj=wintri +datum=WGS84 +no_defs +over"
map_wintri = st_transform(map_6933, wintri)
bbox_wintri = st_transform(bbox_6933, wintri)
bbox_grid_wintri = st_transform(bbox_grid_6933,wintri) %>%
  filter(st_is_valid(geometry))

# 绘图
p_6933 <- ggplot() + 
  geom_sf(data = bbox_grid_6933, aes(fill = units::drop_units(area)), alpha=0.75) +
  geom_sf(data = bbox_6933, fill = NA, color = 'white') +
  geom_sf(data = map_6933, fill = NA, color = 'white')

p_wintri <- ggplot() + 
  geom_sf(data = bbox_grid_wintri, aes(fill = units::drop_units(area)), alpha=0.75) +
  geom_sf(data = bbox_wintri, fill = NA, color = 'white') +
  geom_sf(data = map_wintri, fill = NA, color = 'white')

cowplot::plot_grid(p_6933, p_wintri, nrow = 2,labels = c('6933','wintri'))

关键说明 #

• 高纬度阈值（80°）可根据实际需求调整，阈值越低，过滤的网格越多，投影成功率越高，但会丢失极区网格
• st_wrap_dateline必须在WGS84坐标系下使用，否则无法正确处理跨180°经线的几何
• 全球范围拆分南北半球是避免投影时出现两极拓扑问题的有效方法

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