多边形并未按照 R 中使用 sf st_wrap_dateline 的预期那样在国际日期变更线上分割

Question

我有跨越日期变更线的多边形，但是相反，它们是反转并环绕远离数据变更线的。

作为示例，此处的代码应生成包含 122 个单元格的完整六角形网格，并在日期变更线上断开。我的想法是 st_wrap_dateline 应该通过沿着日期线插入一条线并创建多个多边形来导致倒置的多边形分裂。我不确定接下来要尝试什么。

代码和输出如下。

任何帮助将不胜感激。

install.packages("h3jsr")
install.packages("sf")

library("h3jsr")
library("sf")


# Make number sequence
pts  = as.data.frame(merge(seq (-90,90,1),seq (0,360,1),all=TRUE))
colnames(pts) <- c("LAT","LON")

pts_ <- st_as_sf(x = pts, coords = c("LON","LAT"),crs = 4326)

# resolve H3 index at Resolution
pts_H3 <- point_to_h3(pts_, res = 0, simple = FALSE)

pts_H3 <- as.data.frame(table(pts_H3$h3_resolution_0))

getH3_Poly <- h3_to_polygon(pts_H3$Var1,simple=FALSE)

# timeline wrap
getH3_Poly <- st_wrap_dateline(getH3_Poly, options = c("WRAPDATELINE=YES"))
plot(getH3_Poly)

Answer 1

只是探索一些关于创建较小示例的想法，获得更好的 bbox，但仍然没有解决 st_wrap_dateline 问题：

pts2 <- data.frame(merge(seq(-90,90,10),seq(0,360,10), all = TRUE))
colnames(pts2) <- c('LAT', 'LON')
pts_ <- st_as_sf(x = pts2, coords = c('LON', 'LAT'), crs=4326)
pts_H3 <- point_to_h3(pts_, res = 0, simple = FALSE)
pts_H3_poly <- h3_to_polygon(pts_H3$h3_resolution_0, simple = FALSE)
some_idx <- sapply(unique(pts_H3_poly$h3_address), function(i) which(pts_H3_poly$h3_address %in% i))
pts_poly_122 <- pts_H3_poly[unlist(unname(sapply(some_idx, `[[`, 1))), ]
pts_poly_122
Simple feature collection with 122 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -179.6744 ymin: -87.3647 xmax: 179.2172 ymax: 87.3647
Geodetic CRS:  WGS 84
First 10 features:
        h3_address h3_resolution                       geometry
1  80f3fffffffffff             0 POLYGON ((-179.6744 -73.310...
2  80effffffffffff             0 POLYGON ((-34.4418 -87.3647...
3  80e7fffffffffff             0 POLYGON ((48.29116 -69.3713...
4  80ddfffffffffff             0 POLYGON ((-7.138629 -66.192...
5  80d1fffffffffff             0 POLYGON ((4.435671 -35.7412...
6  80c1fffffffffff             0 POLYGON ((-15.54273 -23.004...
8  8099fffffffffff             0 POLYGON ((-11.66475 -4.4670...
10 8075fffffffffff             0 POLYGON ((-4.013998 11.5453...
11 8059fffffffffff             0 POLYGON ((-2.483865 22.1975...
13 8039fffffffffff             0 POLYGON ((-12.38413 40.8691...

诚然，这并不是那么有用，但可以更快地绘图。而 bbox 更接近预期。 然后尝试一下：

pts_poly_122_wrap2 <- st_wrap_dateline(pts_poly_122, options= c('WRAPDATELINE=YES', 'DATELINEOFFSET=20'))
plot(pts_poly_122_wrap2)

“DATELINEOFFSET=20”的选择只是为了关闭/与默认值不同，这是否更接近您的预期？

Answer 2

创建了一个使用

hextess()

包中的

statstat

的方法。这里的关键是逻辑

trim

选项，这是默认选项。边界上剪切的图块的大小/均匀度显然取决于数据的范围，但这仍然是对

st_wrap_dateline()

噩梦的改进。我确信使用输入参数可以减少这个问题。

请注意，

hextess()

比

st_make_grid()

慢得多，尤其是在更精细的分辨率下，但同样，至少这是有效的。

library(spatstat)
library(sf)
library(dplyr)
library(ggplot2)
library(rnaturalearth)

# Sample data to create tesselation for
countries <- rnaturalearth::ne_countries(returnclass='sf')

# Get extent of data to be 'gridded'
st_bbox(countries)
#       xmin       ymin       xmax       ymax 
# -180.00000  -90.00000  180.00000   83.64513 

# Create 'flat' hexagonal tesselation based on extent of data
temp1 <- hextess(owin(c(-180, 180), c(-90, 90)), # xmin xmax ymin ymax
                 10, # Length of hexagon side
                 trim = TRUE) # Trim to extent, this is the default

# Convert tess object to df
temp2 <- data.frame(temp1)

# Convert df to sf
hex_tess <- temp2 %>%
  st_as_sf(coords = c("x", "y"), crs = 4326) %>%
  group_by(Tile) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")

# Check for invalid geometry
any(st_is_valid(hex_tess) == FALSE)
# [1] FALSE

ggplot() +
  geom_sf(data = countries) +
  geom_sf(data = hex_tess,
          colour = "red",
          fill = NA)

结果：