从 WGS 84 / NSIDC EASE-Grid 2.0 Global sf_transform 到 Winkel-Tripel 后的无效几何

Question

我已经为全局等面积六角网格编写了代码。我可以使用 WGS 84 / NSIDC EASE-Grid 2.0 全局投影创建和绘制网格。我想将此网格重新投影到 Winkel-Tripel 投影上。该绘图适用于底图和边界框的某些组合，但不适用于其他组合。当它不起作用时，它会返回以下错误。

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

我尝试了很多地图组合、从地图派生的边界框和手动边界框，但似乎无法让它适用于整个世界或北纬大国，如美国、加拿大和俄罗斯。当它不起作用时，我发现某些几何图形在转换为 Winkel-Tripel 后变得无效，但 st_make_valid() 对它们不起作用。请参阅下面的有效代码和两个在尝试绘制 ggplot 对象时返回上述错误的代码示例。

感谢您提供的任何帮助！

格伦

#####################

例一

有效的代码

library(rnaturalearth)
library(sf)
library(dplyr)

################
#map
world = ne_countries(scale = 'medium', returnclass = 'sf') %>%
  st_transform(crs=4236)

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

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

# WGS 84 / NSIDC EASE-Grid 2.0 Global
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'))

预期结果...

#####################

示例 2

如果

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

更改为

，代码将不起作用

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

错误：

> cowplot::plot_grid(p_6933, p_wintri, nrow = 2,labels = c('6933','wintri'))
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

错误：

> cowplot::plot_grid(p_6933, p_wintri, nrow = 2,labels = c('6933','wintri'))
Error in CPL_geos_op2(op, x, y) : 
  Evaluation error: IllegalArgumentException: point array must contain 0 or >1 elements.